Pseudospectrum of non square Matrix in Python

Question

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?

Answer 1

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