1 | function [U, B, V] = bidiagonalize2(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 | error(nargchk(1, 1, nargin)), error(nargoutchk(3, 3, nargout)) |
---|
7 | |
---|
8 | % References: |
---|
9 | % |
---|
10 | % Sangwine, S. J. and Le Bihan, N., |
---|
11 | % Quaternion singular value decomposition based on bidiagonalization |
---|
12 | % to a real or complex matrix using quaternion Householder transformations, |
---|
13 | % Applied Mathematics and Computation, 182(1), 1 November 2006, 727-738, |
---|
14 | % DOI:10.1016/j.amc.2006.04.032. |
---|
15 | % |
---|
16 | % Sangwine, S. J. and Le Bihan, N., |
---|
17 | % Quaternion Singular Value Decomposition based on Bidiagonalization |
---|
18 | % to a Real Matrix using Quaternion Householder Transformations, |
---|
19 | % arXiv:math.NA/0603251, 10 March 2006, available at http://www.arxiv.org/ |
---|
20 | |
---|
21 | % This is an improved version of the reference implementation. It does not |
---|
22 | % use explicit Householder matrices but instead computes the equivalent |
---|
23 | % matrix update as detailed in section 5.1.4. [p196] of : |
---|
24 | % |
---|
25 | % Gene H. Golub and Charles F. van Loan, 'Matrix Computations', |
---|
26 | % Johns Hopkins University Press, 2nd edition, 1989. |
---|
27 | % |
---|
28 | % This code has not replaced the code in bidiagonalize.m because it is |
---|
29 | % actually slower! It is intended to integrate this code into |
---|
30 | % bidiagonalize.m at a later date, as it will work for matrices where the |
---|
31 | % explicit Householder matrix would be too big. |
---|
32 | |
---|
33 | [r, c] = size(A); |
---|
34 | |
---|
35 | if prod([r, c]) == 1 |
---|
36 | error('Cannot bidiagonalize a matrix of one element.'); |
---|
37 | end |
---|
38 | |
---|
39 | if c <= r |
---|
40 | [U, B, V] = internal_bidiagonalizer(A); % Gives an upper bidiagonal result. |
---|
41 | else |
---|
42 | % This requires a lower bidiagonal result. We handle this by a recursive |
---|
43 | % call on the Hermitian transpose of A. The results for U and V must be |
---|
44 | % interchanged and B must be transposed to get the correct result for A. |
---|
45 | |
---|
46 | [V, B, U] = internal_bidiagonalizer(A'); B = B.'; |
---|
47 | end |
---|
48 | |
---|
49 | V = V'; % Transpose and conjugate V for compatibility with earlier code. |
---|
50 | B = check(B); % Verify the result and convert to exactly bidiagonal real form. |
---|
51 | |
---|
52 | % ---------------------------------------------------------------------------- |
---|
53 | |
---|
54 | function [U, B, V] = internal_bidiagonalizer(A) |
---|
55 | |
---|
56 | [r, c] = size(A); |
---|
57 | |
---|
58 | % Compute and apply a Householder transformation to the first column of A. |
---|
59 | |
---|
60 | % Old code using explicit Householder matrix: |
---|
61 | % U = householder_matrix(A(:, 1), eye(r, 1)); |
---|
62 | |
---|
63 | [u, zeta] = householder_vector(A(:, 1), eye(r, 1)); % New code; |
---|
64 | |
---|
65 | % Old matrix product using explicit Householder matrix U: |
---|
66 | % B = U * A; |
---|
67 | |
---|
68 | B = (1./zeta) .* (A - u * (A' * u)'); % New code. |
---|
69 | |
---|
70 | V = quaternion(eye(c)); |
---|
71 | |
---|
72 | % If there is more than one column, we now need to transform the first row (excluding the |
---|
73 | % first element). A recursive call on the transposed conjugate matrix does this. The left |
---|
74 | % and right unitary results are interchanged. |
---|
75 | |
---|
76 | if c > 1 |
---|
77 | [V(2 : end, 2 : end), T, W] = internal_bidiagonalizer(B(:, 2 : end)'); |
---|
78 | B(:, 2 : end) = T.'; |
---|
79 | |
---|
80 | % Old code, involving explicit Householder matrix: |
---|
81 | % U = W * U; |
---|
82 | |
---|
83 | W = W .* (1./zeta); % New code, step 1: premultiply zeta into W. |
---|
84 | U = W - (W * u) * u'; % New code. |
---|
85 | else |
---|
86 | % c must be 1. We have to form the Householder matrix explicitly in this case. |
---|
87 | % Since c is 1, the matrix has only one element, so this is easy and fast. |
---|
88 | U = (1 ./ zeta) .* (quaternion(1,0,0,0) - u .* conj(u)); |
---|
89 | end |
---|
90 | |
---|
91 | % --------------------------------------------------------------------------------------------- |
---|
92 | |
---|
93 | function R = check(B) |
---|
94 | |
---|
95 | % Verify results, and convert the result to exactly bidiagonal form with no vector part. |
---|
96 | |
---|
97 | [r, c] = size(B); |
---|
98 | |
---|
99 | if r == 1 || c == 1 |
---|
100 | % The matrix is degenerate (a row or column vector) and we have to deal |
---|
101 | % with it differently because the Matlab diag function in this case |
---|
102 | % constructs a matrix instead of extracting the diagonal (how clever to |
---|
103 | % use the same name for both ideas!). |
---|
104 | D = B(1); % The first element is the diagonal. There is no super-diagonal. |
---|
105 | O = B(2 : end); % The rest is the off-diagonal. |
---|
106 | elseif c <= r |
---|
107 | D = [diag(B); diag(B, +1)]; % Extract the diagonal and super-diagonal. |
---|
108 | O = tril(B, -1) + triu(B, +2); % Extract the off-diagonal part. |
---|
109 | else |
---|
110 | D = [diag(B); diag(B, -1)]; % Extract the diagonal and sub-diagonal. |
---|
111 | O = tril(B, -2) + triu(B, +1); % Extract the off-diagonal part. |
---|
112 | end |
---|
113 | |
---|
114 | T1 = max(max(abs(O))); % Find the modulus of the largest off bidiagonal element. |
---|
115 | T2 = max(abs(D)); % Find the modulus of the largest bidiagonal element. |
---|
116 | |
---|
117 | % Note that T1 and T2 may be complex, if A was a complexified quaternion matrix. Therefore we |
---|
118 | % take the modulus of each before comparing them. This has no effect if T1 and 2 are real. NB |
---|
119 | % T2 and/or T1 could be exactly zero (example, if A was zero, or an identity matrix). |
---|
120 | % Therefore we do not divide one by the other, but instead multiply by eps. |
---|
121 | |
---|
122 | if abs(T1) > abs(T2) * eps * 2 |
---|
123 | warning('Result of bidiagonalization was not accurate.'); |
---|
124 | disp('Information: largest on- and off-bidiagonal moduli were:'); |
---|
125 | disp(abs(T2)); |
---|
126 | disp(abs(T1)); |
---|
127 | end |
---|
128 | |
---|
129 | % Verify that the diagonal elements have neglible vector parts. |
---|
130 | |
---|
131 | T1 = max(abs(s(D))); % The largest scalar modulus of the bidiagonal result |
---|
132 | T2 = max(abs(v(D))); % The largest vector modulus of the bidiagonal result. |
---|
133 | |
---|
134 | if abs(T2) > abs(T1) * eps |
---|
135 | warning('Result of bidiagonalization was not accurate.') |
---|
136 | disp('Information: largest on-diagonal vector and scalar moduli were:'); |
---|
137 | disp(abs(T2)); |
---|
138 | disp(abs(T1)); |
---|
139 | end |
---|
140 | |
---|
141 | if r == 1 || c == 1 |
---|
142 | R = s(B); % The diagonal has only one element, so we can just forget the off-diagonal. |
---|
143 | else |
---|
144 | R = s(B - O); % Subtract the off-diagonal part and take the scalar part of the result. |
---|
145 | end |
---|
146 | |
---|