Repository: scipy/scipy
Pull request: MAINT: linalg.solve: raise when diagonal matrix is exactly singular
#### Reference issue
Closes gh-22263
Closes gh-22265

#### What does this implement/fix?
gh-22263 noted that `solve` no longer raises when the input matrix is all zeros; this is because the the matrix is now treated as diagonal, and the computation does not fail. Typical NumPy warnings were emitted, but no `LinAlgWarning` was emitted because `rcond` is NaN rather than a small number.

This PR ensures that the `LinAlgWarning` is emitted when `rcond` is NaN by checking `not (rcond >= E)` instead of `rcond < lamch('E')`. It also makes the documented reason for raising a `LinAlgError` more precise, and it adds a test for the warning/error produced for an all-zero matrix. (I already added a test for an ill-conditioned but not *exactly* singular in gh-21363.)

gh-22265 noted that with this change, a test would fail due to a warning being emitted. The source of the problem was identified there (missing `lower` arg in calls to `sycon`/`hecon`, and this fixes it).

Problem statement:
BUG: linalg.solve doesn't raise an error when A is a singular matrix
### Describe your issue.

In version 1.15.0, the linalg.solve function no longer raises a `LinAlgError` when the `a` argument is a singular matrix.

### Reproducing Code Example

```python
import scipy
import scipy.linalg
import numpy as np

A = np.array([
    [0, 0],
    [0, 0]
])

I = np.array([
    [1, 0],
    [0, 1]
])

x = scipy.linalg.solve(A, I)
```


### Error message

```shell
# expected (1.14.1)

  File "/Users/ariel/.pyenv/versions/3.12.2/envs/oe312/lib/python3.12/site-packages/scipy/linalg/_basic.py", line 219, in solve
    _solve_check(n, info)
  File "/Users/ariel/.pyenv/versions/3.12.2/envs/oe312/lib/python3.12/site-packages/scipy/linalg/_basic.py", line 40, in _solve_check
    raise LinAlgError('Matrix is singular.')
numpy.linalg.LinAlgError: Matrix is singular.

# actual (1.15.0)

Warning: divide by zero encountered in divide
  x = (b1.T / diag_a).T
/Users/ariel/.pyenv/versions/3.12.2/envs/oe312/lib/python3.12/site-packages/scipy/linalg/_basic.py:295: RuntimeWarning: invalid value encountered in divide
  x = (b1.T / diag_a).T
/Users/ariel/.pyenv/versions/3.12.2/envs/oe312/lib/python3.12/site-packages/scipy/linalg/_basic.py:297: RuntimeWarning: invalid value encountered in scalar divide
  rcond = abs_diag_a.min() / abs_diag_a.max()
```


### SciPy/NumPy/Python version and system information

```shell
1.15.0 2.2.1 sys.version_info(major=3, minor=12, micro=2, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: extraframeworks
    found: true
    include directory: unknown
    lib directory: unknown
    name: Accelerate
    openblas configuration: unknown
    pc file directory: unknown
    version: unknown
  lapack:
    detection method: extraframeworks
    found: true
    include directory: unknown
    lib directory: unknown
    name: Accelerate
    openblas configuration: unknown
    pc file directory: unknown
    version: unknown
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.13.6
Compilers:
  c:
    commands: cc
    linker: ld64
    name: clang
    version: 15.0.0
  c++:
    commands: c++
    linker: ld64
    name: clang
    version: 15.0.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.11
  fortran:
    commands: gfortran
    linker: ld64
    name: gcc
    version: 13.3.0
  pythran:
    include directory: ../../../../../../private/var/folders/95/0ydz4d79163427j3k5crp3fh0000gn/T/pip-build-env-8oen8nl_/overlay/lib/python3.12/site-packages/pythran
    version: 0.17.0
Machine Information:
  build:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
  cross-compiled: false
  host:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
Python Information:
  path: /private/var/folders/95/0ydz4d79163427j3k5crp3fh0000gn/T/cibw-run-ni2eppd2/cp312-macosx_arm64/build/venv/bin/python
  version: '3.12'
```


