[37] | 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 |
---|