# Pseudospectrum of non square Matrix in Python

I have a rectangular matrix $$A \in \mathbb{R}^{m \times n}$$

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import svdvals ,schur

A = np.array([[0, 2, 2,4],
[0, 1, 2,4],
[1, 0, 1,4],
[0, 2, 2,4],
[0, 1, 2,4],
[1, 0, 1,4],
[0, 1, 2,4],
[1, 0, 1,4]])
A

def RECTSPA(A,grid):
m,n = A.shape
I = np.eye(n)
if m >= 2*n:
a = A[n+1:m,:]
Q, R = np.linalg.qr(a)
T, U = schur(a,output='real')
R[0:n]
T = I
else:
S2, T2, Q, Z = linalg.qz(A[m-n+1,:], I[m-n+1,:])
for z in range(1,grid):
q,r = np.linalg.qr(z*T-S2)
psevals = np.min(svdvals(r))
return(psevals)



The code above implements the Pseudospectra which is 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.

The first thing to note for matrices $$A \in \mathbb{C}^{m \times n},$$ $$m \geq 2n$$ is that although we need to compute $$\sigma_{\min}(zI - A)$$ at each grid point, only the upper $$n \times n$$ portion of this matrix changes from point to point:

$$zI-A = \begin{pmatrix} z \cdot I_{n} & \\ 0 & \end{pmatrix} -\begin{pmatrix} A_{1} & \\ A_{2} & \end{pmatrix}$$

Since singular values are invariant under unitary transformations, we can replace $$A_{2}$$ by $$QA_{2}$$ for any unitary matrix $$Q.$$ In particular, we can perform a $$QR$$ factorization $$A_{2} = Q R,$$ and then $$\sigma_{\min}( zI-A) = \begin{pmatrix} z \cdot I_{n} -A_{1}& \\ -A_{2}& \end{pmatrix} -\begin{pmatrix} z \cdot I_{n} - A_{1} & \\ -R & \end{pmatrix}$$

If $$R$$ has any rows of zeros (as it certainly will if $$m > 2n)$$, these will not affect the singular values and can be removed, leaving a matrix $$S$$ of dimension $$2n \times n$$ with the following structure:

$$S= \begin{pmatrix} x & x & x & x \\ x & x & x & x \\ x & x & x & x \\ x & x & x & x \\ \hline x & x & x & x \\ & x & x & x \\ & & x & x \\ & & & x \\ \end{pmatrix}$$

But my code does not run.Any help?

• Well, what's your question exactly? Commented Apr 20, 2021 at 19:50
• Are you looking for singular values of $A$? I mean: is an error in the line that you call eigvals(A)? If yes, it is obvious that you get errors cause non-square matrices do not have eigenvalues. You are looking for scipy.linalg.svdvals to get the singular values. Commented Apr 20, 2021 at 20:26
• You are looking for eigenvalues of what? Eigenvalues of $A$? $A$ doesn't have any eigenvalues because it is not a square matrix. You should know that eigenvalues are defined only for square matrices. Commented Apr 20, 2021 at 20:53
• I understand that non square matrices have no eigen values but the pseudo spectra some some $\varepsilon$ and beyond will find some eigen values
– user38211
Commented Apr 20, 2021 at 21:03
• @Jan the algorithm is from the same author here
– user38211
Commented Apr 24, 2021 at 13:00

A couple of mistakes on sight:

def RECTSPA(A,grid):
""" description: what does this do ? """
m,n = A.shape
I = np.eye(n)
if m >= 2*n:
a = A[n+1:m,:]              # ◀◀◀ A[n:], python is 0-origin
Q, R = np.linalg.qr(a)
T, U = schur(a,output='real')
R[0:n]                      # ◀◀◀ noop
T = I
else:
S2, T2, Q, Z = linalg.qz(A[m-n+1,:], I[m-n+1,:])
for z in range(1,grid):         # ◀◀◀ 1 2 ... ??
q,r = np.linalg.qr(z*T-S2)  # ◀◀◀ S2 undefined if m >= 2n
psevals = np.min(svdvals(r))
return(psevals)


"My code does not run" questions are better asked on stackoverflow with a runnable example, which this is far from. I'd suggest start small: first Matlab with testcases where you can check the steps, then translate to Python or ... with lots of print statements. See also
numpy-for-matlab-users
Norvig, Teach Yourself Programming in Ten Years