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=[1])     # Center columns of X
    evalues, V = eigen(B'B / (n - 1))         # EigenDecomposition of Covariance Matrix     
    PC = V[:, argmax(evalues)]                # Grab principal component and compute loading
    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}
$$

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$.
There is a deep relationship between PCA and SVD. In fact, you can use SVD to compute PCA. See: Relationship between SVD and PCA. How to use SVD to perform PCA?
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$.
