  function P3Grade()
% Checks out BestNby2...

clc
n = 10;
noise = .001;
A = randn(n,n);
[Q,T] = schur(A);
Q0 = Q(:,1:2);
B0 = T(1:2,1:2);
B = B0 + noise*randn(2,2);

XCVL = BestNby2CVL(A,B);
[QCVL,R] = qr(XCVL,0);
s = svd(QCVL'*Q0);
distCVL = sqrt(1 - s(2)^2);
errQCVL = norm(QCVL'*QCVL - eye(2,2));
fprintf('dist(ran(Q0),ran(XCVL))   = %10.3e\n||QCVL''QCVL - I||         = %10.3e\n',distCVL,errQCVL)
disp(' ')
XYOU = BestNby2(A,B);
[QYOU,R] = qr(XYOU,0);
s = svd(QYOU'*QCVL);
distYOU = sqrt(1 - s(2)^2);
errQYOU = norm(QYOU'*QYOU - eye(2,2));
fprintf('dist(ran(XCVL),ran(XYOU)) = %10.3e\n||QYOU''QYOU - I||         = %10.3e\n',distYOU,errQYOU)


  function X = BestNby2CVL(A,B)
% A is nxn, Q is nx2, Q'*Q = eye(2,2).
% B is a 2x2 matrix that minimizes norm(AQ-QB,'fro').
[n,n] = size(A);
I = eye(n,n);
% norm(AX-XB,'fro')^2 = 
%                       || A*X(:,1) - B(1,1)*X(:,1) - B(2,1)*X(:,2) ||^2 +
%                       || A*X(:,2) - B(1,2)*X(:,1) - B(2,2)*X(:,2) ||^2 
%                     = || M*v ||
%  where 

M = [A-B(1,1)*I  -B(2,1)*I ; -B(1,2)*I  A-B(2,2)*I];

% and v = [X(:,1);X(:,2)]. To minimize || Mv|| over all unit vectors we 
% compute 
[~,S,V] = svd(M);
X = [V(1:n,2*n)  V(n+1:2*n,2*n)];
