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. |
---|