  function ShowHessOfSym()
% Illustrates the decomposition U'AU = T where A is symmetric, U is
% orthogonal, and T is tridiagonal. 
% And it show how to compute an eigenvector of a tridiagonal matrix if you
% know the eigenvalue.

clc
% Generate a random symmetric matrix
n = 5;
A = randn(n,n); A = A+A'
% Compute an orthogonal U so that U'AU = T is tridiagonal....
format short
[U,T] = hess(A)
% Let's compute the eigenvalues of T:
d = eig(T);
% Let's compute the eigenvectors by a back substitution process.
% We assume that no subdiagonal entry of T is zero. (If T(k+1,k) = 0,
% T is block diagonal with diagonal blocks T(1:k,1:k) and T(k+1:n,k+1:n).
% Suppose lambda is an eigenvalue of T. We want to find a nonzero x
% so  (T - lambda*I)x = 0.
a = diag(T);
b = diag(T,-1);
% Set x(n) = 1 and visit the ith component of (T - lambdaI)x = 0:
%    b(i-1)x(i-1) + (a(i)-lambda)x(i) b(i)x(i+1) = 0
% i.e.,
%          x(i-1) = (-(a(i)-lambda)x(i) - b(i)x(i+1))/b(i-1)
% Evaluating this for i=n:-1:2 gives us x(1:n-1).

% Let' use this idea to assemble all the eigenvectors of T into X:
X = zeros(n,n);
for k=1:n
    x = zeros(n,1); x(n) = 1;
    lambda = d(k);
    for i=n:-1:2
        if i==n
           x(i-1) = (-(a(i)-lambda)*x(i))/b(i-1);
        else
           x(i-1) = (-(a(i)-lambda)*x(i) - b(i)*x(i+1))/b(i-1);
        end
    end
    X(:,k) = x/norm(x);
end
X = X
D = diag(d)
fprintf('|| X''*X - I ||   = %10.3e\n',norm(X'*X - eye(n,n)))
fprintf('|| X''*T*X - D || = %10.3e\n',norm(X'*T*X - D))

