scipy.linalg.pinvh gives incorrect results

[abhinavnatarajan]

scipy.linalg.pinvh (computing the Moore-Penrose pseduoinverse) gives numerically incorrect results which do not satisfy the Moore-Penrose conditions. The results substantially differ from the output of scipy.linalg.pinv for the following Hermitian matrix.

Reproducing code example:

import numpy as np
from scipy.linalg import pinvh, pinv

x = np.array([[1,-1, 0], [-1, 2, -1], [0, -1, 1]])
y1 = pinvh(x)
y2 = pinv(x)

np.allclose(np.matmul(np.matmul(x, y1), x), x) # true
np.allclose(np.matmul(np.matmul(y1, x), y1), y1) # false
np.allclose(np.conjugate(np.matmul(x, y1)), np.matmul(x, y1)) # true
np.allclose(np.conjugate(np.matmul(y1, x)), np.matmul(y1, x)) #true

np.allclose(np.matmul(np.matmul(x, y2), x), x) # true
np.allclose(np.matmul(np.matmul(y2, x), y2), y2) # true
np.allclose(np.conjugate(np.matmul(x, y2)), np.matmul(x, y2)) # true
np.allclose(np.conjugate(np.matmul(y2, x)), np.matmul(y2, x)) #true

The analytically computed eigenvalues of this matrix are 3, 1, and 0. The reported eigenvalues from eigvals and eigvalsh are 3, 1, and -3.36770206e-17+0.j, or 3, 1, and 3.92505363e-17 respectively.

Scipy/Numpy/Python version information:

1.4.1 1.18.1 sys.version_info(major=3, minor=7, micro=6, releaselevel='final', serial=0)

[ilayn]

Ouch! I can reproduce this. Thanks @abhi01nat for this. The problem is not pinvh but the condition number doesn't cooperate. Even a very small number would return sensible results for example try pinvh(a, rcond = 1e-14). I'll take a look when I can spare some time but if you already spotted the reason and want to work on it please go ahead.

[ilayn]

Yep it's the evr driver of eigh not playing ball and loosing precision compared to svd.

sp.linalg.eigh(a, driver='ev')
Out[55]: 
(array([9.99658224e-17, 1.00000000e+00, 3.00000000e+00]),
 array([[-5.77350269e-01, -7.07106781e-01,  4.08248290e-01],
        [-5.77350269e-01,  9.71445147e-17, -8.16496581e-01],
        [-5.77350269e-01,  7.07106781e-01,  4.08248290e-01]]))

sp.linalg.eigh(a, driver='evr')
Out[56]: 
(array([2.66453526e-15, 1.00000000e+00, 3.00000000e+00]),
 array([[ 5.77350269e-01, -7.07106781e-01, -4.08248290e-01],
        [ 5.77350269e-01,  9.42055475e-16,  8.16496581e-01],
        [ 5.77350269e-01,  7.07106781e-01, -4.08248290e-01]]))

sp.linalg.eigh(a, driver='evx')
Out[57]: 
(array([9.99658224e-17, 1.00000000e+00, 3.00000000e+00]),
 array([[-5.77350269e-01, -7.07106781e-01,  4.08248290e-01],
        [-5.77350269e-01,  9.71445147e-17, -8.16496581e-01],
        [-5.77350269e-01,  7.07106781e-01,  4.08248290e-01]]))

We might want to switch to another driver for pinvh.

[keithbriggs]

You mean "losing precision".

[ev-br]

Why on earth does the eigenvalue accuracy depends on whether eigenvectors are wanted or not:

In [15]: a
Out[15]: 
array([[ 1, -1,  0],
       [-1,  2, -1],
       [ 0, -1,  1]])

In [16]: linalg.eigh(a, driver='evr')[0]
Out[16]: array([2.66453526e-15, 1.00000000e+00, 3.00000000e+00])

In [17]: linalg.eigvalsh(a, driver='evr')
Out[17]: array([3.92505363e-17, 1.00000000e+00, 3.00000000e+00])

[ev-br]

Likewise with other drivers: not asking for eigenvectors improves the zero eigenvalue. What gives?

In [19]: linalg.eigh(a, driver='ev')[0]
Out[19]: array([9.99658224e-17, 1.00000000e+00, 3.00000000e+00])

In [18]: linalg.eigvalsh(a, driver='ev')
Out[18]: array([3.92505363e-17, 1.00000000e+00, 3.00000000e+00])

In [20]: linalg.eigh(a, driver='evx')[0]
Out[20]: array([9.99658224e-17, 1.00000000e+00, 3.00000000e+00])

In [21]: linalg.eigvalsh(a, driver='evx')
Out[21]: array([3.92505363e-17, 1.00000000e+00, 3.00000000e+00])

Netlib docs have a cryptic note (https://netlib.org/lapack/explore-html-3.6.1/d2/d8a/group__double_s_yeigen_ga2ad9f4a91cddbf67fe41b621bd158f5c.html), no idea what this is meant to mean:

    If the eigenvectors are desired, the algorithm attains full
    accuracy of the computed eigenvalues only right before
    the corresponding vectors have to be computed, see steps c) and d).

[ev-br]

Opened Reference-LAPACK/lapack#997 to see what upstream thinks of this.