0
$\begingroup$

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}$$

enter image description here

But my code does not run.Any help?

$\endgroup$
10
  • 1
    $\begingroup$ Well, what's your question exactly? $\endgroup$ – Alone Programmer Apr 20 at 19:50
  • 2
    $\begingroup$ 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. $\endgroup$ – Alone Programmer Apr 20 at 20:26
  • 1
    $\begingroup$ 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. $\endgroup$ – Alone Programmer Apr 20 at 20:53
  • 1
    $\begingroup$ I understand that non square matrices have no eigen values but the pseudo spectra some some $\varepsilon$ and beyond will find some eigen values $\endgroup$ – Tsiantakis Apr 20 at 21:03
  • 1
    $\begingroup$ @Jan the algorithm is from the same author here $\endgroup$ – Tsiantakis Apr 24 at 13:00
0
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.