2013年10月15日,星期二

另一种随机技术

在一个 以前的帖子 我讨论了随机特征图,它可以结合内核的功能和原始线性方法的速度。我最近还使用了另一种随机技术来实现正义, 随机SVD 。这是具有许多应用程序的出色原始语言,例如,您可以与随机化的特征图混搭以获得快速的内核PCA,用作双线性潜在因子模型(又称为矩阵分解)的快速初始化器,或利用杠杆作用来计算 巨型CCA .

基本思想是用随机向量探测矩阵,以发现矩阵的低维顶部范围,然后在该空间中执行更便宜的计算。对于平方矩阵,这是直观的:特征向量构成了矩阵的作用仅是按比例缩放的基础,而顶部特征向量具有与之相关的较大比例因子,因此由矩阵缩放的随机向量将成比例地变大。最高的特征方向。这种直觉表明,如果本征谱几乎是平坦的,那么用随机探针捕获顶部本征空间将真的很困难。一般来说,这是对的,如果存在“large spectral gap”,即连续的特征值之间存在较大差异时。但是,即使这有点悲观,因为在机器学习中,有时您并不关心是否获得了子空间“wrong”,例如,如果您要最小化平方重建误差,则几乎相等的特征值会产生较低的后悔。

这是mnist上随机PCA的示例:您需要 以Matlab格式下载mnist 运行这个。
rand('seed',867);
randn('seed',5309);

tic
fprintf('loading mnist');

% get mnist from http://cs.nyu.edu/~roweis/data/mnist_all.mat
load('mnist_all.mat');

trainx=single([train0; train1; train2; train3; train4; train5; train6; train7; train8; train9])/255.0;
testx=single([test0; test1; test2; test3; test4; test5; test6; test7; test8; test9])/255.0;
st=[size(train0,1); size(train1,1); size(train2,1); size(train3,1); size(train4,1); size(train5,1); size(train6,1); size(train7,1); size(train8,1); size(train9,1)];
ss=[size(test0,1); size(test1,1); size(test2,1); size(test3,1); size(test4,1); size(test5,1); size(test6,1); size(test7,1); size(test8,1); size(test9,1)];
paren = @(x, varargin) x(varargin{:});
yt=[]; for i=1:10; yt=[yt; repmat(paren(eye(10),i,:),st(i),1)]; end
ys=[]; for i=1:10; ys=[ys; repmat(paren(eye(10),i,:),ss(i),1)]; end

clear i st ss
clear train0 train1 train2 train3 train4 train5 train6 train7 train8 train9
clear test0 test1 test2 test3 test4 test5 test6 test7 test8 test9

fprintf(' finished: ');
toc

[n,k]=size(yt);
[m,p]=size(trainx);

tic
fprintf('estimating top 50 eigenspace of (1/n) X”X using randomized technique');

d=50;
r=randn(p,d+5);                % NB: we add an extra 5 dimensions here
firstpass=trainx'*(trainx*r);  % this can be done streaming in O(p d) space
q=orth(firstpass);
secondpass=trainx'*(trainx*q); % this can be done streaming in O(p d) space
secondpass=secondpass/n;
z=secondpass'*secondpass;      % note: this is small, i.e., O(d^2) space
[v,s]=eig(z);
pcas=sqrt(s);
pcav=secondpass*v*pinv(pcas);
pcav=pcav(:,end:-1:6);         % NB: and we remove the extra 5 dimensions here
pcas=pcas(end:-1:6,end:-1:6);  % NB: the extra dimensions make the randomized
                               % NB: algorithm more accurate.

fprintf(' finished: ');
toc

tic
fprintf('estimating top 50 eigenspace of (1/n) X”X using  艾格斯 ');

opts.isreal = true; 
[fromeigsv,fromeigss]=eigs(double(trainx'*trainx)/n,50,'LM',opts);

fprintf(' finished: ');
toc


% relative accuracy of eigenvalues
%
% plot((diag(pcas)-diag(fromeigss))./diag(fromeigss))

% largest angle between subspaces spanned by top eigenvectors
% note: can't be larger than pi/2 ~ 1.57
%
% plot(arrayfun(@(x) subspace(pcav(:,1:x),fromeigsv(:,1:x)),linspace(1,50,50))); xlabel('number of eigenvectors'); ylabel('largest principal angle'); set(gca,'YTick',linspace(0,pi/2,5)); 
当我在笔记本电脑上运行时,我得到
>> randsvd
loading mnist finished: Elapsed time is 6.931381 seconds.
estimating top 50 eigenspace of (1/n) X'X using randomized technique finished: Elapsed time is 0.505763 seconds.
estimating top 50 eigenspace of (1/n) X'X using  艾格斯  finished: Elapsed time is 1.051971 seconds.
运行时间的差异不是很大,但是在较大的矩阵上,差异可能是几分钟,而“比您等待的时间更长”。好吧,那有多好?即使假设 艾格斯 是地面真理,有几种方法可以回答这个问题,但是假设我们想从随机技术中获得与 艾格斯 (同样,这在机器学习中通常过于严格)。在这种情况下,我们可以测量最大 主角 在发现的前$ k $个子空间之间 艾格斯 和by the randomized PCA, as a function of $k$. Values near zero indicate that the two subspaces are very nearly identical, whereas values near $\pi / 2$ indicate one subspace contains vectors which are orthogonal to vectors in the other subspace.

In general we see the top 6 or so extracted eigenvectors are spot on, and then it gets worse, better, and worse again. Note it is not monotonic, because if two eigenvectors are reordered, once we have both of them the subspaces will have a small largest 主角 . Roughly speaking anywhere there is a 光谱间隙大 we can expect to get the subspace up to the gap correct, i.e., if there is a flat plateau of eigenvalues followed by a drop than 在 the end of the plateau the largest 主角 should decrease.

Redsvd 提供两遍随机SVD和PCA的开源实现。

2条评论:

  1. 好贴。它'值得注意的是,随机SVD是经典SVD算法(如子空间迭代算法)的早期停止版本。

    回复 删除
  2. 您能否阐明Halko等人的多种算法中的哪一种(可能是4.4的变体)?另外,做'pcas' and 'pcav'这里分别指代奇异值和奇异向量(svd中的S和V),如果是,那么有理由不直接在第二遍而不是第二遍eig()进行svd()'*第二遍,据说SVD在数值上更稳定。谢谢。

    回复 删除