  function P3Grade()
% Checks out TwoMatTrace
% 10 points
% There is a closed form solution that involves the schur decomp of A2-A1.
% Around 5/10 points for methods that rely on interesting ideas but which
% are incorrect.

clc

% Example 1
n = 1000;
A1 = randn(n,n); A1 = A1+A1';
A2 = randn(n,n); A2 = A2+A2';
[QCVL,rCVL,sigmaCVL] = TwoMatTraceCVL(A1,A2);
[QYOU,rYOU,sigmaYOU] = TwoMatTrace(A1,A2);
err = norm(QYOU'*QYOU - eye(n,n));
fprintf('||QYOU''*QYOU - I || = %10.3e\n',err)
d1 = sqrt(1-min(svd(QCVL(:,1:rCVL)'*QYOU(:,1:rYOU)))^2);
d2 = sqrt(1-min(svd(QCVL(:,rCVL+1:n)'*QYOU(:,rYOU+1:n)))^2);
fprintf('n = %1d\nrYOU = %1d\nrCVL = %1d\nsigmaYOU = %15.12f\nsigmaCVL = %15.12f\n',...
                              n,rYOU,rCVL,sigmaYOU,sigmaCVL)
fprintf('d1 = %10.6f\n',d1)
fprintf('d2 = %10.6f\n',d2)

% Example 2
n = 100;
A1 = randn(n,n); A1 = A1*A1';
A2 = randn(n,n); A2 = -A2*A2';
[QCVL,rCVL,sigmaCVL] = TwoMatTraceCVL(A1,A2);
[QYOU,rYOU,sigmaYOU] = TwoMatTrace(A1,A2);
fprintf('\n\nn = %1d\nrYOU = %1d\nrCVL = %1d\nsigmaYOU = %15.12f\nsigmaCVL = %15.12f\n',...
                              n,rYOU,rCVL,sigmaYOU,sigmaCVL)
                          
% Example 3
n = 100;
A1 = randn(n,n); A1 = -A1*A1';
A2 = randn(n,n); A2 = A2*A2';
[QCVL,rCVL,sigmaCVL] = TwoMatTraceCVL(A1,A2);
[QYOU,rYOU,sigmaYOU] = TwoMatTrace(A1,A2);
fprintf('\n\nn = %1d\nrYOU = %1d\nrCVL = %1d\nsigmaYOU = %15.12f\nsigmaCVL = %15.12f\n',...
                              n,rYOU,rCVL,sigmaYOU,sigmaCVL)




function [Q,r,sigma] = TwoMatTraceCVL(A1,A2)
% A1 and A2 are nxn and symmetric
% Q is orthogonal and r (0<=r<=n) maximize
%   phi(Q,r) = tr(Q(:,1:r)'*A1Q(:,1:r)) + tr(Q(:,r+1:n)'*A1Q(:,r+1:n))
% and sigma is the maximized value.

% For a given r we write Q1 = Q(:,1:r) and Q2 = Q(:r+1:n) and get
%   phi(Q,r) = tr(Q1'*A1*Q1) + tr(Q2'*A2*Q2) 
%            = tr(A1*Q1*Q1') + tr(A2*Q2*Q2')
%            = tr(A1*Q1*Q1') + tr(A2*(I - Q1*Q1')
%            = tr(A1*Q1*Q1')) - tr(A2*Q1*Q1') + tr(A2)
%            = tr((A1-A2)Q1*Q1') + tr(A2)
%           <= lambda_1(A1-A2) + ... + lambda_r(A1-A2) + tr(A2)

% Thus, the columns of Q1 should span the space associated with the r
% largest eigenvalues of A1-A2. 
[Q,D] = schur(A1-A2);
d = diag(D);
n = length(d);
% Sort from big to little to resonate with the above analyis..
[d,idx] = sort(d,'descend');
Q = Q(:,idx);
% For a given r, this Q matrix is optimal since the columns of Q(:,1:r) 
% span the subspace associated with the r largest eigenvalues of A1-A2.
% Now consider the "best r" given Q. From the inequality
%
%      phi(Q,r) = lambda_1(A1-A2) + ... + lambda_r(A1-A2) + tr(A2)   (*)
%
% we see that the number of positive eigenvalues is important
nPos = sum(d>0);
if nPos == 0
    % The eigenvalues in the sum (*) are all negative and so
    % they detract from  tr(A2). Do not include them.
    r = 0;
    sigma = sum(diag(A2));   % = tr(A2)
elseif nPos == n
    % The eigenvalues in the sum (*) are all positive and so
    % they add to  tr(A2). Include them all.
    r = n;
    sigma = sum(diag(A1));   % = tr(A1)
else
    % Just include the positive eigenvalues...
    r = nPos;
    sigma = sum(d(1:r)) + sum(diag(A2));
end

    



