How can I find first principal component (and loadings) fast without using covariance matrix?

user_6068731 0
Viewed 112

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' X many times to stretch it towards the principal component. In other words r_{t+1} = X' 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

r_{t+1} = normalize(X^\top \ell_t)\\
\ell_{t+1} = X r_{t+1}
julia pca
2 Answers

In general, you may find singular value decompositions more useful for this.

The definition of the singular value decomposition is

B = U Σ V'

This means that

B'B = V Σ² V'

As a result, you code can avoid the computation of B'B. More importantly, the singular values are always real and thus you don't have to worry about whether B'B will be exactly symmetric.

Even better, Arpack.svds allows you to compute just the largest few singular values.

Here is a version of your code that uses SVD instead of eigen decomposition:

using LinearAlgebra: eigen
using Statistics: mean
using Arpack: svds

function find_principal_component(X)
    n = size(X, 1)
    # Center columns of X
    B = X .- mapslices(mean, X, dims=[1]) 
    # Decomposition of Covariance Matrix
    svd,_ = svds(B / (n - 1), nsv=1)        
    # Grab principal component and compute loading
    PC = svd.V[:, 1]
    return B * PC, PC
end

Running this on a large sparse matrix (100k x 1k, 1M non-zeros) gives this speed:

julia> @time find_principal_component(sprandn(100_000, 1_000, 0.01))
 25.529426 seconds (18.45 k allocations: 3.015 GiB, 0.02% gc time)
([0.014242904195824286, 0.10635817357717596, -0.010142643763158442, ...])

and on a large non-sparse example (1M x 100 entries):

julia> @time find_principal_component(randn(1_000_000, 100))
  4.922949 seconds (1.31 k allocations: 2.280 GiB, 0.02% gc time)
([-0.06629858174095249, 0.6996443876327108, -1.1783642870384952, ...])
user_3995510
user_51413280
Related