# How to plot the $\varepsilon$-pseudospectrum of a matrix

Let $A\in\mathbb{C}^{n\times n}$ and $\varepsilon>0$. The $\varepsilon$-pseudospectrum of $A$ is the set $$\sigma_{\varepsilon}(A)=\left\{z\in\mathbb{C} \middle| \|(A-zI)^{-1}\|\geq\frac{1}{\varepsilon}\right\}.$$ If, for some $z\in\mathbb{C}$, $A-zI$ is not invertible (that is $z$ is an eigenvalue of $A$), we put $\|(A-zI)^{-1}\|:=\infty$. So the spectrum is always contained in the pseudospectrum. The norm here is the operator norm or induced norm of a matrix (see wiki).

I would like to plot the level sets $$\left\{z\in\mathbb{C} \middle| \|(A-zI)^{-1}\|=\frac{1}{\varepsilon}\right\}$$ for various values of $\varepsilon>0$ and given $A\in\mathbb{C}^{n\times n}$ of large order $n$.

There is a package called EigTool for MATLAB which works nicely. However, I did not find anything similar for Mathematica. Is there any (hopefully straightforward) way to effectively plot the level sets in Mathematica?

I still haven't gotten around to implementing the Arnoldi-based algorithm of Toh and Trefethen, so what follows is straightforward, but not very efficient code for visualizing $\varepsilon$-pseudospectra:

mat = ExampleData[{"Matrix", "WEST0067"}];

With[{c = Norm[mat], n = Length[mat]},
DensityPlot[-c Log10[First[
SingularValueList[mat - (x + I y) IdentityMatrix[n, SparseArray],
-1, Method -> "Arnoldi", Tolerance -> 0]]],
{x, -2, 2}, {y, -2, 2}, ColorFunction -> "SolarColors"]]


Compare with the "spectral portrait" here.

It is quite easy to modify the code above to use ContourPlot[] instead of DensityPlot[], since all that matters is the plotting of the (logarithm of the) greatest singular value of the resolvent matrix $(\mathbf A-\lambda\mathbf I)^{-1}$ (equivalently, the least singular value of its inverse). Here's the example used in the EigTool download page:

(* Grcar matrix, https://arxiv.org/abs/1203.2390 *)
grcar[r : _Integer?Positive : 3, n_Integer?Positive] :=
SparseArray[{{j_, k_} /; j == k + 1 :> -1, {j_, k_} /; 0 <= k - j <= r :> 1}, {n, n}]

mat = grcar[32]; eigs = Eigenvalues[N[mat]];

With[{c = Norm[mat], n = Length[mat]},
ContourPlot[-c Log10[First[
SingularValueList[mat - SparseArray[Band[{1, 1}] -> x + I y, {n, n}],
-1, Method -> "Arnoldi", Tolerance -> 0]]],
{x, -2, 4}, {y, -4, 4}, AspectRatio -> Automatic,
Contours -> Range[1, 6, 1/2], ContourShading -> None,
ContourStyle -> Black,
Epilog -> {Directive[Black, AbsolutePointSize[5]], Point[ReIm /@ eigs]}]]


Coloring and other sundry styling is left to the interested reader.

• Notes: 1. The matrix norm (c in the code above) is merely a scaling factor, and can be left out if desired. 2. Jeff Bryant did some previous work here. – J. M. will be back soon Jan 28 '17 at 18:41