BUG: linalg: `hecon` returns NaN incorrectly with some lower triangular data?
### Describe your issue.

In `scipy.linalg.tests.test_basic.py::TestSolve::test_hermitian`, we have the case below.

The assertions currently pass without warning, but it turns out that under the hood, `rcond` is NaN in the second `solve`. This needs further investigation.

### Reproducing Code Example

```python
import numpy as np
from numpy.testing import assert_array_almost_equal
from scipy.linalg import solve
a = np.array([[-1.84, 0.11-0.11j, -1.78-1.18j, 3.91-1.50j],
              [0, -4.63, -1.84+0.03j, 2.21+0.21j],
              [0, 0, -8.87, 1.58-0.90j],
              [0, 0, 0, -1.36]])
b = np.array([[2.98-10.18j, 28.68-39.89j],
              [-9.58+3.88j, -24.79-8.40j],
              [-0.77-16.05j, 4.23-70.02j],
              [7.79+5.48j, -35.39+18.01j]])
res = np.array([[2.+1j, -8+6j],
                [3.-2j, 7-2j],
                [-1+2j, -1+5j],
                [1.-1j, 3-4j]])
x = solve(a, b, assume_a='her')
assert_array_almost_equal(x, res)
# Also conjugate a and test for lower triangular data
x = solve(a.conj().T, b, assume_a='her', lower=True)
assert_array_almost_equal(x, res)
```


### Error message

None right now; will warn after gh-22266 merges


### SciPy/NumPy/Python version and system information

```shell
main
```

Hints from issue discussion:
I expect this is a knock-on effect of #12824, which enabled the auto-detection of matrix structure and dispatching to different solvers depending on that structure. Different solvers go through different codepaths and may detect singular matrix inputs in different ways. Detection of singularity is not a pre-computation check but rather a translation of the LAPACK error codes that are emitted during the computation. In your case, `solve()` detected that the matrix was diagonal, so it went down the `assume_a='diagonal'` codepath, which doesn't use LAPACK at all but basic pure-numpy computations. It might be worth implementing our own check in this case, but it will really _only_ trigger in this one case where you pass an all-0s matrix; a pretty contrived occurrence that may not be worth the performance hit in all the other cases.

You can use `assume_a='general'` to get the old behavior and use the general-matrix solver always.
Ok, that makes sense. Thank you for your help. Much appreciated!
https://github.com/scipy/scipy/issues/22263#issuecomment-2573509403 is right. I knew that more general banded matrices currently do not raise an error when it's singular, but we're working on it (gh-21844). ~~I thought we discussed doing this with diagonal matrices, and I assume we didn't just forget the discussion and that we implemented what we decided... but I might be wrong. @ilayn should we raise for singular diagonal matrices for consistency with the rest?~~

IIRC, this was supposed to warn - that's why we calculated `rcond` here:
https://github.com/scipy/scipy/blob/6e246d0b54dd55dc69232a0caae6772228a7ac25/scipy/linalg/_basic.py#L297
~~so I'm not sure why it *isn't* raising.~~

But because we have `0/0`, `rcond` is NaN here. In other cases, like 
```python3
A = np.array([
    [1e-300, 0],
    [0, 1]
])
```
it warns as expected. I think errors only happen when an underlying LAPACK routine raises them; otherwise we warn.

We can fix that.

Yeah, the diagonal case is LAPACK-free.

I suspect that we should also change the documented reason for `Raises` `LinAlgError` from "If the matrix is singular" to "If the computation fails because of matrix singularity". We don't (and can't) really guarantee that we always raise if the matrix is singular (even to some numerical tolerance). We get a lot of questions like this.

Commit message: DOC: linalg.solve: more precise reason for raising LinAlgError

