A6 Scripts

  function [Q,sep] = Invariant22(T)
% T is a real n-by-n upper quasi-triangular matrix. 
% Assume that T11 = T(1:n-2,1:n-2) is upper triangular.
% Assume that T22 = T(n-1:n,n-1:n) has complex eigenvalues.
% Q is a real n-by-2 matrix with Q'*Q = eye(2,2) such that T*Q = Q*S where S (2-by-2)
%    has the same eigenvalues as T22.
% sep = sep(T11,T22).






% P1 Test script for Invariant22
clc
randn('seed',0)

% Example 1

n = 8;
T = triu(randn(n,n));
T(n-1:n,n-1:n) = [3 4 ; -10 1];

[Q,sep] = Invariant22(T);
disp('Example 1.')
disp(' ')
disp(sprintf('sep(T11,T22)   = %10.3e',sep))
disp(sprintf('|| Q''*Q - I || = %10.3e',norm(Q'*Q - eye(2,2))))
disp(' ')
disp('If S = Q''*T*Q =')
S = Q'*T*Q;    % norm(T*Q - Q*S,'fro') is minimized by S = Q'*T*S
disp(' ')
disp(S)
disp(sprintf('then   || T*Q - Q*S ||_F = %10.3e',norm(T*Q - Q*S,'fro')))


% Example 2

n = 8;
T = triu(-ones(n,n),1)+ .025*diag(randn(n,1));
T(n-1:n,n-1:n) = .001*[3 4 ; -10 1];

[Q,sep] = Invariant22(T);
disp(' ')
disp(' ')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('Example 2.')
disp(' ')
disp(sprintf('sep(T11,T22)   = %10.3e',sep))
disp(sprintf('|| Q''*Q - I || = %10.3e',norm(Q'*Q - eye(2,2))))
disp(' ')
disp('If S = Q''*T*Q =')
S = Q'*T*Q;    % norm(T*Q - Q*S,'fro') is minimized by S = Q'*T*S
disp(' ')
disp(S)
disp(sprintf('then   || T*Q - Q*S ||_F = %10.3e',norm(T*Q - Q*S,'fro')))

% Example 3

% In the Example 2 problem, let's see how Q changes if we perturb T.
disp(' ')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('Example 3')
disp(' ')
disp('Perturb the Example 2 matrix by O(eps) and compute its real schur form:')
A = T + eps*randn(n,n);
[Q1,T1] = schur(A)
disp(' ')
disp('Note that the complex eigenvalue pair "moves up" the diagonal.')
disp('Let''s recompute the associated invariant subspace ran(QA)')
disp('and check the residual.')
[U,sepA] = Invariant22(T1(1:6,1:6));
QA = Q1(:,1:6)*U;
disp(' ')
disp('If SA = QA''*A*QA =')
SA = QA'*A*QA;    
disp(' ')
disp(SA)
disp(sprintf('then   || A*QA - QA*SA ||_F = %10.3e',norm(A*QA - QA*SA,'fro')))
disp(' ')
dist = sqrt(1 - min(svd(QA'*Q))^2);
disp(sprintf('However, dist(ran(Q),ran(QA)) = %10.3e  which is about eps/sep(T11,T22).',dist))

  function [Q,H] = SpecialHess(A)
% A is an n-by-n real matrix
% Q is orthogonal so that Q'*A*Q = H is upper Hessenberg
% and e'*Q(:,2:n) = 0 where e = ones(n,1).



% P2 Test script for SpecialHess
clc
randn('seed',0)
n = 8;
A = randn(n,n);
[Q,H] = SpecialHess(A)
disp('Residuals:')
disp(' ')
disp(sprintf('        || Q''*Q - I ||    = %10.3e',norm(Q'*Q - eye(n,n))))
disp(sprintf('        || Q''*A*Q - H ||  = %10.3e',norm(H - Q'*A*Q)))
disp(sprintf('        || e''*Q(:,2:n) || = %10.3e',norm(ones(1,n)*Q(:,2:n))))

P3-P4

     function [x,Lambda2,condLambda2] = GoogleMatrix(G,p,k)
   % G is an N-by-N 0-1 matrix represented in sparse form and 0 < p < 1.
   % k is a positive integer, the number of orthogonal iteration steps
   % Let A = A_p be the Google matrix. 
   % x is a column N-vector with positive entries so that Ax = x (approximately)
   % Lamda2 is the absolute value of A's second largest eigenvalue in 
   %    absolute value
   % Lambda2 is the condition number of Lambda2, i.e., the secant of the
   % angle between the left and right eigenvector associated with Lambda2.



  function G = MakeG(N,q)
% Constructs an N-by-N random 0-1 matrix with up to q nonzeros per
% column and returns the result (in sparse format) in G. Assumes q<N.
% For modest N, G_explicit = full(G) gives you a look at G in conventional
% form.

i = [];
j = [];
for k = 1:N
    [vals,rows] = sort(rand(N,1));  % rows is a random permutation of 1:N
    nk = floor((q+1)*rand(1,1));    % number of outlinks on page k
    if nk>0
       j = [j;k*ones(nk,1)];
       i = [i;rows(1:nk)];
    end
end
G = sparse(i,j,ones(length(i),1),N,N);




% Script P3and4

clc
disp('Example 1.')
disp(' ')
randn('seed',0)
N = 100; q = 10; p = .85;
G = MakeG(N,q);
[x,Lambda2,condLambda2] = GoogleMatrix(G,p,50);
x = x/sum(x);
[xsort,idx] = sort(-x);
disp(sprintf('Lambda2 = %10.6f  ',Lambda2))
disp(sprintf('condLambda2 = %10.3e',condLambda2))
disp(' ')
disp('Rank    Page     PageRank')
disp('--------------------------------')
for i=1:10
    disp(sprintf(' %2d  %6d    %10.6f',i,idx(i),x(idx(i))))
end

disp(' ')
disp('***********************************************************************')
disp('Example 2.')
disp(' ')
randn('seed',0)
N = 500; q = 10; p = .85;
G = MakeG(N,q);
[x,Lambda2,condLambda2] = GoogleMatrix(G,p,50);
x = x/sum(x);
[xsort,idx] = sort(-x);
disp(sprintf('Lambda2 = %10.6f  ',Lambda2))
disp(sprintf('condLambda2 = %10.3e',condLambda2))
disp(' ')
disp('Rank    Page     PageRank')
disp('--------------------------------')
for i=1:10
    disp(sprintf(' %2d  %6d    %10.6f',i,idx(i),x(idx(i))))
end