我有一个矩阵$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
我想要类似幂法的方法，但不需要计算协方差矩阵，这可能非常昂贵。
可能的解决方案
我注意到一些有趣的事情。设CCD_ 1是在时间CCD_。幂方法的思想是从随机的r_t开始，并将其乘以X' X多次，以将其向主分量拉伸。换句话说，r_{t+1} = X' X r_t一旦我们有了主分量r_t，那么载荷就是简单的ell_t = X r_t。这意味着我们可以编写r_{t+1} = X^top ell_t因此，可以从随机初始化的r_t和r_t0开始，然后进行
r_{t+1} = normalize(X^top ell_t)\
ell_{t+1} = X r_{t+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, ...])
尝试使用KrylovKit.jl。具体而言，eigsolve(X, howmany=1, which=:LM])将为您提供幅度最大的特征值和相关的特征向量。文档位于https://jutho.github.io/KrylovKit.jl/stable/man/eig/