对于大小为 n×p 的数据矩阵,PRINCOMP
将返回大小为 p×p 的系数矩阵,其中每一列都是使用原始维度表示的主成分,因此在您的情况下,您将创建一个大小为的输出矩阵:
1036800*1036800*8 bytes ~ 7.8 TB
考虑使用PRINCOMP(X,'econ')
仅返回具有显着差异的 PC
或者,考虑执行SVD 的主成分分析 http://en.wikipedia.org/wiki/Singular_value_decomposition#Relation_to_eigenvalue_decomposition: 在你的情况下n<<p
,并且协方差矩阵无法计算。因此,不要分解 p×p 矩阵XX'
,仅分解较小的 n×n 矩阵就足够了X'X
。参考这张纸 http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/4000/pdf/imm4000.pdf以供参考。
EDIT:
这是我的实现,该函数的输出与PRINCOMP http://www.mathworks.com/help/stats/princomp.html(无论如何前三个):
function [PC,Y,varPC] = pca_by_svd(X)
% PCA_BY_SVD
% X data matrix of size n-by-p where n<<p
% PC columns are first n principal components
% Y data projected on those PCs
% varPC variance along the PCs
%
X0 = bsxfun(@minus, X, mean(X,1)); % shift data to zero-mean
[U,S,PC] = svd(X0,'econ'); % SVD decomposition
Y = X0*PC; % project X on PC
varPC = diag(S'*S)' / (size(X,1)-1); % variance explained
end
我刚刚在我的 4GB 机器上尝试过,运行得很好:
» x = rand(16,1036800);
» [PC, Y, varPC] = pca_by_svd(x);
» whos
Name Size Bytes Class Attributes
PC 1036800x16 132710400 double
Y 16x16 2048 double
varPC 1x16 128 double
x 16x1036800 132710400 double
Update:
The princomp
函数已被弃用,取而代之的是pca http://www.mathworks.com/help/stats/pca.htmlR2012b 中引入,其中包含更多选项。