Final Exam Solutions

Test Script:

% Tests Final Exam Solutions

randn('state',0)
clc

n = 10;
A = randn(n,n);
H = PosHess(A);
disp(sprintf('Problem 1. Test PosHess.\n\nSimilarity Check = %10.3e',norm(eig(A)-eig(H))))
disp(' ')
disp('The subdiagonal of H = ')
disp(' ')
disp(sprintf('%6.2f',diag(H,-1)))
n = 100;
A = randn(n,n);
flops(0)
H = PosHess(A);
f = flops;
disp(sprintf('\nn = %3d, flops = %1d\n',n,f))
disp('----------------------------------------------------------------------------')

% Test reflect
disp(sprintf('\nProblem 2. Test Reflect.\n'))

disp('----------')
x = 0; y = 0; [c,s] = reflect(x,y);
disp(sprintf('x = %20.16f   y = %20.16f',x,y))
disp(sprintf('c = %20.16f   s = %20.16f\n',c,s))
disp(sprintf('sqrt(c^2 + s^2) = %20.16f',sqrt(c^2+s^2)))
disp(sprintf('top = -c*x+s*y  = %20.16f',-c*x+s*y))
disp(sprintf('bot =  s*x+c*y  = %20.16f\n',s*x+c*y))

disp('----------')
x = 4; y = 3; [c,s] = reflect(x,y);
disp(sprintf('x = %20.16f   y = %20.16f',x,y))
disp(sprintf('c = %20.16f   s = %20.16f\n',c,s))
disp(sprintf('sqrt(c^2 + s^2) = %20.16f',sqrt(c^2+s^2)))
disp(sprintf('top = -c*x+s*y  = %20.16f',-c*x+s*y))
disp(sprintf('bot =  s*x+c*y  = %20.16f\n',s*x+c*y))


disp('----------')
x = -4; y = 3; [c,s] = reflect(x,y);
disp(sprintf('x = %20.16f   y = %20.16f',x,y))
disp(sprintf('c = %20.16f   s = %20.16f\n',c,s))
disp(sprintf('sqrt(c^2 + s^2) = %20.16f',sqrt(c^2+s^2)))
disp(sprintf('top = -c*x+s*y  = %20.16f',-c*x+s*y))
disp(sprintf('bot =  s*x+c*y  = %20.16f\n',s*x+c*y))

disp('----------')
x = 1; y = .0000000001 ; [c,s] = reflect(x,y);
disp(sprintf('x = %20.16f   y = %20.16f',x,y))
disp(sprintf('c = %20.16f   s = %20.16f\n',c,s))
disp(sprintf('sqrt(c^2 + s^2) = %20.16f',sqrt(c^2+s^2)))
disp(sprintf('top = -c*x+s*y  = %20.16f',-c*x+s*y))
disp(sprintf('bot =  s*x+c*y  = %20.16f\n',s*x+c*y))

disp('----------')
x = .0000000001; y = 1 ; [c,s] = reflect(x,y);
disp(sprintf('x = %20.16f   y = %20.16f',x,y))
disp(sprintf('c = %20.16f   s = %20.16f\n',c,s))
disp(sprintf('sqrt(c^2 + s^2) = %20.16f',sqrt(c^2+s^2)))
disp(sprintf('top = -c*x+s*y  = %20.16f',-c*x+s*y))
disp(sprintf('bot =  s*x+c*y  = %20.16f\n',s*x+c*y))
disp('----------------------------------------------------------------------------')


% Test OrthoHessQR

disp(sprintf('\nProblem 3. Test OrthoHessQR.\n\nAn n = 6 example.\n'))
disp('H - G_1*G_2*G_3*G_4*G_5*diag(1,1,1,1,1,-1) = ')
n = 6;
[A,H] = GenOrth(n);
[c,s] = OrthoHessQR(H);
H0 = eye(n,n);
H0(n,n) = -1;
for k=n-1:-1:1
   H0(k:k+1,:) = [-c(k) s(k);s(k) c(k)]*H0(k:k+1,:);
end
disp(' ')
disp(H0-H)
disp(' ')
disp('       c        s  ')
disp([c s ])


disp(' ')
n = 100;
[A,H] = GenOrth(n);
flops(0)
[c,s] = OrthoHessQR(H);
f = flops;
disp(sprintf('If n = %1d then flops = %1d',n,f))

disp(' ')
disp('----------------------------------------------------------------------------')
disp(' ')
disp('Problem 4. Test TriDiagProd.')
disp(' ')

