function TestP2
% Test script for Problem 2

clc
% A random example...
m = 10;
B = randn(m,m); B = B'*B ; C = randn(m,m); C = C'*C ; A0 = kron(B,C);
% Multiples of this noise matrix will be added....
E = randn(m^2,m^2); E = E + E';  E = E/norm(E);

disp('CVL Results:')
disp(' ')
disp('   A-Error   min_phiVal   norm(gradient)  EXITFLAG  Iterations  phiEvals')
disp('----------------------------------------------------------------------')
for Aerror = linspace(0,1,11)*norm(A0)
   A = A0 + Aerror*E;
   [y,z,EXITFLAG,OUTPUT] = KronMinCV(A);
   [phiVal,g] = phi([y;z],A);
   its = OUTPUT.iterations;
   phiEvals = OUTPUT.funcCount;
   fprintf('% 10.3f   %12.3e  %10.1e       %2d          %4d     %4d\n',Aerror,phiVal,norm(g),EXITFLAG,its,phiEvals)
end

disp(' ')
disp(' ')
disp('Your Results:')
disp(' ')
disp('   A-Error   min_phiVal   norm(gradient)  EXITFLAG  Iterations  phiEvals')
disp('----------------------------------------------------------------------')
for Aerror = linspace(0,1,11)*norm(A0)
   A = A0 + Aerror*E;
   [y,z,EXITFLAG,OUTPUT] = KronMin(A);
   [phiVal,g] = phi([y;z],A);
   its = OUTPUT.iterations;
   phiEvals = OUTPUT.funcCount;
   fprintf('% 10.3f   %12.3e  %10.1e       %2d          %4d     %4d\n',Aerror,phiVal,norm(g),EXITFLAG,its,phiEvals)
end

disp(' ')
disp('Gradient and phi (5pts):')
disp(' ')
disp('Correct use of fminunc (5pts):')
disp(' ')
disp('Starting value effort (5pts):')


%%% Sample Solution...............


  function [y,z,EXITFLAG,OUTPUT] = KronMinCV(A)
% A is an n-by-n symmetric matrix with n = m^2
% y and z are column m-vectors with the property that x = kron(y,z)
% minimizes (x'*A*x)/(x'*x)

 [n,n] = size(A);
 m = round(sqrt(n));
 x0 = InitialGuess(A);
 options = optimset('Display','off','GradObj','on','tolx',10^-12,'tolfun',10^-12);
 [x,phiVal,EXITFLAG,OUTPUT] = fminunc(@(x) phi(x,A),x0,options);
 y = x(1:m); z = x(m+1:2*m);
  
 
  function x = InitialGuess(A)
  %Provides an initial guess for fminunc
  [n,n] = size(A);
  m = round(sqrt(n));
  
  % Use eigs to get the minimum eigenvector of A. This will be O(m^4)...
  opts = struct('disp',0);
  [u,d] = eigs(A,1,'SA',opts);
 
  % If we assume that u = kron(y,z), then u is a length-m block vector
  % with m-vector components. Every component is a multiple of z. So
  % we can estimate a multiple of z this way...
  z = zeros(m,1);
  for i=1:m
     rows = 1+(i-1)*m:i*m;
     z = z + u(rows);
  end
  
  % Similarly for y...
  y = zeros(m,1);
  for i=1:m
     rows = i:m:n;
     y = y + u(rows);
  end
  % WLOG can assume z has unit 2-norm. This nails down both y and z...
  z = z/norm(z);
  y = y*norm(u)/norm(y);
  x = [y;z];
  
 
  function [phiVal,g] = phi(x,A)
  % phiVal is the value of the objective function and g is its
  % gradient.
  n = length(x); m = n/2;
  y = x(1:m);     yTy = y'*y;
  z = x(m+1:n);   zTz = z'*z;
  F = MakeF(A,z); Fy  = F*y;  yTFy = y'*Fy;
  G = MakeG(A,y); Gz  = G*z;  zTGz = z'*Gz;
  phiVal = (yTFy)/(yTy*zTz);
  g = (2/(yTy*zTz))*[ Fy-(yTFy/yTy)*y;... 
                     Gz-(zTGz/zTz)*z];
  
  function F = MakeF(A,z)
% A is an m^2 -by- m^2 symmetric matrix. 
% F is an m-by-m symmetric matrix with F(i,j) = z'*A(rows,cols)*z 
% where rows = 1+(i-1)*m:i*m and cols = 1+(j-1)*m:j*m
  [n,n] = size(A); m = sqrt(n); F = zeros(m,m);
  for j=1:m
     for i=1:j
        cols = 1+(j-1)*m:j*m; rows = 1+(i-1)*m:i*m;
        F(i,j) = z'*A(rows,cols)*z; F(j,i) = F(i,j);
     end
  end
  
  function G = MakeG(A,y)
% A is an m^2 -by- m^2 symmetric matrix. 
% G is an m-by-m symmetric matrix with G(i,j) = y'*A(rows,cols)*y 
% where rows = i:m:n and cols = j:m:n
  [n,n] = size(A);  m = sqrt(n); G = zeros(m,m);
  for j=1:m
     for i=1:j
        cols = j:m:n; rows = i:m:n;
        G(i,j) = y'*A(rows,cols)*y; G(j,i) = G(i,j);
     end
  end