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