julia - 如何在不使用协方差矩阵的情况下快速找到第一个主成分(和载荷)？

Question

我有一个矩阵 $X$，我想找到它的第一个主成分和相应的载荷。我想在不计算 $X$ 的协方差矩阵的情况下执行此操作。我怎样才能这样做？

这是标准版本，它使用协方差矩阵的特征分解。

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

或者，可以使用幂方法，它仍然使用协方差矩阵

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     

我想要类似幂方法的方法，但不需要计算协方差矩阵，这可能会非常昂贵。

可能的解决方案

我注意到一些有趣的事情。让r_t为 t 时刻的主分量向量。幂法的思想是从随机r_t开始并将其乘以 X' X多次将其拉伸(stretch)向主成分。换句话说r_{t+1} = X' X r_t 一旦我们有了主成分r_t那么负载就是\ell_t = X r_t 。这意味着我们可以写 r_{t+1} = X^\top \ell_t 因此，可以从 r_t 开始和\ell_t随机初始化然后执行

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

Answer 1

一般来说，您可能会发现奇异值分解对此更有用。

奇异值分解的定义为

B = U Σ V'

这意味着

B'B = V Σ² V'

因此，您的代码可以避免 B'B 的计算。更重要的是，奇异值始终是实数，因此您不必担心 B'B 是否完全对称。

更好的是，Arpack.svds 允许您计算最大的几个奇异值。

这是使用 SVD 而不是特征分解的代码版本:

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

在一个大型稀疏矩阵(100k x 1k，1M 非零)上运行它可以得到这样的速度:

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, ...])

在一个大型非稀疏示例(1M x 100 个条目)上:

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, ...])