pr_spm_svd

 computationally efficient SVD (that can handle sparse arguments)
 FORMAT [U,S,V] = spm_svd(X,u);
 X    - {m x n} matrix
 u    - threshold for normalized eigenvalues (default = 1e-6)

 U    - {m x p} singular vectors
 V    - {m x p} singular variates
 S    - {p x p} singular values
___________________________________________________________________________
 @(#)spm_svd.m    2.2 Karl Friston 00/10/10

0001 function [U,S,V] = spm_svd(X,U)
0002 % computationally efficient SVD (that can handle sparse arguments)
0003 % FORMAT [U,S,V] = spm_svd(X,u);
0004 % X    - {m x n} matrix
0005 % u    - threshold for normalized eigenvalues (default = 1e-6)
0006 %
0007 % U    - {m x p} singular vectors
0008 % V    - {m x p} singular variates
0009 % S    - {p x p} singular values
0010 %___________________________________________________________________________
0011 % @(#)spm_svd.m    2.2 Karl Friston 00/10/10
0012 
0013 
0014 % default threshold
0015 %---------------------------------------------------------------------------
0016 if nargin < 2
0017     U = 1e-6;
0018 end
0019 
0020 % deal with sparse matrices
0021 %---------------------------------------------------------------------------
0022 [M N] = size(X);
0023 p     = find(any(X,2));
0024 q     = find(any(X,1));
0025 X     = X(p,q);
0026 
0027 % SVD
0028 %---------------------------------------------------------------------------
0029 [i j s] = find(X);
0030 [m n]   = size(X);
0031 if any(i - j)
0032 
0033     % off-leading diagonal elements - full SVD
0034     %-------------------------------------------------------------------
0035     X     = full(X);
0036     if m > n
0037 
0038         [v S v] = svd(spm_atranspa(X),0);
0039         S       = sparse(S);
0040         s       = diag(S);
0041         s       = s*length(s)/sum(s);
0042         j       = find(s >= U);
0043         v       = v(:,j);
0044         u       = spm_en(X*v);
0045         S       = sqrt(S(j,j));
0046 
0047     elseif m < n
0048 
0049         [u S u] = svd(spm_atranspa(X'),0);
0050         S       = sparse(S);
0051         s       = diag(S);
0052         s       = s*length(s)/sum(s);
0053         j       = find(s >= U);
0054         u       = u(:,j);
0055         v       = spm_en(X'*u);
0056         S       = sqrt(S(j,j));
0057 
0058     else
0059 
0060         [u S v] = svd(X,0);
0061         S       = sparse(S);
0062         s       = diag(S).^2;
0063         s       = s*length(s)/sum(s);
0064         j       = find(s >= U);
0065         v       = v(:,j);
0066         u       = u(:,j);
0067         S       = S(j,j);
0068     end
0069 
0070 else
0071     S     = sparse(1:n,1:n,s,m,n);
0072     u     = speye(m,n);
0073     v     = speye(m,n);
0074     [i j] = sort(-s);
0075     S     = S(j,j);
0076     v     = v(:,j);
0077     u     = u(:,j);
0078 
0079     s     = diag(S).^2;
0080     s     = s*length(s)/sum(s);
0081     j     = find(s >= U);
0082     v     = v(:,j);
0083     u     = u(:,j);
0084     S     = S(j,j);
0085 
0086 end
0087 
0088 % replace in full matrices
0089 %---------------------------------------------------------------------------
0090 j      = length(j);
0091 U      = sparse(M,j);
0092 V      = sparse(N,j);
0093 U(p,:) = u;
0094 V(q,:) = v;