n = 6;
M = magic(n);
Y = M(:,1:3)
Z = M(:,4:6)
X = TriDiagProd(Y,Z)
B = diag(Y(1:n-1,1),-1) + diag(Y(:,2)) + diag(Y(1:n-1,3),1);
C = diag(Z(1:n-1,1),-1) + diag(Z(:,2)) + diag(Z(1:n-1,3),1);
A = B*C;
disp('Exact X = ')
disp(' ')
disp([[diag(A,-2);0;0] [diag(A,-1);0] diag(A) [diag(A,1);0] [diag(A,2);0;0]])

n = 100;
Y = randn(n,3);
Z = randn(n,3);
flops(0)
X = TriDiagProd(Y,Z);
f = flops;
disp(sprintf('\n\nIf n = %1d, then flops = %1d\n',n,f))
disp('----------------------------------------------------------------------------')

disp(' ')
disp('Problem 5. Test OrthoEig.')
disp(' ')


n = 8;
[A,H] = GenOrth(n);
[c,s] = OrthoEig(A);
[z,idx] = sort(real(c));
ev1 = zeros(n,1);
ev1(1:2:n) = c(idx)+ sqrt(-1)*s(idx);
ev1(2:2:n) = c(idx)- sqrt(-1)*s(idx);

ev2 = eig(A);
[z,jdx] = sort(real(ev2));
ev2 = ev2(jdx);


format long
disp(' ')
disp('An n = 8 example.')
disp(' ')
disp('Eigenvalues via OrthoEig:')
disp(' ')
disp(ev1)
disp(' ')

disp('Eigenvalues via eig:')
disp(' ')
disp(ev2)
disp(' ')

format short

n = 100;
[A,H] = GenOrth(n);
flops(0)
[c,s] = OrthoEig(A);
f = flops;
flops(0)
H = hess(A);
fHess = flops;
disp(sprintf('\nFor an n = %1d example flops = %1d (excluding Hess reduction)',n,f-fHess))

Problem 1 (5 points)

  function H = PosHess(A)
% A is an n-by-n real matrix.
% H is an n-by-n upper Hessenberg matrix with nonnegative
% subdiagonal elements that is orthogonally similar to A. 

[n,n] = size(A);
H = hess(A);
for k=1:n-1
   if H(k+1,k) < 0
      % Multiple row k+1 and column k+1 by -1.
      % This is a similarity transformation.
      H(:,k+1) = -H(:,k+1);
      H(k+1,:) = -H(k+1,:);
   end
end

Problem 2 (5 points)

  function [c,s] = reflect(x,y)
% x and y are real scalars and y>=0.
% c and s are a cosine-sine pair so that (a) s>=0
% and (b) [-c s; s c]*[x;y] = [sqrt(x^2 + y^2);0].

if y==0
   if x>=0
      c = -1; s = 0;
   else
      c = 1; s = 0;
   end
else
   if y>=abs(x)
      tau = -x/y; s = 1/sqrt(1+tau^2); c = s*tau;
   else
      tau = -y/x; c = 1/sqrt(1+tau^2); s = c*tau;
      if s<0
         c = -c; s = -s;
      end
   end
end

Problem 3 (10 points)

  function [c,s] = OrthoHessQR(H)
% H is an orthogonal, n-by-n, upper Hessenberg matrix with
% nonnegative subdiagonal entries. Assume that n is even
% and that H has no real eigenvalues.
% c and s are column (n-1)-vectors so that
%
%           H = (G_1*...G_n-1)diag(1,...,1,-1)
%
% where G_k is the identity except that G_k(k:k+1,k:k+1) = [-c(k) s(k);s(k) c(k)]
% where s(k)>=0 and c(k)^2+s(k)^2 = 1 for k=1:n-1.

[n,n] = size(H);
c = zeros(n-1,1);
s = zeros(n-1,1);

for k=1:n-1
   
   % Determine a cosine-sine pair (c(k),s(k)) so that 
   % [-c(k) s(k);s(k) c(k) ]*H(k:k+1,k) = [1;0]
   
   [c(k),s(k)] = reflect(H(k,k),H(k+1,k));
   
   % H(k:k+1,k:n) = [-c(k) s(k);s(k) c(k) ]*H(k:k+1,k:n) 
   % Notice that the first row of the update is the first row of eye(n-k+1,n-k+1)
   
   H(k+1,k:n) = [s(k) c(k)]*H(k:k+1,k:n);

end

Problem 4 (15 points)

  function X = TriDiagProd(Y,Z)
% Suppose Y and Z are n-by-3 matrices and define
% the tridiagonal matrices B and C by
%
%       B = diag(Y(1:n-1,1),-1) + diag(Y(:,2)) + diag(Y(1:n-1,3),1)
%       C = diag(Z(1:n-1,1),-1) + diag(Z(:,2)) + diag(Z(1:n-1,3),1).
%
% X is a n-by-5 matrix with the property that if
%
%   A = diag(X(1:n-2,1),-2) + diag(X(1:n-1,2),-1) + diag(X(:,3)) +
%       diag(X(1:n-1,4), 1) + diag(X(1:n-2,5), 2)
%
% then A = B*C.

