# Find first Principal Component (and loading) using a fast iterative algorithm without covariance matrix

I have a matrix $$X$$ and I would like to find its first principal component and the corresponding loadings. I would like to do this without computing the covariance matrix of $$X$$. How can I do so?

This is the standard version, which uses the eigendecomposition of the covariance matrix.

using LinearAlgebra: eigen
using Statistics: mean

function find_principal_component(X)
n = size(X, 1)
B = X .- mapslices(mean, X, dims=)     # Center columns of X
evalues, V = eigen(B'B / (n - 1))         # EigenDecomposition of Covariance Matrix
return B * PC, PC
end


Alternatively, one could use the power method, which still uses the covariance matrix

function power_method(X, niter=50)
pc = randn(size(X, 2))
pc /= norm(pc)
M = X'X
for i in 1:niter
pc = M * pc
pc /= norm(pc)
end
return X * pc, pc
end


I would like something like the power method, but without needing to compute the covariance matrix, which can be quite costly.

# Possible solution

I noticed something interesting. Let $$r_t$$ be the principal component vector at time $$t$$. The idea of the power method is to start with a random $$r_t$$ and multiply it by $$X^\top X$$ many times to stretch it towards the principal component. In other words $$r_{t+1} = X^\top X r_t$$ Once we have the principal component $$r_t$$ then the loadings are simply $$\ell_t = X r_t$$ This means we can write $$r_{t+1} = X^\top \ell_t$$ One could therefore start with $$r_t$$ and $$\ell_t$$ initialized randomly and then do \begin{align} r_{t+1} &= \widehat{X^\top \ell_t} \\ \ell_{t+1} &= X r_{t+1} \end{align}

• Why is your proposed solution an answer to your question? Your method involves computing a product with X’X at each iteration, which is a scalar multiple of the covariance. I thought you wanted to avoid forming the covariance matrix.
– Sycorax
Dec 2, 2021 at 17:21
• Just use truncated SVD: it already exists in Julia: github.com/JuliaLinearAlgebra/TSVD.jl Dec 2, 2021 at 17:37
• @Physics_Student but certainly doing one matrix-vector product would be better than doing two, right?
– Sycorax
Dec 2, 2021 at 17:41
• Literature on this should be around what is called "Lanczos Bidiagonalization" methods. You might want to invest time to reimplement it just so you understand what it does. Dec 2, 2021 at 17:48
• Finding structure with randomness: Stochastic algorithms for constructing approximate matrix decompositions Halko, et al., 2009 (arXiv:909) arxiv.org/pdf/0909.4061.pdf
– Sycorax
Dec 2, 2021 at 18:01

Suppose $$X$$ is mean-centered (you have subtracted the mean of each column) with the columns storing the features and the rows storing the observations. PCA is the eigendecomposition of the covariance matrix $$\Sigma = \frac{1}{n-1}X^T X$$.
SVD does not require forming $$\Sigma$$, and you can use power iteration to compute the singular vector to the largest singular value from $$X$$ directly. Note that this will work best if the largest singular value is much larger than all other singular values.
This is simply a naive approach using tools that you've already outlined. I'm sure there are better ones, perhaps exploiting specific knowledge about $$X$$.