Patch to review:
diff --git a/scipy/linalg/_basic.py b/scipy/linalg/_basic.py
index 38f4776e4977..69f4d0b02a37 100644
--- a/scipy/linalg/_basic.py
+++ b/scipy/linalg/_basic.py
@@ -45,7 +45,7 @@ def _solve_check(n, info, lamch=None, rcond=None):
     if lamch is None:
         return
     E = lamch('E')
-    if rcond < E:
+    if not (rcond >= E):  # `rcond < E` doesn't handle NaN
         warn(f'Ill-conditioned matrix (rcond={rcond:.6g}): '
              'result may not be accurate.',
              LinAlgWarning, stacklevel=3)
@@ -135,7 +135,7 @@ def solve(a, b, lower=False, overwrite_a=False,
     ValueError
         If size mismatches detected or input a is not square.
     LinAlgError
-        If the matrix is singular.
+        If the computation fails because of matrix singularity.
     LinAlgWarning
         If an ill-conditioned input a is detected.
     NotImplementedError
diff --git a/scipy/linalg/tests/test_basic.py b/scipy/linalg/tests/test_basic.py
index 782fd0865046..cbeacc8ad76f 100644
--- a/scipy/linalg/tests/test_basic.py
+++ b/scipy/linalg/tests/test_basic.py
@@ -734,6 +734,7 @@ def test_simple_complex2(self):
                                                 [-4-5.j, -3+4.j],
                                                 [6.j, 2-3.j]]))
 
+    @pytest.mark.xfail(reason="See gh-22265")
     def test_hermitian(self):
         # An upper triangular matrix will be used for hermitian matrix a
         a = np.array([[-1.84, 0.11-0.11j, -1.78-1.18j, 3.91-1.50j],
@@ -789,6 +790,25 @@ def test_ill_condition_warning(self, structure):
         with pytest.warns(LinAlgWarning, match=message):
             solve(A, b, assume_a=structure)
 
+    @pytest.mark.thread_unsafe
+    @pytest.mark.parametrize('structure',
+                             ('diagonal', 'tridiagonal', 'lower triangular',
+                              'upper triangular', 'symmetric', 'hermitian',
+                              'positive definite', 'general', None))
+    def test_exactly_singular_gh22265(self, structure):
+        n = 10
+        A = np.zeros((n, n))
+        b = np.ones(n)
+        message = "Ill-conditioned matrix..."
+        if structure in {'diagonal', None}:
+            with (pytest.warns(LinAlgWarning, match=message),
+                  np.errstate(all='ignore')):
+                solve(A, b, assume_a=structure)
+            return
+
+        with pytest.raises(LinAlgError, match="singular."):
+            solve(A, b, assume_a=structure)
+
     def test_multiple_rhs(self):
         a = np.eye(2)
         b = np.random.rand(2, 12)


Code review comments:
diff_hunk: @@ -789,6 +790,25 @@ def test_ill_condition_warning(self, structure):
         with pytest.warns(LinAlgWarning, match=message):
             solve(A, b, assume_a=structure)
 
+    @pytest.mark.thread_unsafe
+    @pytest.mark.parametrize('structure',
+                             ('diagonal', 'tridiagonal', 'lower triangular',
+                              'upper triangular', 'symmetric', 'hermitian',
+                              'positive definite', 'general', None))
+    def test_exactly_singular_gh22265(self, structure):
path: scipy/linalg/tests/test_basic.py
text: @user1:
```suggestion
    def test_exactly_singular_gh22263(self, structure):
```

diff_hunk: @@ -734,6 +734,7 @@ def test_simple_complex2(self):
                                                 [-4-5.j, -3+4.j],
                                                 [6.j, 2-3.j]]))
 
+    @pytest.mark.xfail(reason="See gh-22265")
path: scipy/linalg/tests/test_basic.py
text: @user1:
Might be worth fixing this in this PR. Up to you.

@author:
If you have a simple solution for it, perhaps. I meant to sign up to fix gh-22263, but not every pre-existing bug that I exposed while fixing it : )

@author:
Now fixed here.