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 |
---|