[37] | 1 | function [U,S,V] = svd(X, econ) |
---|
| 2 | % SVD Singular value decomposition. |
---|
| 3 | % (Quaternion overloading of standard Matlab function.) |
---|
| 4 | |
---|
| 5 | % Copyright © 2005 Stephen J. Sangwine and Nicolas Le Bihan. |
---|
| 6 | % See the file : Copyright.m for further details. |
---|
| 7 | |
---|
| 8 | % Reference: |
---|
| 9 | % |
---|
| 10 | % S. J. Sangwine and N. Le Bihan, |
---|
| 11 | % Quaternion Singular Value Decomposition based on Bidiagonalization |
---|
| 12 | % to a Real Matrix using Quaternion Householder Transformations, |
---|
| 13 | % arXiv:math.NA/0603251, 10 March 2006. Available at http://www.arxiv.org/. |
---|
| 14 | % |
---|
| 15 | % Note that the title of the article above refers to a real bidiagonal |
---|
| 16 | % matrix, but the algorithm is valid even when X is a complex quaternion |
---|
| 17 | % matrix, in which case the bidiagonal matrix will be complex. The later |
---|
| 18 | % paper below covers both the real and complex cases. |
---|
| 19 | % |
---|
| 20 | % S. J. Sangwine and N. Le Bihan, |
---|
| 21 | % Quaternion Singular Value Decomposition based on Bidiagonalization |
---|
| 22 | % to a Real or Complex Matrix using Quaternion Householder Transformations, |
---|
| 23 | % Applied Mathematics and Computation, in press, 2006. |
---|
| 24 | % DOI:10.1016/j.amc.2006.04.032. Available via http://www.doi.org/. |
---|
| 25 | |
---|
| 26 | error(nargchk(1, 2, nargin)), error(nargoutchk(0, 3, nargout)) |
---|
| 27 | |
---|
| 28 | if nargin == 2 |
---|
| 29 | if econ ~= 0 & econ ~= 'econ' |
---|
| 30 | error('Use svd(X,0) or svd(X,''econ'') for economy size decomposition.'); |
---|
| 31 | end |
---|
| 32 | end |
---|
| 33 | |
---|
| 34 | if numel(X) == 1 |
---|
| 35 | |
---|
| 36 | % The argument is a single quaternion. This case could be handled by |
---|
| 37 | % using the standard Matlab svd() function on the complex adjoint of |
---|
| 38 | % X, but there are problems if X is a complexified quaternion, since |
---|
| 39 | % we cannot make a complex value with complex parts. |
---|
| 40 | % |
---|
| 41 | % For this reason we output an error message and leave it to the user |
---|
| 42 | % to use the appropriate adjoint. |
---|
| 43 | |
---|
| 44 | X = inputname(1); if X == '' X = 'X'; end |
---|
| 45 | |
---|
| 46 | disp('The svd function is not implemented for single quaternions.'); |
---|
| 47 | disp(sprintf('Try using svd(adjoint(%s, ''real'')) or', X)); |
---|
| 48 | disp(sprintf(' svd(adjoint(%s, ''complex'')', X)); |
---|
| 49 | error('Implementation restriction - see advice above'); |
---|
| 50 | end |
---|
| 51 | |
---|
| 52 | % X is not a single quaternion, so proceed ... |
---|
| 53 | |
---|
| 54 | % The method used is to bidiagonalize X using Householder transformations to obtain |
---|
| 55 | % H, B, and G where B is a real or complex bidiagonal matrix, and H and G are the |
---|
| 56 | % products of the Householder matrices used to compute B. Then: H' * B * G' = X. |
---|
| 57 | |
---|
| 58 | [H, B, G] = bidiagonalize(X); |
---|
| 59 | |
---|
| 60 | % The second step is to compute the SVD of B using the standard Matlab routine, |
---|
| 61 | % for a real or complex matrix (which happens to be bidiagonal). |
---|
| 62 | |
---|
| 63 | if nargout == 0 |
---|
| 64 | |
---|
| 65 | % The econ parameter is passed to the Matlab svd routine, even though |
---|
| 66 | % it appears to have no effect if there is no output parameter, since |
---|
| 67 | % only the singular values are output, and they are output as a column |
---|
| 68 | % vector. |
---|
| 69 | |
---|
| 70 | if nargin == 1 |
---|
| 71 | svd(B) |
---|
| 72 | else |
---|
| 73 | svd(B, econ) |
---|
| 74 | end |
---|
| 75 | |
---|
| 76 | elseif nargout == 1 |
---|
| 77 | |
---|
| 78 | % The econ parameter is passed to the Matlab svd routine, even though |
---|
| 79 | % it appears to have no effect if there is only one output parameter, |
---|
| 80 | % because the output is a column vector of the singular values. |
---|
| 81 | |
---|
| 82 | if nargin == 1 U = svd(B); else U = svd(B, econ); end |
---|
| 83 | |
---|
| 84 | else |
---|
| 85 | |
---|
| 86 | % In this case, the econ parameter has been given, but if we were to |
---|
| 87 | % pass it to the Matlab svd routine, the returned U and V matrices |
---|
| 88 | % would not be compatible with H and G. Therefore, we don't pass it, |
---|
| 89 | % but after multiplying the results by H and G, we truncate them to |
---|
| 90 | % give the economy mode result, so that the result is compatible with |
---|
| 91 | % that produced by the standard Matlab svd function. |
---|
| 92 | |
---|
| 93 | [US, S, VS] = svd(B); |
---|
| 94 | |
---|
| 95 | U = H' * US; % Multiply together the intermediate unitary |
---|
| 96 | V = (VS' * G')'; % matrices to construct U and V. |
---|
| 97 | % Notice that US and VS may be complex and that we must |
---|
| 98 | % use ' and not .' here. (H and G are quaternion-valued, |
---|
| 99 | % and therefore H' means a quaternion transpose conjugate.) |
---|
| 100 | |
---|
| 101 | if nargin == 2 |
---|
| 102 | |
---|
| 103 | % The economy size decomposition has been demanded, so we have |
---|
| 104 | % to truncate U or V and S accordingly. The description of how |
---|
| 105 | % the truncation is done can be found in the Matlab help page |
---|
| 106 | % for the standard svd function. |
---|
| 107 | |
---|
| 108 | [m, n] = size(X); |
---|
| 109 | |
---|
| 110 | % In what follows, we need to use subscripted references, but |
---|
| 111 | % these do not work inside a class method (which this is). A |
---|
| 112 | % brief mention of this can be found in the Matlab documentation |
---|
| 113 | % under Programming/Classes and Objects/Designing User Classes |
---|
| 114 | % under the title Object Indexing within Methods. The equally |
---|
| 115 | % obscure solution is to use the substruct function, as below. |
---|
| 116 | % For each statement where it is used, the comment gives the |
---|
| 117 | % normal Matlab notation (which won't work here). Since S is real, |
---|
| 118 | % the normal notation works. |
---|
| 119 | |
---|
| 120 | if econ == 0 |
---|
| 121 | if m > n % If m <= n, there is nothing to do. |
---|
| 122 | U = subsref(U, substruct('()', {':', 1:n})); % U = U(:,1:n) |
---|
| 123 | S = S(1 : n, 1 : n); |
---|
| 124 | end |
---|
| 125 | elseif econ == 'econ' |
---|
| 126 | if m == n |
---|
| 127 | return |
---|
| 128 | elseif m > n % The Matlab documentation says >= here, but if m |
---|
| 129 | % and n are equal all three matrices are square, |
---|
| 130 | % so we don't need to do anything. |
---|
| 131 | U = subsref(U, substruct('()', {':', 1:n})); % U = U(:,1 n) |
---|
| 132 | S = S(1 : n, 1 : n); |
---|
| 133 | else % m < n, since we eliminated m == n above. |
---|
| 134 | S = S(1 : m, 1 : m); |
---|
| 135 | V = subsref(V, substruct('()', {1:n, 1:m})); % V = V(1:n,1:m) |
---|
| 136 | end |
---|
| 137 | else |
---|
| 138 | error('econ parameter has incorrect value'); |
---|
| 139 | end |
---|
| 140 | end |
---|
| 141 | end |
---|