A4 Solutions

  function m = Sturm(a,b,lambda)
% a is a column n-vector
% b is a column (n-1)-vector with non-zero entries
% lambda is a scalar
% m is the number of eigenvalues of T that are less than or equal to lambda where
% T is the symmetric tridiagonal matrix with T(i,i) = a(i) for i=1:n and 
% T(i,i+1) = b(i) for i=1:n-1
  
pkm2 = 1;
pkm1 = a(1) - lambda;
if pkm2*pkm1<=0
    m = 1;
else
    m = 0;
end
n = length(a);
for k=2:n
    pk = (a(k) - lambda)*pkm1 - b(k-1)^2*pkm2;
    
    % If pk and pkm1 have opposite sign, increment m...
    
    % Method 1
    % if pk*pkm1 <=0 , m = m + 1; end
    
    % Method 2
    % if sign(pk) ~= sign(pkm1) , m = m + 1; end
    
    % Method 3
     if sign(pk)*sign(pkm1) <=0 , m = m + 1; end
    
    pkm2 = pkm1;
    pkm1 = pk;
end
 
% For even modest problems, the pk can get huge. For example, an
% 2m-by-2m matrix with eigenvalues -m, -m+1 , ..., -1 , 1 2 ,..., m
% has the property that p_n(0) = (m!)^2. The floating point work
% cannot represent such quantities. Thus, if pk = 1000! (for example)
% pk is assigned the value of inf and this can complicate the sign checking.

% one point for t he explanation

% You should also note that if pk is supposed to be exactly 0 that will
% not generally happen. Altho it is interesting to explore the imapact
% of more liberal definitions of zero, i.e., |pk| <= EPS *norm(a) there is
% marginal value in oncorporating such a standard. "Off by 1" in the value
% returned by sturm when mu is an "exact" eigenvalue is perfectly OK given 
% fuzzy data and inexact arithetic

    function [mu,SturmCalls] = Locate(a,b,L,R,tol,m)
  % a is a column n-vector
  % b is a column (n-1)-vector with non-zero entries.
  % Let T be the symmetric tridiagonal matrix 
  % with T(i,i) = a(i) for i=1:n and T(i,i+1) = b(i) for i=1:n-1.
  % tol is a positive scalar.
  % mu is an approximation to the mth smallest eigenvalue of T
  % SturmCalls is the number of calls to Sturm that are required.
  % Assume that the m-th smallest eigenvalue of T is in [L,R]
  
  LeftVal  = Sturm(a,b,L);    % Should be <= m-1
  RightVal = Sturm(a,b,R);    % Should be >= m
  SturmCalls = 2;
 
  % Invariant: LeftVal <=m-1   RightVal > m-1
  while (R-L > tol)
      mid = L + (R-L)/2;
      midVal = Sturm(a,b,mid);
      SturmCalls = SturmCalls+1;
      if midVal <=m-1
         L = mid;
         LeftVal = midVal;
      else
         R = mid;
         RightVal = midVal;
      end
  end
  mu = (L+R)/2;
  
  
  
      function [mu,SturmCalls] = MinEig(a,b,tol)
  % a is a column n-vector
  % b is a column (n-1)-vector with non-zero entries.
  % tol is a positive scalar.
  % mu is an approximation to lambda_min, the eigenvalue of T that has 
  % minimum absolute value where T is the symmetric tridiagonal matrix 
  % with T(i,i) = a(i) for i=1:n and T(i,i+1) = b(i) for i=1:n-1.
  % SturmCalls is the number of calls to Sturm that are required.
  % It should be the case that |mu - lambda_min| <= tol.
  
  n = length(a);
  GershRadii = [0 ; abs(b)] + [abs(b); 0];
  L = min(a-GershRadii);
  R = max(a+GershRadii);  
  % All eigenvalues in union of intervals [a(k)-GershRadii(k),a(k)+GershRadii(k)], 
    
  m = Sturm(a,b,0);  % The number of eigenvalues <= 0
  if m==0
      [mu,sCallsL] = Locate(a,b,0,R,tol,n);
      SturmCalls = sCallsL + 1;
  elseif m==n
      [mu,sCallsR] = Locate(a,b,L,0,tol,1);
      SturmCalls = sCallsR + 1;
  else
      % Look for eigenvalue m and m+1, which straddle the origin
      
      
  
  % Let's find the smallest negative eigenvalue..

     [muL,sCallsL] = Locate(a,b,L,0,tol,m);
  
  % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      
  % We could just do the same thing to the right of the origin...
  
     %   [muR,sCallsR] = Locate(a,b,0,R,tol,m+1);
     %   if abs(muL)<= abs(muR)
     %      mu = muL;
     %   else
     %      mu = muR;
     %   end
     %   SturmCalls = 1 + sCallsL + sCallsR;

  % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
  % But why not check to see how many eigenvalues are to the left of -muL?
  
     mplus = Sturm(a,b,-muL);
     if mplus == m
        mu = muL;
        SturmCalls = 2 + sCallsL;
     else
        % We know that the smallest positive eigenvalue is between 0 and -muL
        [muR,sCallsR] = Locate(a,b,0,-muL,tol,m+1);
        if abs(muL)<= abs(muR)
           mu = muL;
        else
           mu = muR;
        end
        SturmCalls = 2 + sCallsL + sCallsR;
     end 
end

% - 3 points for not exploiting the fact that once you know muL, you can cut
% some corners getting muL. On random examples, half the time |muL| < |muR| and
% you only need one call to Locate Sturm to discover that. And when that is not the
% case, you go after muR with a smaller initial interval, i.e., [0,-muL] instead of
% [0,R]

    function [Q,D,mu,x] = SchurPlus(A,v)
  % A is an n-by-n symmetric positive definite matrix.
  % v is a column n-vector
  % Q'*A*Q = D is the Schur decomposition of A.
  % mu is the largest eigenvalue of A + v*v' and x is an associated eigenvector, 
  % normalized so  norm(x,2) = 1. Assume that the largest eigenvalue of A and
  % A + v*v' are unique
  
  n = length(v);
  [Q,D] = schur(A);
  
  % (A + v*v')x = mu* x  transforms to (D + w*w')y = mu*y  where
  % D = Q'*A*Q, w = Q'*v, and y = Q'*x
  % So go after the largest eigenvalue/eigenvector of D + w*w' using power
  % method.
  
  w = Q'*v;
  d = diag(D);
  
  % Apply power method to get mu and x.
  
  %Initial guess (1 point for something better than a random guess)
  
  [mu,j] = max(d+w.*w);  % mu = largest diagonal value of D + w*w'
  y = zeros(n,1); y(j) = 1;
  
  tol = 10^(-14)*norm(A + v*v','inf');
  
  z = d.*y + (w'*y)*w;   % z = (D + w*w')*(current y)
  r = z - mu*y;          % r = z - (current mu)*(current y)
  while norm(r) > tol
     y = z/norm(z);         % New y is the normalized old z
     z = d.*y + (w'*y)*w;   % z = (D + w*w')*(current y)
     mu = y'*z;
     r = z - mu*y;          % r = z - (current mu)*(current y)
  end
  
  % About 11n flops per iteration
  % -2 if O(n^2) per iteration
  
  x = Q*y;