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