[37] | 1 | function [P, B] = tridiagonalize(A) |
---|
| 2 | % Tridiagonalize Hermitian matrix A, such that P * A * P' = B and P' * B * P = A. |
---|
| 3 | % B is real, P is unitary, and B has the same eigenvalues as A. |
---|
| 4 | |
---|
| 5 | % Copyright © 2005 Stephen J. Sangwine and Nicolas Le Bihan. |
---|
| 6 | % See the file : Copyright.m for further details. |
---|
| 7 | |
---|
| 8 | error(nargchk(1, 1, nargin)), error(nargoutchk(2, 2, nargout)) |
---|
| 9 | |
---|
| 10 | % NB: This is a reference implementation. It uses explicit Householder |
---|
| 11 | % matrices. Golub and van Loan, 'Matrix Computations', 2e, 1989, section |
---|
| 12 | % 5.1.4 discusses efficient calculation of Householder transformations. |
---|
| 13 | % This has been tried for bidiagonalization (see bidiagonalize2.m), but |
---|
| 14 | % found to be slower than the use of explicit matrices as used here. This |
---|
| 15 | % requires further study. |
---|
| 16 | |
---|
| 17 | [r, c] = size(A); |
---|
| 18 | |
---|
| 19 | if r ~= c |
---|
| 20 | error('Cannot tridiagonalize a non-square matrix'); |
---|
| 21 | end |
---|
| 22 | |
---|
| 23 | if r < 2 |
---|
| 24 | error('Cannot tridiagonalize a matrix smaller than 2 by 2.'); |
---|
| 25 | end |
---|
| 26 | |
---|
| 27 | if ~ishermitian(A) |
---|
| 28 | error('Matrix to be tridiagonalized is not (accurately) Hermitian'); |
---|
| 29 | end |
---|
| 30 | |
---|
| 31 | [P, B] = internal_tridiagonalizer(A); |
---|
| 32 | |
---|
| 33 | B = check(B); % Verify the result and convert to exactly tridiagonal real form. |
---|
| 34 | |
---|
| 35 | % ----------------------------------------------------------------------------- |
---|
| 36 | |
---|
| 37 | function [P, B] = internal_tridiagonalizer(A) |
---|
| 38 | |
---|
| 39 | [r, c] = size(A); |
---|
| 40 | |
---|
| 41 | if r ~= c |
---|
| 42 | error('Internal tridiagonalize error - non-square matrix'); |
---|
| 43 | end |
---|
| 44 | |
---|
| 45 | % Compute and apply a Householder transformation to the first row and column of A (these are |
---|
| 46 | % conjugates of each other because A is Hermitian). We omit the first element of the first |
---|
| 47 | % column in computing the Householder matrix, because it is already real (A is Hermitian). |
---|
| 48 | % We apply P on both sides of A so that the first row and column are nulled out (apart from |
---|
| 49 | % the first two elements in each case). Note that if B is 2 by 2, this is all we have to do. |
---|
| 50 | |
---|
| 51 | P = quaternion(eye(r)); |
---|
| 52 | P(2 : end, 2 : end) = householder_matrix(A(2 : end, 1), eye(r - 1, 1)); |
---|
| 53 | B = P * A * P'; |
---|
| 54 | |
---|
| 55 | % Apply the algorithm recursively to sub-matrices of B, except when the sub-matrix is 2 by 2. |
---|
| 56 | |
---|
| 57 | if r > 2 |
---|
| 58 | Q = quaternion(eye(r)); |
---|
| 59 | [Q(2 : end, 2 : end), B(2 : end, 2 : end)] = internal_tridiagonalizer(B(2 : end, 2 : end)); |
---|
| 60 | P = Q * P; |
---|
| 61 | end |
---|
| 62 | |
---|
| 63 | % --------------------------------------------------------------------------------------------- |
---|
| 64 | |
---|
| 65 | function R = check(B); |
---|
| 66 | |
---|
| 67 | % Verify results, and convert the result to exactly tridiagonal form with no vector part. |
---|
| 68 | |
---|
| 69 | r = length(B); % We know B is square, so we don't need to check this again. |
---|
| 70 | |
---|
| 71 | D = [diag(B); diag(B, +1); diag(B, -1)]; % Extract the three diagonals. |
---|
| 72 | O = tril(B, -2) + triu(B, +2); % Extract the off-diagonal part. |
---|
| 73 | |
---|
| 74 | % Find the largest on and off tridiagonal elements. We use abs twice to |
---|
| 75 | % allow for the case where B is complex (this occurs when A is a complexified |
---|
| 76 | % quaternion matrix). |
---|
| 77 | |
---|
| 78 | T1 = max(max(abs(abs(O)))); % Find the largest off tridiagonal element. |
---|
| 79 | T2 = max(abs(abs(D))); % Find the largest tridiagonal element. |
---|
| 80 | |
---|
| 81 | % NB T2 and/or T1 could be exactly zero (example, if A was zero, or an identity matrix). |
---|
| 82 | % Therefore we do not divide one by the other, but instead multiply by a tolerance. |
---|
| 83 | |
---|
| 84 | tolerance = eps .* 1.0e3; % This is empirically determined. |
---|
| 85 | |
---|
| 86 | if T1 > T2 .* tolerance |
---|
| 87 | warning('Result of tridiagonalization was not accurately tridiagonal.'); |
---|
| 88 | disp('Information: largest on- and off-tridiagonal moduli were:'); |
---|
| 89 | disp(T2); |
---|
| 90 | disp(T1); |
---|
| 91 | end |
---|
| 92 | |
---|
| 93 | % Verify that the diagonal elements have neglible vector parts. |
---|
| 94 | |
---|
| 95 | T1 = max(abs(abs(s(D)))); % The largest scalar modulus of the tridiagonal result |
---|
| 96 | T2 = max(abs(abs(v(D)))); % The largest vector modulus of the tridiagonal result. |
---|
| 97 | |
---|
| 98 | if T2 > T1 .* tolerance |
---|
| 99 | warning('Result of tridiagonalization was not accurately scalar.') |
---|
| 100 | disp('Information: largest on-tridiagonal vector and scalar moduli were:'); |
---|
| 101 | disp(T2); |
---|
| 102 | disp(T1); |
---|
| 103 | end |
---|
| 104 | |
---|
| 105 | R = s(B - O); % Subtract the off-tridiagonal part and take the scalar part of the result. |
---|