BUG: np.linalg.eigvalsh and np.linalg.eigh returns arbitrary values for nan input, but they make different mistakes

[fdraxler]

Describe the issue:

np.linalg.eigh and np.linalg.eigvalsh produce unexpected results for certain inputs containing nan (see code example with outputs). They do so inconsistently (first example).

I would argue that any nan in a matrix should result in at least one eigenvalue set to nan. To be more precise, if the matrix is a block diagonal matrix, only the blocks containing nan should return nan.
To avoid block identification, it might be sufficient to check if a matrix contains nan and return all eigenvalues set to nan in this case. This would be an improvement over the current outputs.

Reproduce the code example:

import numpy as np

matrices = np.array([
    [
        # Expected: [nan, nan]
        [float("nan"), 0],
        [0, float("nan")]
    ],
    [
        # Expected: [nan, nan]
        [float("nan"), .1],
        [.1, 1]
    ],
    [
        # Expected: [nan, nan]
        [float("inf"), 0],
        [0, 1]
    ],
    [
        # Expected: [0, 0]
        [0, 0],
        [0, 0]
    ]
])

print(np.linalg.eigvalsh(matrices))
# [[ 0.         -0.        ]
#  [-0.14142136  0.14142136]
#  [        nan         nan]
#  [ 0.          0.        ]]


values, _ = np.linalg.eigh(matrices)
print(values)
# [[        nan         nan]
#  [-0.14142136  0.14142136]
#  [        nan         nan]
#  [ 0.          0.        ]]

Error message:

No response

NumPy/Python version information:

1.21.3 3.9.7 | packaged by conda-forge | (default, Sep 29 2021, 20:33:18)
[Clang 11.1.0 ]

[bashtage]

Misunderstood the issue. Interesting that eigvals refuses to run with LinAlgError: Array must not contain infs or NaNs. Perhaps MKL does strange things when it encounters NaNs? It used to hang forever in some circumstances when it encountered a NaN.

Edit: Identical results with OpenBLAS.

[bashtage]

LAPACK routines assume that input matrices do not contain IEEE 754 special values
such as INF or NaN values. Using these special values may cause LAPACK to return
unexpected results or become unstable.

Probably needs a nan detector for correctness. From https://spec.oneapi.io/versions/0.7/elements/oneMKL/source/domains/lapack/lapack.html

[WarrenWeckesser]

@fdraxler, thanks for reporting this inconsistency in how eigh and eigvalsh handle nan. I agree that a nan in the input should result in nan in the output, or raise an exception. (I don't think we should attempt to recognize a block structure in the input. Checking for nan will add a little overhead to every call; checking for a block structure would probably add a lot more overhead.)

As noted above, eig and eigvals will raise an exception if any value in the input is nan.

eigh, however, has several possible behaviors when the input contains nan:

1. It might raise LinAlgError: Eigenvalues did not converge.
2. It might return meaningless or nan eigenvalues and an array of eigenvectors that is all nan.
3. It might return a result that corresponds to replacing the nans with the values from the transpose of the array. (That could happen if eigh only uses the lower or upper triangular part of the array, and any nans are in the unused part.)

See the concrete examples below.

A simple fix is to follow the example of eig and eigvals and raise an error if any value in the input to eigh or eigvalsh is nan.

An alternative that is more in line with how ufuncs behave is to return all nans when there is a nan in one of the core inputs, as @fdraxler suggested above. That is, just propagate the nan with no errors or warnings. I say "core inputs" because these functions allow "stacked" inputs. The nans from different stacked inputs propagate separately; the output for any stacked input the does not contain nan is computed normally. We could do this for eig, eigvals, eigh, eigvalsh and any other functions in linalg where it seems appropriate. (EDIT: Removed an implementation detail that isn't important right now.)

I guess a third option is to just put a warning in the documentation that passing nan to these functions results in undefined behavior, so if there is a chance the user's input contains nan, they should check it before passing it to the function.

What do folks think?

---

Example of 1.

In [254]: a = np.array([[10, 0, 1.0], [0, 3, 0], [1.0, np.nan, 4.0]])

In [255]: a
Out[255]: 
array([[10.,  0.,  1.],
       [ 0.,  3.,  0.],
       [ 1., nan,  4.]])

In [256]: eigh(a)
---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-256-826195f44a67> in <module>
----> 1 eigh(a)
<snip>
LinAlgError: Eigenvalues did not converge

Example of 2.

In [257]: a = np.array([[9, 1], [np.nan, 1]])

In [258]: a
Out[258]: 
array([[ 9.,  1.],
       [nan,  1.]])

In [259]: eigh(a)
Out[259]: 
(array([nan, nan]),
 array([[nan, nan],
        [nan, nan]]))

In [260]: a = np.array([[9, 1], [1, np.nan]])

In [261]: eigh(a)
Out[261]: 
(array([-1.41421356,  1.41421356]),
 array([[nan, nan],
        [nan, nan]]))

Example of 3.

In [262]: a = np.array([[9, np.nan], [1, 1]])

In [263]: a
Out[263]: 
array([[ 9., nan],
       [ 1.,  1.]])

In [264]: eigh(a)
Out[264]: 
(array([0.87689437, 9.12310563]),
 array([[ 0.12218326, -0.99250756],
        [-0.99250756, -0.12218326]]))

In [265]: aa = np.array([[9, 1], [1, 1]])

In [266]: eigh(aa)
Out[266]: 
(array([0.87689437, 9.12310563]),
 array([[ 0.12218326, -0.99250756],
        [-0.99250756, -0.12218326]]))

[bashtage]

I would rank the options 2, 1, and then 3.

2

I think would could modify 2 so that it returns NaNs when ever the triangular section of an input array has NaN or inf values, and does no check on other triangular area. This seems to be most consistent with what is done with numeric arrays with standard values.

1

This is easy and is closest to what eig and related non-Hermetian functions do. TBH I don't know why these raise rather than simply returning NaN.

3

3 feels a bit wrong to me since the _h functions don't check or impose symmetry in the matrix that is input AFAICT.