# Pseudospectra in R implementation

Pseudospectra are typically computed by establishing a grid with $$N$$ points on a region of the complex plane, computing the resolvent norm $$||(zI − A)^{−1}||$$ at each grid point z, and visualizing with a contour plotter. Letting $$\sigma_{max}(·)$$ and $$\sigma_{min}(·)$$ denote the largest and smallest singular values of an input matrix, respectively, we remark that the resolvent norm satisfies

$$||(zI − A)^{−1}||_{2} =\sigma_{max}( (zI − A)^{−1})= \frac{1}{\sigma_{min}(zI − A)}$$ Thus, one could naively compute pseudospectra by computing the SVD of $$zI−A$$ for each grid point $$z$$ and reporting the reciprocal of the smallest singular value. Doing so in R I have a problem for implement it.


B = matrix(c(6, 6,
-1, 1),2,2,byrow = TRUE);B
I = diag(2)
resolution = 10^(-5)
z = seq(-4,0,by=resolution);z
min(svd(z*I - B)$d)  Any help? ## 1 Answer If you want to compute the minimum singular value for each $$z$$, you need to compute the SVD for each $$z$$, with a loop. For instance: sapply( z, function(u) min(svd(u*I - B)$d) )


Unless you know the eigenvalues are real (as in your example), you may prefer a grid in the complex plane.

xs <- seq(2.5,4.5,.01)
ys <- seq(-.5,.5,.01)
zs <- expand.grid( xs, ys )
zs <- apply( zs, 1, function(u) u[1] + 1i * u[2] )
S <- sapply( zs, function(z) min(svd(z*I - B)$d) ) S <- matrix( S, nr = length(xs), nc = length(ys) ) image( xs, ys, S ) contour( xs, ys, S, add = TRUE )  • If I don't know the eigen values if they are real?Maybe$B+B^{T}/2\$ – Tsiantakis Mar 31 at 5:45
• > image( xs, ys ) Error in image.default(xs, ys) : argument must be matrix-like > contour( xs, ys, S, add = TRUE ) Error in contour.default(xs, ys, S, add = TRUE) : plot.new has not been called yet – Tsiantakis Mar 31 at 5:49
• @Tsiantakis: this is fixed -- I had forgotten one argument. – VDZ Mar 31 at 7:43
• If you do not know the eigenvalues are real, and if you really want that matrix, you can use a grid in the complex plane (as in my example, and as in your textual description) instead of one on the real line (as in your initial code). Making the matrix symmetric also makes the eigenvalues real, but it is a different matrix. – VDZ Mar 31 at 7:45