Support singular matrices?

[TWilliamBell]

It seems like the code for non-normal matrices requires the matrix to be non-singular. Am I missing something? If not, is there any plans to support non-singular matrices in the near future?

[andrenarchy]

I'd have to look a bit deeper into this. I don't remember putting anything in that requires the matrix to be non-singular. Did you try it? Can you provide a minimal example?

[TWilliamBell]

I tried a simple example and you're right that it seems like it works for some singular matrices, but the way the NonNormalAuto is written, my matrix was passed to np.linalg.solve_triangular, which raises an error for singular matrices. Here is a minimal example:

from pseudopy import NonnormalAuto
import numpy as np

mat = np.array([[1, 0, 1],[0, 0, 0], [0, 0, 1]])
pseudo = NonnormalAuto(mat, 1e-5, 1)

numpy.linalg.LinAlgError: singular matrix: resolution failed at diagonal 0

[andrenarchy]

So it's not your input matrix that is passed to solve_triangular but a matrix derived from it with the Schur decomposition. From looking at the code, it seems like the problem arises if you have eigenvalues with multiplicity >1. I need more time to look into this so figuring this out could take a while.

In the meantime, does NonnormalMeshgrid work for you?

[TWilliamBell]

Apologies for the delay.

For NonnormalMeshgrid I'm running into problems with the plot method. I get the following messages:

/opt/anaconda3/lib/python3.7/site-packages/pseudopy/nonnormal.py:122: UserWarning: No contour levels were found within the data range.
  colors=pyplot.rcParams['axes.prop_cycle'].by_key()['color']
/opt/anaconda3/lib/python3.7/site-packages/pseudopy/utils.py:56: UserWarning: Attempting to set identical left == right == 0.020408163265306145 results in singular transformations; automatically expanding.
  pyplot.xlim(numpy.min(vertices.real), numpy.max(vertices.real))

I'm hoping to look at your code more closely and figure out if I can plot the information returned by your function myself.

If you don't mind me asking, is this something you expect you'll have the opportunity to work on? I'm reading Trefethen and Embree currently, and so in a few weeks/months I might feel that I can contribute, but obviously I'd prefer it if that was not necessary.

EDIT: Am I write to believe that if I do the following with the matrix matrix, that I'm plotting the norm of the resolvent of matrix?

pseudo = ps.NonnormalMeshgrid(matrix, real_min=-50, real_max=50, real_n=101,
                 imag_min=-50, imag_max=50, imag_n=101)
pyplot.imshow(pseudo.Vals, extent = (-50, 50, -50, 50))
pyplot.show()

[TWilliamBell]

I got NonnormalMeshgrid working, but I had to plot it myself using a method like in my edit.