I’ve implemented the QR algorithm to find the eigenvalues and eigenvectors of a symmetric matrix using CUBLAS, which works correctly. However, it is very slow to converge. Therefore, I decided to reduce the symmetric matrix to tridiagonal form before running the QR algorithm. The eigenvalues of the original symmetric matrix and the tridiagonal matrix are the same, but how can I transform the eigenvectors of the tridiagonal matrix back to those of the original symmetric matrix? I’m hoping to do it without having to solve a system of equations, because the size of the rows for the CUBLAS solvers are limited to 4070, according to the CUBLAS Programming Guide.
This is the correct step to take.
Normally, tri-diagonalization of a symmetric matrix is done by the Householder transformation. Then, this same transformation can be used to transform the eigen vectors back to the original coordinate frame.
Here’s my octave code, which performs the two transformations. Also find references to the relevant algorithm literature therein. The missing step in the full chain of the eigen value decomposition is finding eigen values of the triangular matrix generated by triDiagonalizeSymmetric(…) function, but you said you know how to do this. As you can see from this code, the application of the Householder transform once it’s known is a fairly trivial task.
[codebox]
function [M, householder] = triDiagonalizeSymmetric(M);
% Use householder transformation to tri-diagonalize a real symmetric
% matrix. This algorithm is coded after algorithm 2.1 (page 82) of "Matrix
% Algorithms" By Gilbert W. Stewart, available for preview on Google Books.
dim = rows(M);
if columns(M) ~= dim
error(‘Matrix must be square’);
end
if 10 * norm(M) * eps < norm(M-M.')
error(‘Matrix must be symmetric’);
end
householder = zeros(dim, dim - 1);
for col = 1 : dim - 1
% Pick the relevant column from the matrix.
x = M(col + 1 : dim, col);
% Figure out Householder reflector for this column. This is modified off
% the algorithm by Stewart, in an attempt to save some flops.
u = x;
u(1) += sign(x(1)) * norm(x, 2);
% NOTE: contrary to the book, we L2-normalize u to 1 (not to sqrt(2), as in
% the book, see formula 2.19 and the expression just below it).
u *= sign(x(1)) / norm(u, 2);
H = eye(dim - col + 1);
% NOTE: since we L2-normalize u to 1, rather than to sqrt(2), as in the
% book, the Householder transformation matrix involves a factor of 2.
H(2 : dim - col + 1, 2 : dim - col + 1) -= 2 * u * u.';
M(col : dim, col : dim) = H * M(col : dim, col : dim) * H.';
householder(col + 1 : dim, col) = u;
end % for col = 1 : dim - 1
[/codebox]
[codebox]
function V = applyHouseholderReflectors(V, householder);
% Apply Householder reflectors to a set of vectors. Both are stored
% column-by-column. Reflectors are assumed to be L2-normalized to one.
if nargin ~= 2
error(‘wrong number of arguments, two expected’);
end
if rows(householder) ~= rows(V)
error(‘Dimensions of Householder matrix and vectors mismatch’);
end
% Apply reflectors one by one, in the reverse order.
for i = columns(householder): -1 : 1
reflector = householder(:, i);
V -= 2 * reflector * (reflector.’ * V);
end % of for i = columns(householder): -1 : 1
[/codebox]
Ah, okay. I thought it was possible to use the saved transformations to obtain the original eigenvectors; I just wasn’t sure how to do it. Thanks for your help.