1 | function [U, B, V] = bidiagonalize(A) |
---|
2 | % Bidiagonalize A, such that U * A * V = B and U' * B * V' = A. B is the |
---|
3 | % same size as A, has no vector part, and is upper or lower bidiagonal |
---|
4 | % depending on its shape. U and V are unitary quaternion matrices. |
---|
5 | |
---|
6 | % Copyright © 2005 Stephen J. Sangwine and Nicolas Le Bihan. |
---|
7 | % See the file : Copyright.m for further details. |
---|
8 | |
---|
9 | error(nargchk(1, 1, nargin)), error(nargoutchk(3, 3, nargout)) |
---|
10 | |
---|
11 | % References: |
---|
12 | % |
---|
13 | % Sangwine, S. J. and Le Bihan, N., |
---|
14 | % Quaternion singular value decomposition based on bidiagonalization |
---|
15 | % to a real or complex matrix using quaternion Householder transformations, |
---|
16 | % Applied Mathematics and Computation, 182(1), 1 November 2006, 727-738, |
---|
17 | % DOI:10.1016/j.amc.2006.04.032. |
---|
18 | % |
---|
19 | % Sangwine, S. J. and Le Bihan, N., |
---|
20 | % Quaternion Singular Value Decomposition based on Bidiagonalization |
---|
21 | % to a Real Matrix using Quaternion Householder Transformations, |
---|
22 | % arXiv:math.NA/0603251, 10 March 2006, available at http://www.arxiv.org/ |
---|
23 | |
---|
24 | % NB: This is a reference implementation. It uses explicit Householder |
---|
25 | % matrices. Golub and van Loan, 'Matrix Computations', 2e, 1989, section |
---|
26 | % 5.1.4 discusses efficient calculation of Householder transformations. |
---|
27 | % This has been tried (see bidiagonalize2.m), but found to be slower than |
---|
28 | % the use of explicit matrices as used here. This requires further study. |
---|
29 | |
---|
30 | [r, c] = size(A); |
---|
31 | |
---|
32 | if prod([r, c]) == 1 |
---|
33 | error('Cannot bidiagonalize a matrix of one element.'); |
---|
34 | end |
---|
35 | |
---|
36 | if c <= r |
---|
37 | [U, B, V] = internal_bidiagonalizer(A); % Gives an upper bidiagonal result. |
---|
38 | else |
---|
39 | % This requires a lower bidiagonal result. We handle this by a recursive |
---|
40 | % call on the Hermitian transpose of A. The results for U and V must be |
---|
41 | % interchanged and B must be transposed to get the correct result for A. |
---|
42 | |
---|
43 | [V, B, U] = internal_bidiagonalizer(A'); B = B.'; |
---|
44 | end |
---|
45 | |
---|
46 | V = V'; % Transpose and conjugate V for compatibility with earlier code. |
---|
47 | B = check(B); % Verify the result and convert to exactly bidiagonal form. |
---|
48 | |
---|
49 | % ---------------------------------------------------------------------------- |
---|
50 | |
---|
51 | function [U, B, V] = internal_bidiagonalizer(A) |
---|
52 | |
---|
53 | [r, c] = size(A); |
---|
54 | |
---|
55 | % Compute and apply a Householder transformation to the first column of A. |
---|
56 | |
---|
57 | U = householder_matrix(A(:, 1), eye(r, 1)); B = U * A; |
---|
58 | V = quaternion(eye(c)); |
---|
59 | |
---|
60 | % If there is more than one column, we now need to transform the first row (excluding the |
---|
61 | % first element). A recursive call on the transposed conjugate matrix does this. The left |
---|
62 | % and right unitary results are interchanged. |
---|
63 | |
---|
64 | if c > 1 |
---|
65 | [V(2 : end, 2 : end), T, W] = internal_bidiagonalizer(B(:, 2 : end)'); |
---|
66 | B(:, 2 : end) = T.'; |
---|
67 | U = W * U; |
---|
68 | end |
---|
69 | |
---|
70 | % --------------------------------------------------------------------------------------------- |
---|
71 | |
---|
72 | function R = check(B) |
---|
73 | |
---|
74 | % Verify results, and convert the result to exactly bidiagonal form with no vector part. |
---|
75 | |
---|
76 | [r, c] = size(B); |
---|
77 | |
---|
78 | if r == 1 || c == 1 |
---|
79 | % The matrix is degenerate (a row or column vector) and we have to deal |
---|
80 | % with it differently because the Matlab diag function in this case |
---|
81 | % constructs a matrix instead of extracting the diagonal (how clever to |
---|
82 | % use the same name for both ideas!). |
---|
83 | D = B(1); % The first element is the diagonal. There is no super-diagonal. |
---|
84 | O = B(2 : end); % The rest is the off-diagonal. |
---|
85 | elseif c <= r |
---|
86 | D = [diag(B); diag(B, +1)]; % Extract the diagonal and super-diagonal. |
---|
87 | O = tril(B, -1) + triu(B, +2); % Extract the off-diagonal part. |
---|
88 | else |
---|
89 | D = [diag(B); diag(B, -1)]; % Extract the diagonal and sub-diagonal. |
---|
90 | O = tril(B, -2) + triu(B, +1); % Extract the off-diagonal part. |
---|
91 | end |
---|
92 | |
---|
93 | % Find the largest on and off bidiagonal elements. We use abs twice to |
---|
94 | % allow for the case where B is complex (this occurs when A is a complexified |
---|
95 | % quaternion matrix). |
---|
96 | |
---|
97 | T1 = max(max(abs(abs(O)))); % Find the largest off bidiagonal element. |
---|
98 | T2 = max(abs(abs(D))); % Find the largest bidiagonal element. |
---|
99 | |
---|
100 | % NB T2 and/or T1 could be exactly zero (example, if A was zero, or an identity matrix). |
---|
101 | % Therefore we do not divide one by the other, but instead multiply by a tolerance. |
---|
102 | |
---|
103 | tolerance = eps .* 1.0e4; % This is empirically determined. |
---|
104 | |
---|
105 | if T1 > T2 * tolerance |
---|
106 | warning('Result of bidiagonalization was not accurately diagonal.'); |
---|
107 | disp('Information: largest on- and off-bidiagonal moduli were:'); |
---|
108 | disp(T2); |
---|
109 | disp(T1); |
---|
110 | end |
---|
111 | |
---|
112 | % Verify that the diagonal elements have neglible vector parts. |
---|
113 | |
---|
114 | T2 = max(abs(abs(s(D)))); % The largest scalar modulus of the bidiagonal result |
---|
115 | T1 = max(abs(abs(v(D)))); % The largest vector modulus of the bidiagonal result. |
---|
116 | |
---|
117 | if T1 > T2 * tolerance |
---|
118 | warning('Result of bidiagonalization was not accurately scalar.') |
---|
119 | disp('Information: largest on-diagonal scalar and vector moduli were:'); |
---|
120 | disp(T2); |
---|
121 | disp(T1); |
---|
122 | end |
---|
123 | |
---|
124 | if r == 1 || c == 1 |
---|
125 | R = s(B); % The diagonal has only one element, so we can just forget the off-diagonal. |
---|
126 | else |
---|
127 | R = s(B - O); % Subtract the off-diagonal part and take the scalar part of the result. |
---|
128 | end |
---|
129 | |
---|