function [upper, lower, z] = seplambda(A, B, abstol, prtlevel)
% call: [upper, lower, z] = seplambda(A, B, abstol, prtlevel)
% bisection method to compute seplambda(A,B): the largest epsilon such that the
% epsilon-pseudospectra of A and B do not overlap, or equivalently
% the global minimum of max{sigmin(A - zI), sigmin(B - zI) over all
% z in the complex plane
% input:
% A: square complex matrix
% B: square complex matrix, not necessarily same dimension as A
% abstol: compute upper and lower bounds within abstol of each other
% prtlevel: print level (0 none, 1 print upper and lower bounds as they
% are determined by bisection, 2 also print info from initial phase)
% output:
% upper: strict upper bound on seplambda(A,B)
% lower: lower bound on seplambda(A,B)
% z: complex number certifying upper, i.e.,
% with upper = sigmin(A-zI) = sigmin(B-zI) within a tolerance
% caveat: lower bound may be invalid if an error is made
% deciding whether an eigenvalue is real or imaginary
% (tolerances are used in this decision). To get any guarantees would need to
% use specialized eigenvalue software instead of "eig".
% Upper bound is guaranteed as it is certified by z (within a tolerance).
% However, since upper bound is obtained by bisection, not optimization,
% it may not be a good upper bound.
% method:
% (1) initially: check along lines through eigenvalues to ensure epsilon is small
% enough that a component of A's epsilon-pseudospectrum can't be strictly
% inside a component of B's
% (2) bisection: at each step, do real/imaginary eigenvalue check to find whether
% eps-pseudospectra of A and B intersect for a given eps
% Written by Michael Overton, overton@cs.nyu.edu, Jan 2004, for the SVG meeting.
% Last revision Jan 2005.
% Reference: An Algorithm to Compute Sep_Lambda, Ming Gu and Michael L. Overton
% Submitted to SIAM J. Matrix Anal Appl, Jan 2005.
% http://www.cs.nyu.edu/overton/papers/pdffiles/seplambda.pdf
% M-files called: seplam_initial, seplam_bisect, linecheck, realimagcheck
if nargin < 4
prtlevel = 1;
end
if nargin < 3
abstol = 1e-6; % absolute tolerance for bisection
end
if size(A,1) ~= size(A,2) | size(B,1) ~= size(B,2)
error('A and B must be square')
end
maxnorm = max([norm(A,inf); norm(B,inf)]);
if abstol < 1e-16*maxnorm
abstol = 1e-16*maxnorm;
if prtlevel > 0
fprintf('\abstol too small, resetting it to %e', abstol)
end
end
% initial phase: checking along lines
[f,z] = seplam_initial(A, B, prtlevel);
if f < inf
upper = f;
else
upper = max(min(svd(A)), min(svd(B)));
end
if prtlevel > 0
fprintf('\nstarting bisection with upper bound equal to %e', upper)
end
lower = 0;
% bisection
[lower, upper, z] = seplam_bisect(A, B, lower, upper, z, abstol, prtlevel);
if prtlevel > 0
fprintf('\nfinished: lower = %e, upper = %e\n', lower, upper)
end