forked from ossadtchi/PSIICOS
-
Notifications
You must be signed in to change notification settings - Fork 0
/
spm_svd.m
95 lines (80 loc) · 2.51 KB
/
spm_svd.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
function [U,S,V] = spm_svd(X,U,T)
% computationally efficient SVD (that can handle sparse arguments)
% FORMAT [U,S,V] = 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) 2008 Wellcome Trust Centre for Neuroimaging
% Karl Friston
% $Id: spm_svd.m 4124 2010-11-18 16:56:53Z karl $
% default thresholds
%---------------------------------------------------------------------------
if nargin < 2, U = 1e-6; end
if nargin < 3, T = 0; end
% deal with sparse matrices
%---------------------------------------------------------------------------
[M N] = size(X);
p = find(any(X,2));
q = find(any(X,1));
X = X(p,q);
% SVD
%---------------------------------------------------------------------------
[i j s] = find(X);
[m n] = size(X);
if any(i - j)
% off-leading diagonal elements - full SVD
%-------------------------------------------------------------------
X = full(X);
if m > n
[v S v] = svd(X'*X,0);
S = sparse(S);
s = diag(S);
j = find(s*length(s)/sum(s) >= U & s >= T);
v = v(:,j);
u = spm_en(X*v);
S = sqrt(S(j,j));
elseif m < n
[u S u] = svd(X*X',0);
S = sparse(S);
s = diag(S);
j = find(s*length(s)/sum(s) >= U & s >= T);
u = u(:,j);
v = spm_en(X'*u);
S = sqrt(S(j,j));
else
[u S v] = svd(X,0);
S = sparse(S);
s = diag(S).^2;
j = find(s*length(s)/sum(s) >= U & s >= T);
v = v(:,j);
u = u(:,j);
S = S(j,j);
end
else
S = sparse(1:n,1:n,s,m,n);
u = speye(m,n);
v = speye(m,n);
[i j] = sort(-s);
S = S(j,j);
v = v(:,j);
u = u(:,j);
s = diag(S).^2;
j = find(s*length(s)/sum(s) >= U & s >= T);
v = v(:,j);
u = u(:,j);
S = S(j,j);
end
% replace in full matrices
%---------------------------------------------------------------------------
j = length(j);
U = sparse(M,j);
V = sparse(N,j);
if j
U(p,:) = u;
V(q,:) = v;
end