我之前在 StackOverflow 上问过这个问题，但似乎在这里可能更合适，因为它没有得到任何关于 SO 的答案。它有点像统计和编程之间的交集。
我需要编写一些代码来进行 PCA（主成分分析）。我浏览了著名的算法并实现了这个，据我所知，它相当于 NIPALS 算法。它可以很好地找到前 2-3 个主成分，但随后收敛似乎变得非常缓慢（大约数百到数千次迭代）。以下是我需要的详细信息：
在处理大量特征（大约 10,000 到 20,000 个）和数百个样本大小时，该算法必须高效。
它必须在没有像样的线性代数/矩阵库的情况下可以合理实现，因为目标语言是 D，它还没有，即使有，我也不希望将它作为依赖项添加到相关项目中.
附带说明一下，在同一个数据集上，R 似乎可以非常快地找到所有主成分，但它使用奇异值分解，这不是我想自己编写的代码。
我已经实现了“Halko, N., Martinsson, PG, Shkolnisky, Y., & Tygert, M. (2010) 中给出的随机 SVD。一种用于大型数据集主成分分析的算法。Arxiv 预印本 arXiv： 1007.5510, 0526。检索于 2011 年 4 月 1 日，来自http://arxiv.org/abs/1007.5510。”。如果你想获得截断的 SVD，它确实比 MATLAB 中的 svd 变体快得多。你可以在这里得到它：
function [U,S,V] = fsvd(A, k, i, usePowerMethod)
% FSVD Fast Singular Value Decomposition
%
% [U,S,V] = FSVD(A,k,i,usePowerMethod) computes the truncated singular
% value decomposition of the input matrix A upto rank k using i levels of
% Krylov method as given in [1], p. 3.
%
% If usePowerMethod is given as true, then only exponent i is used (i.e.
% as power method). See [2] p.9, Randomized PCA algorithm for details.
%
% [1] Halko, N., Martinsson, P. G., Shkolnisky, Y., & Tygert, M. (2010).
% An algorithm for the principal component analysis of large data sets.
% Arxiv preprint arXiv:1007.5510, 0526. Retrieved April 1, 2011, from
% http://arxiv.org/abs/1007.5510.
%
% [2] Halko, N., Martinsson, P. G., & Tropp, J. A. (2009). Finding
% structure with randomness: Probabilistic algorithms for constructing
% approximate matrix decompositions. Arxiv preprint arXiv:0909.4061.
% Retrieved April 1, 2011, from http://arxiv.org/abs/0909.4061.
%
% See also SVD.
%
% Copyright 2011 Ismail Ari, http://ismailari.com.
if nargin < 3
i = 1;
end
% Take (conjugate) transpose if necessary. It makes H smaller thus
% leading the computations to be faster
if size(A,1) < size(A,2)
A = A';
isTransposed = true;
else
isTransposed = false;
end
n = size(A,2);
l = k + 2;
% Form a real n×l matrix G whose entries are iid Gaussian r.v.s of zero
% mean and unit variance
G = randn(n,l);
if nargin >= 4 && usePowerMethod
% Use only the given exponent
H = A*G;
for j = 2:i+1
H = A * (A'*H);
end
else
% Compute the m×l matrices H^{(0)}, ..., H^{(i)}
% Note that this is done implicitly in each iteration below.
H = cell(1,i+1);
H{1} = A*G;
for j = 2:i+1
H{j} = A * (A'*H{j-1});
end
% Form the m×((i+1)l) matrix H
H = cell2mat(H);
end
% Using the pivoted QR-decomposiion, form a real m×((i+1)l) matrix Q
% whose columns are orthonormal, s.t. there exists a real
% ((i+1)l)×((i+1)l) matrix R for which H = QR.
% XXX: Buradaki column pivoting ile yapılmayan hali.
[Q,~] = qr(H,0);
% Compute the n×((i+1)l) product matrix T = A^T Q
T = A'*Q;
% Form an SVD of T
[Vt, St, W] = svd(T,'econ');
% Compute the m×((i+1)l) product matrix
Ut = Q*W;
% Retrieve the leftmost m×k block U of Ut, the leftmost n×k block V of
% Vt, and the leftmost uppermost k×k block S of St. The product U S V^T
% then approxiamtes A.
if isTransposed
V = Ut(:,1:k);
U = Vt(:,1:k);
else
U = Ut(:,1:k);
V = Vt(:,1:k);
end
S = St(1:k,1:k);
end
要测试它，只需在同一个文件夹中创建一个图像（就像一个大矩阵，你可以自己创建矩阵）
% Example code for fast SVD.
clc, clear
%% TRY ME
k = 10; % # dims
i = 2; % # power
COMPUTE_SVD0 = true; % Comment out if you do not want to spend time with builtin SVD.
% A is the m×n matrix we want to decompose
A = im2double(rgb2gray(imread('test_image.jpg')))';
%% DO NOT MODIFY
if COMPUTE_SVD0
tic
% Compute SVD of A directly
[U0, S0, V0] = svd(A,'econ');
A0 = U0(:,1:k) * S0(1:k,1:k) * V0(:,1:k)';
toc
display(['SVD Error: ' num2str(compute_error(A,A0))])
clear U0 S0 V0
end
% FSVD without power method
tic
[U1, S1, V1] = fsvd(A, k, i);
toc
A1 = U1 * S1 * V1';
display(['FSVD HYBRID Error: ' num2str(compute_error(A,A1))])
clear U1 S1 V1
% FSVD with power method
tic
[U2, S2, V2] = fsvd(A, k, i, true);
toc
A2 = U2 * S2 * V2';
display(['FSVD POWER Error: ' num2str(compute_error(A,A2))])
clear U2 S2 V2
subplot(2,2,1), imshow(A'), title('A (orig)')
if COMPUTE_SVD0, subplot(2,2,2), imshow(A0'), title('A0 (svd)'), end
subplot(2,2,3), imshow(A1'), title('A1 (fsvd hybrid)')
subplot(2,2,4), imshow(A2'), title('A2 (fsvd power)')
当我在桌面上运行它以获得大小为 635*483 的图像时，我得到
Elapsed time is 0.110510 seconds.
SVD Error: 0.19132
Elapsed time is 0.017286 seconds.
FSVD HYBRID Error: 0.19142
Elapsed time is 0.006496 seconds.
FSVD POWER Error: 0.19206
如您所见，对于 的低值k，它比使用 Matlab SVD 快 10 倍以上。顺便说一句，您可能需要以下简单的功能来测试功能：
function e = compute_error(A, B)
% COMPUTE_ERROR Compute relative error between two arrays
e = norm(A(:)-B(:)) / norm(A(:));
end
我没有添加 PCA 方法，因为使用 SVD 实现起来很简单。您可以查看此链接以查看他们的关系。
你可以尝试使用几个选项。
1-惩罚矩阵分解。您对 u 和 v 应用一些惩罚约束以获得一些稀疏性。已用于基因组数据的快速算法
见惠顿提布希拉尼。他们也有一个 R-pkg。“一种惩罚矩阵分解，适用于稀疏主成分和典型相关分析。”
2-随机 SVD。由于 SVD 是一种主要算法，因此可能需要找到一个非常快速的近似值，尤其是对于探索性分析。使用随机 SVD，您可以对庞大的数据集进行 PCA。
请参阅 Martinsson、Rokhlin 和 Tygert “矩阵分解的随机算法”。Tygert 有用于非常快速地实现 PCA 的代码。
下面是 R 中随机 SVD 的简单实现。
ransvd = function(A, k=10, p=5) {
n = nrow(A)
y = A %*% matrix(rnorm(n * (k+p)), nrow=n)
q = qr.Q(qr(y))
b = t(q) %*% A
svd = svd(b)
list(u=q %*% svd$u, d=svd$d, v=svd$v)
}
听起来您可能想使用Lanczos 算法。如果做不到这一点，您可能需要咨询Golub & Van Loan。我曾经从他们的文本中编写了一个 SVD 算法（在所有语言的 SML 中），它运行得相当好。
我建议尝试内核 PCA，它的时间/空间复杂度取决于示例数 (N) 而不是特征数 (P)，我认为这更适合您的设置 (P>>N)。内核 PCA 基本上与 NxN 内核矩阵（数据点之间的相似性矩阵）一起使用，而不是 PxP 协方差矩阵，后者对于大 P 来说很难处理。内核 PCA 的另一个好处是它可以学习非线性投影以及如果您将它与合适的内核一起使用。请参阅有关内核 PCA 的这篇论文。