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