-
Notifications
You must be signed in to change notification settings - Fork 0
/
symmetric_orthogonalise_custom.m
108 lines (91 loc) · 3.68 KB
/
symmetric_orthogonalise_custom.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
96
97
98
99
100
101
102
103
104
105
106
107
108
function [L, ignore, ignore2, W] = symmetric_orthogonalise_custom(A, maintainMagnitudes)
%SYMMETRIC_ORTHOGONALISE closest orthogonal matrix
%
% L = SYMMETRIC_ORTHOGONALISE(A) returns orthonormal matrix L which
% is closest to A, as measured by the Frobenius norm of (L-A).
%
% The orthogonal matrix is constructed from a singular value decomposition
% of A.
%
% L = SYMMETRIC_ORTHOGONALISE(A, KEEP_MAGNITUDES) returns the orthogonal
% matrix L, whose columns have the same magnitude as the respective
% columns of A, and which is closest to A, as measured by the Frobenius
% norm of (L-A), if KEEP_MAGNITUDES is TRUE.
%
% The orthogonal matrix is constructed from a singular value decomposition
% of A.
%
% [L, ~, ~, W] = SYMMETRIC_ORTHOGONALISE(...) also returns a weighting matrix
% such that L = A * W;
%
% See also: ROINETS.HOUSEHOLDER_ORTHOGONALISE, ORTH, SVD.
% References: Naidu, A. R. "Centrality of Lowdin Orthogonalizations",
% arXiv 1105.3571v1, May 2011.
% Available at: http://arxiv-web3.library.cornell.edu/pdf/1105.3571v1.pdf
% Copyright 2013 OHBA
% This program is free software: you can redirstribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% $LastChangedBy: [email protected] $
% $Revision: 239 $
% $LastChangedDate: 2014-08-15 14:58:49 +0100 (Fri, 15 Aug 2014) $
% Contact: giles.colclough 'at' magd.ox.ac.uk
% Originally written on: GLNXA64 by Giles Colclough, 31-Oct-2013 13:30:05
if nargin < 2 || ~exist('maintainMagnitudes', 'var'),
maintainMagnitudes = false;
end
[ignore, ignore2] = deal([]); % to match up with other functions
if maintainMagnitudes,
D = diag(sqrt(diag(A' * A)));
if nargout > 1,
% call function again
[Lnorm, ~, ~, W] = symmetric_orthogonalise_custom(A * D, false);
% scale result
L = Lnorm * D;
W = D * W * D;
else
% call function again
Lnorm = symmetric_orthogonalise_custom(A * D, false);
% scale result
L = Lnorm * D;
end%if
else
[U, S, V] = svd(A, 'econ');
if ~isempty(S),
% we need to check that we have sufficient rank
S = diag(S);
tol = max(size(A)) * S(1) * eps(class(A));
r = sum(S > tol);
isFullRank = (r >= size(A,1));
if isFullRank,
% polar factors of A
L = U * conj(transpose(V));
if nargout > 1,
W = V * diag(1.0 ./ S) * conj(transpose(V));
end%if
else % not enough degrees of freedom
error_message = ['The ROI time-course matrix is not full rank. \n', ...
' This prevents you from using an all-to-all ', ...
'orthogonalisation method. \n', ...
' Your data have rank %d, and you are looking ', ...
'at %d ROIs. \n', ...
' You could try reducing the number of ROIs, or ', ...
'using an alternative orthogonalisation method. \n'];
error('ROInets:RankError',error_message,r,ROInets.cols(A));
end%if
else
L = [];
W = [];
end%if
end%if
end%symmetric_orthogonalise
% [EOF]