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


*

*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

*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.
 A: 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.