[n,m] = size(Y);

X = zeros(n,5);
X(1:n-2,1) = Y(2:n-1,1).*Z(1:n-2,1);
X(1:n-1,2) = Y(2:n,2).*Z(1:n-1,1)   + Y(1:n-1,1).*Z(1:n-1,2);
X(1:n,3) = [0;Y(1:n-1,1).*Z(1:n-1,3)] + Y(:,2).*Z(:,2) + [Y(1:n-1,3).*Z(1:n-1,1);0];
X(1:n-1,4) = Y(1:n-1,2).*Z(1:n-1,3) + Y(1:n-1,3).*Z(2:n,2);
X(1:n-2,5) = Y(1:n-2,3).*Z(2:n-1,3);

Problem 5 (15 points)

  function [lambda,mu] = OrthoEig(A)
% H is an n-by-n unreduced upper Hessenberg matrix that is orthogonal.
% Assume that n is even and that 1 and -1 are not eigenvalues.
% lambda and mu are column m-vectors where m = n/2 and the eigenvalues of
% A are lambda(k)(plus or minus) sqrt(-1)*mu(k) for k=1:m

[n,n] = size(A);
m = n/2;
H = PosHess(A);
[c_phi,s_phi] = OrthoHessQR(H);

% Get the cosines and sines of the half-angles:
   
c_phi_half = sqrt((1+c_phi)/2);
s_phi_half = sqrt((1-c_phi)/2);
   
% Represent  Q_o in Y
   
% Q_o = zeros(n,n);
% for k=1:2:n-1
%    Q_o(k:k+1,k:k+1) = [-c_phi_half(k) s_phi_half(k);s_phi_half(k) c_phi_half(k)];
% end

Y = zeros(n,3);
Y(1:2:n-1,1) =  s_phi_half(1:2:n-1);
Y(1:2:n-1,2) = -c_phi_half(1:2:n-1);
Y(2:2:n,2)   =  c_phi_half(1:2:n-1);
Y(1:2:n-1,3) =  s_phi_half(1:2:n-1);

% Represent Q_e in Z

% Q_e = zeros(n,n);
% Q_e(1,1) = 1;
% for k=2:2:n-1
%    Q_e(k:k+1,k:k+1) = [-c_phi_half(k) s_phi_half(k);s_phi_half(k) c_phi_half(k)];
% end
% Q_e(n,n) = -1;


Z = zeros(n,3);
Z(2:2:n-1,1) =  s_phi_half(2:2:n-1);
Z(1,2)       =  1;
Z(2:2:n-1,2) = -c_phi_half(2:2:n-1);
Z(3:2:n-1,2) =  c_phi_half(2:2:n-1);
Z(n,2)       = -1;
Z(2:2:n-1,3) =  s_phi_half(2:2:n-1);

% Represent Q = Q_o*Q_e in X.

X = TriDiagProd(Y,Z);
   
% Define D by
%
% D = zeros(n,n);
% for k=1:n
%     D(k,k) = (-1)^k;
% end

% Note: D = D_o = - D_e .

% Consider the matrix C1 = (D*Q - Q*D)/2.
% Apart from pre and post multiplication by diagonal matrices
% of plus and minus ones, C1 = diag(diag(Q,-1),-1) + diag(Q,1),1).
%(We can ignore these diagonal scalings because they do not effect
% the singular values.

% Set up the Bidiagonal BC defined by
%         C1 = (D*Q - Q*D)/2
%         C1 = C1([2:2:n 1:2:n],[1:2:n 2:2:n]);
%         BC = C1(1:m,1:m);

BC = diag(X(1:2:n-1,2)) + diag(X(2:2:n-1,4),1);
chalf = svd(BC);


% Consider the matrix S1 = (D*Q + Q*D)/2. Apart from pre and post 
% multiplication by diagonal matricesof plus and minus ones, 
% S1 = diag(diag(Q,-2),-2) + diag(diag(Q)) + diag(diag(Q,2),2)
%(We can ignore these diagonal scalings because they do not effect
% the singular values.

% Set up the bidiagonal BS defined by
%         S1 = (D*Q + Q*D)/2;
%         S1 = S1([1:2:n 2:2:n],[1:2:n 2:2:n]);
%         BS = S1(1:m,1:m)

BS = diag(X(1:2:n-1,3)) + diag(X(1:2:n-2,5),1); 
shalf = svd(BS);

% Since the cosines are ordered big to little we must order the
% sines from little to big.

shalf = shalf(m:-1:1);

% Use double angle formulae to get the cosines and sines

lambda = chalf.*chalf - shalf.*shalf;
mu = 2*chalf.*shalf;