1 | function [U,S,V] = svdj(A, tol) |
---|
2 | % SVDJ Singular value decomposition using Jacobi algorithm. |
---|
3 | |
---|
4 | % Copyright © 2005, 2006 Nicolas Le Bihan and Stephen J. Sangwine. |
---|
5 | |
---|
6 | % This function works for real, complex or quaternion matrices. |
---|
7 | % |
---|
8 | % Arguments: A - a real, complex or quaternion matrix. |
---|
9 | % tol - a tolerance, defaults to eps if omitted. |
---|
10 | % Returns: U, V - singular vectors (unitary matrices). |
---|
11 | % S - singular values. |
---|
12 | % |
---|
13 | % The singular value decomposition computed by this function is the same as |
---|
14 | % that computed by the function svd.m, but this code, being based on a |
---|
15 | % cyclic Jacobi algorithm, is more accurate. However, it is also slower. |
---|
16 | % It is intended as a reference implementation since the Jacobi algorithm |
---|
17 | % is known to be the most accurate SVD algorithm. |
---|
18 | % |
---|
19 | % This function will work on real, complex and quaternion matrices. |
---|
20 | % |
---|
21 | % References: |
---|
22 | % |
---|
23 | % N. Le Bihan and S. J. Sangwine, "Jacobi Method for Quaternion Matrix |
---|
24 | % Singular Value Decomposition", Applied Mathematics and Computation, 2006. |
---|
25 | % DOI:10.1016/j.amc.2006.09.055. |
---|
26 | % |
---|
27 | % S. J. Sangwine and N. Le Bihan, "Computing the SVD of a quaternion matrix", |
---|
28 | % in Seventh International Conference on Mathematics in Signal Processing, |
---|
29 | % 17-20 December 2006, The Royal Agricultural College, Cirencester, UK. |
---|
30 | % Institute of Mathematics and its Applications, 2006. |
---|
31 | |
---|
32 | error(nargchk(1, 2, nargin)), error(nargoutchk(0, 3, nargout)) |
---|
33 | |
---|
34 | if nargin == 1 |
---|
35 | tol = eps; % Default value for the tolerance. |
---|
36 | end |
---|
37 | |
---|
38 | if nargout == 2 |
---|
39 | error('The number of output parameters must be 0, 1 or 3'); |
---|
40 | end |
---|
41 | |
---|
42 | % We should verify here that A is a real/complex (i.e. double) or |
---|
43 | % quaternion matrix. The code cannot work for arbitrary datatypes. |
---|
44 | |
---|
45 | if ~isreal(A) && isa(A,'quaternion') |
---|
46 | error('svd_jacobi does not work with complex quaternion matrices'); |
---|
47 | end |
---|
48 | |
---|
49 | [M,N] = size(A); K = min(M,N); % K is the number of singular values. |
---|
50 | |
---|
51 | % In what follows we need to be able to construct a quaternion or real or |
---|
52 | % complex matrix according to the type of the actual supplied for the |
---|
53 | % parameter A. This is a tricky bit of programming, but the key to it is |
---|
54 | % the concept of function handles (see Matlab Help). This permits us to |
---|
55 | % call the constructor function for the type of A, which is double in the |
---|
56 | % case of real or complex A, quaternion if A is quaternion. |
---|
57 | |
---|
58 | F = str2func(class(A)); % F is a function handle. |
---|
59 | |
---|
60 | V = F(eye(N)); |
---|
61 | |
---|
62 | % Calculate the sum of the moduli of the diagonal elements of the |
---|
63 | % implicit matrix B = A' * A. This sum is invariant and we do not |
---|
64 | % need to calculate it again. We normalise it by the matrix size. |
---|
65 | |
---|
66 | On = 0; for c = A, On = On + sum(abs(c).^2); end; On = On ./ N; |
---|
67 | |
---|
68 | Previous_Off = Inf; % We test on each sweep to make sure the Off diagonal |
---|
69 | % sum is reducing. If it does not reduce we stop. We |
---|
70 | % use infinity as the initial value. |
---|
71 | |
---|
72 | while true % This is really a repeat .. until, but since Matlab does not |
---|
73 | % provide this statement, we use an if .. break at the end of |
---|
74 | % the loop. |
---|
75 | |
---|
76 | % Sweep through the upper triangular part of the implicit matrix B. |
---|
77 | |
---|
78 | R = 0; % We count the rotations so that we know if we have not done any |
---|
79 | % during a whole sweep. |
---|
80 | |
---|
81 | for r = 1 : N - 1 |
---|
82 | for c = r + 1 : N |
---|
83 | |
---|
84 | % Calculate the three elements of the implicit matrix B that are |
---|
85 | % needed to calculate a Jacobi rotation. Since B is Hermitian, the |
---|
86 | % fourth element (b_cr) is not needed. |
---|
87 | |
---|
88 | b_rr = sum(abs(A(:,r)).^2); % Real value. |
---|
89 | b_cc = sum(abs(A(:,c)).^2); % Real value. |
---|
90 | b_rc = A(:,r)' * A(:,c); % Same type as A. |
---|
91 | |
---|
92 | % Calculate a Jacobi rotation (four elements of G). The two values |
---|
93 | % that we calculate are a real value, C = cos(theta) and S, a value |
---|
94 | % of the same type as A, such that |S| = sin(theta). |
---|
95 | |
---|
96 | m = abs(b_rc); |
---|
97 | |
---|
98 | if m ~= 0 % If the off-diagonal element is zero, we don't rotate. |
---|
99 | |
---|
100 | tau = (b_cc - b_rr) / (2 * m); % tau is real and will be zero if |
---|
101 | % the two on-diagonal elements are |
---|
102 | % equal. In this case G will be an |
---|
103 | % identity matrix, and there is no |
---|
104 | % point in further calculating it. |
---|
105 | if tau ~= 0 |
---|
106 | |
---|
107 | R = R + 1; % Count the rotation we are about to perform. |
---|
108 | |
---|
109 | t = sign(tau) ./ (abs(tau) + sqrt(1 + tau .^ 2)); |
---|
110 | C = 1 ./ sqrt(1 + t .^ 2); |
---|
111 | S = (b_rc .* t .* C) ./ m; |
---|
112 | |
---|
113 | % Initialize the rotation matrix, which is the same size as the |
---|
114 | % implicit matrix B. |
---|
115 | |
---|
116 | % We have to create an identity matrix here of the same type as A, |
---|
117 | % that is, quaternion if A is a quaternion, double if A is double. |
---|
118 | % To do this we use a function handle (q.v.) constructed from the |
---|
119 | % class type of A. This was done before the loop, since the type |
---|
120 | % of A is invariant. |
---|
121 | |
---|
122 | G = F(eye(N)); |
---|
123 | |
---|
124 | G(r,r) = F(C); |
---|
125 | G(c,c) = F(C); |
---|
126 | G(r,c) = S; |
---|
127 | G(c,r) =-conj(S); |
---|
128 | |
---|
129 | % Update of A and V. This is apparently terrible as it performs a |
---|
130 | % full matrix multiplication and most of G is zero. An alternative |
---|
131 | % is to multiply only the row and column of A that are affected, |
---|
132 | % but this turns out to work slower by about a factor of two. See |
---|
133 | % the (commented) code below. Clearly Matlab exploits the near |
---|
134 | % identity structure of G in some way, perhaps because it uses a |
---|
135 | % special (sparse?) structure to store a near identity matrix. |
---|
136 | % |
---|
137 | % Temp = A(:,r) .* G(r,r) + A(:,c) .* G(c,r); |
---|
138 | % A(:,c) = A(:,r) .* G(r,c) + A(:,c) .* G(c,c); |
---|
139 | % A(:,r) = Temp; |
---|
140 | |
---|
141 | A = A * G; |
---|
142 | V = V * G; |
---|
143 | end |
---|
144 | end |
---|
145 | end |
---|
146 | end |
---|
147 | |
---|
148 | if R == 0 % then there were no rotations on the last sweep... |
---|
149 | |
---|
150 | % This condition can occur with pathological matrices and there must |
---|
151 | % be some fix to enable the SVD to be calculated. For the moment at |
---|
152 | % least the problem is detected - prior code would enter an infinite |
---|
153 | % loop under this condition. |
---|
154 | |
---|
155 | error('No rotations performed during sweep.') |
---|
156 | end |
---|
157 | |
---|
158 | % Calculate the sum of the off-diagonal elements of the matrix B. |
---|
159 | |
---|
160 | B = A' * A; |
---|
161 | |
---|
162 | Off = sum(sum(abs(triu(B, 1))))/(N.^2); % Normalise by the matrix size! |
---|
163 | |
---|
164 | if (Off/On) < tol |
---|
165 | break; % The off-diagonal sum is small enough for us to stop. |
---|
166 | end |
---|
167 | |
---|
168 | if Previous_Off < Off |
---|
169 | warning('Terminating sweeps: off diagonal sum increased on last sweep.') |
---|
170 | break; |
---|
171 | end; |
---|
172 | |
---|
173 | Previous_Off = Off; |
---|
174 | |
---|
175 | end |
---|
176 | |
---|
177 | % Extract and sort the singular values. The vector T may be longer than the |
---|
178 | % number of singular values (K) in cases where A is not square. |
---|
179 | |
---|
180 | [T,IX] = sort(sqrt(abs(diag(B))),'descend'); |
---|
181 | |
---|
182 | if nargout == 0 || nargout == 1 % .. only the singular values are needed. |
---|
183 | U = T(1:K); |
---|
184 | end |
---|
185 | |
---|
186 | if nargout == 3 % .. the singular vectors and singular values are needed. |
---|
187 | |
---|
188 | A = A(:, IX); % Map the columns of A and V into the same order as the |
---|
189 | V = V(:, IX); % singular values, using the sort indices in IX. |
---|
190 | |
---|
191 | % Construct the left singular vectors. These are in A but we need |
---|
192 | % to divide each column by the corresponding singular value. This |
---|
193 | % calculation is done by replicating T to make a matrix which can |
---|
194 | % then be divided into A element-by-element (vectorized division). |
---|
195 | |
---|
196 | U = A ./ repmat(T',M,1); |
---|
197 | |
---|
198 | S = diag(T); % Construct a diagonal matrix of singular values from |
---|
199 | % the vector T, because when there are three output |
---|
200 | % parameters, S is required to be a matrix. |
---|
201 | end |
---|