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