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

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)
    return X * pc, pc

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} $$

  • 1
    $\begingroup$ 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. $\endgroup$
    – Sycorax
    Dec 2, 2021 at 17:21
  • 3
    $\begingroup$ Just use truncated SVD: it already exists in Julia: github.com/JuliaLinearAlgebra/TSVD.jl $\endgroup$
    – usεr11852
    Dec 2, 2021 at 17:37
  • 2
    $\begingroup$ @Physics_Student but certainly doing one matrix-vector product would be better than doing two, right? $\endgroup$
    – Sycorax
    Dec 2, 2021 at 17:41
  • 2
    $\begingroup$ 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. $\endgroup$
    – usεr11852
    Dec 2, 2021 at 17:48
  • 2
    $\begingroup$ Finding structure with randomness: Stochastic algorithms for constructing approximate matrix decompositions Halko, et al., 2009 (arXiv:909) arxiv.org/pdf/0909.4061.pdf $\endgroup$
    – Sycorax
    Dec 2, 2021 at 18:01

1 Answer 1


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$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.