# Evaluate approximation of PCA from randomized algorithms

I have been comparing different PCA implementations (some via explicit calculation of the covariance matrix, some with randomized/truncated SVD) in terms of speed, and now wanted to compare how good different randomized algorithms approximate the exact (to the numerical precision) solution via covariance matrix.

Are there any suggested metrics to evaluate the similarity of two resulting eigenvector bases, i.e. the $k$ eigenvectors returned?

My intuition would be to assert that the vectors have the same "general direction" (i.e. in case the first dimension does not match directions, one vector would be flipped), and then computing either

1. the Frobenius norm of the absolute value matrix $abs(A)=abs(X-V)$, where $X$ is the exact solution and $V$ the resulting matrix from a randomized algorithm containing the eigenvectors, or
2. A modified/weighted version of the Euclidean distance, on the pairwise vectors. I.e., $$\mathcal{L}(X,V)=\sum_{j=0}^k \omega_j \sum_{i=0}^n \sqrt{(a_{ji}-x_{ji})^2}$$ where $\omega_j$ is a weight for the $j$-th PC, $n$ the dimension of the eigenvectors, and $a_{ji}$ and $x_{ji}$ the elements of the eigenvector matrices, respectively.

The reason for the second version could be that I want to have a better approximation for the first few components (as they explain more of the variance; maybe even weigh them according to that explained variance), and not so much care about the last percent of remaining variance spread across multiple dimensions, potentially accumulating a high loss.

Generally with things like this you really care about the similarity of the subspaces rather than the particular basis vectors, since you could keep the subspace the same yet change the particular eigenvectors by quite a lot.

Here's an example of that: suppose your features are uncorrelated so you end up with a covariance matrix of $C :=\text{diag}(1,1,0)$. Taking the top two eigenvectors is reasonable here. But what are they? This is a two dimensional eigenspace so e.g. $$B_1 = \left(\begin{array}{cc}1 & 0 \\ 0 & 1 \\ 0 & 0\end{array}\right), B_2 = \frac{1}{\sqrt 2}\left(\begin{array}{cc}1 & 1 \\ -1 & 1 \\ 0 & 0\end{array}\right)$$ both provide orthonormal bases. So if you're just looking at the entries, you could compute $\|B_1-B_2\|_F = 2$, but that's unfair because really they're both exactly correct bases for the same space.

So that's the motivation for this distance: for two matrices $A$ and $B$ with orthonormal columns, consider $$d(A, B) := \|P_A - P_B\|$$ where $\|\cdot\|$ denotes the operator norm and for a matrix $X$, $P_X = X(X^TX)^{-1}X^T$ is the matrix projecting into the column space of $X$. Note that $$AA^T = A(A^TA)^{-1}A^T = P_A$$ and similarly for $B$, so since we have orthonormal columns we get $$d(A,B) = \|AA^T-BB^T\|.$$

So suppose $A$ and $B$ have exactly the same column spaces. For a matrix $M$, $\|M\|^2 = \lambda_{\max}(M^TM)$ so consider $$(P_A - P_B)^T(P_A - P_B) = P_A + P_B - P_BP_A - P_AP_B$$ by the projection matrices being symmetric and idempotent. $P_A$ is completely in the column space of $A$, and similarly for $P_B$, so if $A$ and $B$ have the same column space then $P_BP_A = P_B$ and $P_AP_B = P_A$ in which case $$P_A + P_B - P_BP_A - P_AP_B = 0.$$ To check this, you can compute that $B_1B_1^T -B_2B_2^T = 0$.

Furthermore, let $Q$ be a rotation of the columns and apply it to $B$, say, via $BQ$. Then $$P_A - P_{BQ} = AA^T - BQ(Q^TB^TBQ)^{-1}Q^TB^T = P_A - P_B$$ so this metric is invariant to rotations of the basis vectors. That's a nice property for this sort of thing.

This $d$ appears as a metric on the Grassmanian (see here in particular) although I know very little about the general theory of these things, but that could give you somewhere to look for more.

Another reference that would probably be useful for working with this kind of thing is Stewart and Sun's Matrix Perturbation Theory, since the effect of noise on projections is exactly the kind of thing that matrix perturbation theory studies.

• +1. Interesting answer. If one takes Frobenius norm for $d()$ instead of the operator norm, then this is the same measure as considered here stats.stackexchange.com/questions/141611. Do you have any intuition about the norm choice here? – amoeba says Reinstate Monica Jun 18 '18 at 19:04
• @amoeba thanks for that link, very interesting. I'll have to think about it and will hopefully update this evening. My first thought is that since $$\|P_A - P_B\| = \lambda_{\max}(AA^T-BB^T) = \lambda_{\max}(AA^T+BB^T-AA^TBB^T-BB^TAA^T)$$ I can bound my $d$ as $$\lambda_{\max}(AA^T-BB^T) \leq 2 -2\inf_{x : \|x\|=1} x^TAA^TBB^Tx$$ so that's a bound in terms of just an extreme eigenvalue of $AA^TBB^T$ instead of $\text{tr}(AA^TBB^T)$ which uses all eigenvalues of $AA^TBB^T$. But I bet there's a lot more to be said and I'd really like to see this $d$ better connected to principal angles – jld Jun 18 '18 at 19:37
• just noticed there's a typo, it should be $\|P_A-P_B\| = \lambda_{\max}\left((AA^T-BB^T)^T(AA^T-BB^T)\right)$ – jld Jun 18 '18 at 19:42
• Before I accept the answer (thanks for the extensive reasoning, btw!), I wanted to poke the issue of "less important" PCs (i.e. the eigenvectors corresponding to the smaller eigenvalues). If I get that correctly, the distance you are proposing is more a general distance metric (although the seemingly more correct version). I think I could probably also compute the distance on a subset of the vectors to get what i want, i guess... – dennlinger Jun 19 '18 at 5:46
• @dennlinger I'd have to take some time to investigate the behavior of this function, but maybe it's interesting. Let $\lambda_1\geq\dots\geq\lambda_m$ be the distinct eigenvalues of the covariance matrix in question, and let $A_i$ and $B_i$ be eigenbases of $\lambda_i$ for $i=1,\dots,m$. For a real-valued and symmetric matrix, eigenvectors of different eigenvalues are orthogonal so I think it's sensible enough to consider the projections separately. Then consider $$\sum_{i=1}^m \lambda_i \|P_{A_i} - P_{B_i}\|$$as the metric, where it's the eigenspace differences weighted by the eigenvalues – jld Jun 20 '18 at 3:11