Repository: scipy/scipy
Pull request: MAINT: `linalg.logm`: avoid `ZeroDivisionError`, emit warning properly, avoid redundant squareness check
#### Reference issue
Closes gh-17136

#### What does this implement/fix?
gh-17136 noted a few issues in `linalg.logm`:

- If the norm of `A` is exactly zero and `lange` is used to evaluate it, the error estimate check raises a `ZeroDivisionError` because `lange` returns a Python `float` instead of a NumPy `float` of the appropriate dtype.
- There is a redundant shape check - it was being checked at the beginning of the function and immediately afterwards in `_logm`.

While working on this, I noticed that any warnings are `print`ed to the console rather than emitted via `warnings.warn`. 

This PR fixes these things.

#### Additional information
There are still issues here related to the an old-style API in which a `disp` parameter controls whether a warning is emitted or a tuple is returned... let's not get into that here. In the meantime, I am respecting that API: it only printed the warning when `disp=True`, so I only emit the warning when `disp=True`. Presumably emitting the warning with `warnings` rather than `print`ing is an acceptable change in behavior that can simply be documented in the release notes.

Comments in the issue suggest that if the error is high, the returned array should include NaNs rather than inaccurate values. It looks like this is because the function is not confident in its error estimate or tolerance. Since the result may still be useful to the user in some cases, I'm leaving it alone.

Problem statement:
BUG: logm sometimes raises ZeroDivisionError for matrices with zero norm
### Describe your issue.

Running `scipy.linalg.logm` on a matrix with zero norm leads to inconsistent behaviour. For example, running
```python
logm(np.zeros((2, 2), dtype=int))
```
leads to the following result:
```
array([[-46.05170186,   0.        ],
       [  0.        , -46.05170186]])
```
with an appropriate warning  displayed ("logm result may be inaccurate, approximate err = inf"). 

However, the same matrix with floats leads to the ZeroDivisionError shown below. Note that this error actually occurs after the `logm` is calculated. The actual algorithm is consistent irrespective of data type. It happens while estimating the error of the algorithm (`errtest`).

I would suggest that the behaviour for the integer matrix is better. It would be easy to get that behaviour consistently across all zero-norm matrices: namely to check that `norm(A,1)` is zero before computing `errtest`. If `norm(A,1)` is zero `errtest` can be set to infinity (and the warning message will normally be displayed).

On a closely related note: I think returning an matrix of `NaN`s is a better behaviour for these zero-norm matrices (i.e., if  `errtest` is infinite) rather than the random numbers returned now.

### Reproducing Code Example

```python
from scipy.linalg import logm
import numpy as np
logm(np.zeros((2, 2), dtype=int))
logm(np.zeros((2, 2), dtype=float))
```


### Error message

```shell
ZeroDivisionError                         Traceback (most recent call last)
Cell In [7], line 1
----> 1 logm(np.zeros((2, 2), dtype=float))

File <python-dir>/python3.10/site-packages/scipy/linalg/_matfuncs.py:204, in logm(A, disp)
    202 errtol = 1000*eps
    203 #TODO use a better error approximation
--> 204 errest = norm(expm(F)-A,1) / norm(A,1)
    205 if disp:
    206     if not isfinite(errest) or errest >= errtol:

ZeroDivisionError: float division by zero
```


### SciPy/NumPy/Python version information

1.9.1 1.23.3 sys.version_info(major=3, minor=10, micro=6, releaselevel='final', serial=0)

Hints from issue discussion:
I did some playing around with the code, and apparently python gives a runtime warning when dtype = int and a ZeroDivisionError if dtype = float, even though the values for norm(expm(F)-A,1) and norm(A,1) are the exact same. Very weird.
The problem is due to different behavior of `np.linalg.norm` and `scipy.linalg.norm`.  To reproduce the problem just consider the following (if we use `norm` from scipy package): 
```python
from scipy.linalg import norm
>>> 0/norm(np.zeros((2,2), dtype=int), 1)
nan
>>> 0/norm(np.zeros((2,2), dtype=float), 1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ZeroDivisionError: float division by zero
```
and, if we use `norm` from `numpy` directly, we get the following:
```python
>>> 0/np.linalg.norm(np.zeros((2,2), dtype=int), 1)
nan
>>> 0/np.linalg.norm(np.zeros((2,2), dtype=float), 1)
nan
```

SciPy uses `blas` or `lapack` to compute matrix [norm ](https://github.com/scidam/scipy/blob/5be8bc48fe959e8f12d99495d0323f16c7d4ea4d/scipy/linalg/_misc.py#L17) and this yields to different types of results:

```python
>>> type(norm(np.zeros((2,2), dtype=float), 1))
<class 'float'>
>>> type(norm(np.zeros((2,2), dtype=int), 1))
<class 'numpy.float64'>
```

Actually, to handle this edge case I would add an extra check for non-zero norm (somewhere in [these lines](https://github.com/scidam/scipy/blob/5be8bc48fe959e8f12d99495d0323f16c7d4ea4d/scipy/linalg/_matfuncs.py#L204) of the code):
```python
if  np.allclose(norm(A, 1), 0):
    errest = np.inf
    F = np.empty_like(A)
    F.fill(np.nan)
else:
    errest = norm(expm(F) - A, 1) / norm(A, 1)
```

Also, function [_logm](https://github.com/scidam/scipy/blob/5be8bc48fe959e8f12d99495d0323f16c7d4ea4d/scipy/linalg/_matfuncs_inv_ssq.py#L857) uses check for matrix squareness.  The same check is performed in [this ](https://github.com/scidam/scipy/blob/5be8bc48fe959e8f12d99495d0323f16c7d4ea4d/scipy/linalg/_matfuncs.py#L197) line of code. Actually, we have a double check for matrix squareness in `logm` function.

Commit message: MAINT: linalg.logm: avoid ZeroDivisionError, emit warning properly, avoid redundant squareness check

Patch to review:
diff --git a/scipy/linalg/_matfuncs.py b/scipy/linalg/_matfuncs.py
index 122581ec4053..29d24acce08b 100644
--- a/scipy/linalg/_matfuncs.py
+++ b/scipy/linalg/_matfuncs.py
@@ -1,6 +1,7 @@
 #
 # Author: Travis Oliphant, March 2002
 #
+import warnings
 from itertools import product
 
 import numpy as np
@@ -151,7 +152,7 @@ def logm(A, disp=True):
     A : (N, N) array_like
         Matrix whose logarithm to evaluate
     disp : bool, optional
-        Print warning if error in the result is estimated large
+        Emit warning if error in the result is estimated large
         instead of returning estimated error. (Default: True)
 
     Returns
@@ -194,17 +195,19 @@ def logm(A, disp=True):
            [ 1.,  4.]])
 
     """
-    A = _asarray_square(A)
+    A = np.asarray(A)  # squareness checked in `_logm`
     # Avoid circular import ... this is OK, right?
     import scipy.linalg._matfuncs_inv_ssq
     F = scipy.linalg._matfuncs_inv_ssq._logm(A)
     F = _maybe_real(A, F)
     errtol = 1000*eps
     # TODO use a better error approximation
-    errest = norm(expm(F)-A, 1) / norm(A, 1)
+    with np.errstate(divide='ignore', invalid='ignore'):
+        errest = norm(expm(F)-A, 1) / np.asarray(norm(A, 1), dtype=A.dtype).real[()]
     if disp:
         if not isfinite(errest) or errest >= errtol:
-            print("logm result may be inaccurate, approximate err =", errest)
+            message = f"logm result may be inaccurate, approximate err = {errest}"
+            warnings.warn(message, RuntimeWarning)
         return F
     else:
         return F, errest
diff --git a/scipy/linalg/tests/test_matfuncs.py b/scipy/linalg/tests/test_matfuncs.py
index c6614961a271..b9e271213176 100644
--- a/scipy/linalg/tests/test_matfuncs.py
+++ b/scipy/linalg/tests/test_matfuncs.py
@@ -18,6 +18,7 @@
                           expm, expm_frechet, expm_cond, norm, khatri_rao)
 from scipy.linalg import _matfuncs_inv_ssq
 from scipy.linalg._matfuncs import pick_pade_structure
+from scipy.linalg._matfuncs_inv_ssq import LogmExactlySingularWarning
 import scipy.linalg._expm_frechet
 
 from scipy.optimize import minimize
@@ -249,6 +250,15 @@ def test_empty(self, dt):
         assert log_a.shape == (0, 0)
         assert log_a.dtype == log_a0.dtype
 
+    @pytest.mark.parametrize('dtype', [int, float, np.float32, complex, np.complex64])
+    def test_no_ZeroDivisionError(self, dtype):
+        # gh-17136 reported inconsistent behavior in `logm` depending on input dtype:
+        # sometimes it raised an error, and sometimes it printed a warning message.
+        # check that this is resolved and that the warning is emitted properly.
+        with (pytest.warns(RuntimeWarning, match="logm result may be inaccurate"),
+              pytest.warns(LogmExactlySingularWarning)):
+            logm(np.zeros((2, 2), dtype=dtype))
+
 
 class TestSqrtM:
     def test_round_trip_random_float(self):


Code review comments:
diff_hunk: @@ -194,17 +195,19 @@ def logm(A, disp=True):
            [ 1.,  4.]])
 
     """
-    A = _asarray_square(A)
+    A = np.asarray(A)  # squareness checked in `_logm`
     # Avoid circular import ... this is OK, right?
     import scipy.linalg._matfuncs_inv_ssq
     F = scipy.linalg._matfuncs_inv_ssq._logm(A)
     F = _maybe_real(A, F)
     errtol = 1000*eps
     # TODO use a better error approximation
-    errest = norm(expm(F)-A, 1) / norm(A, 1)
+    with np.errstate(divide='ignore', invalid='ignore'):
+        errest = norm(expm(F)-A, 1) / np.asarray(norm(A, 1), dtype=A.dtype).real[()]
     if disp:
         if not isfinite(errest) or errest >= errtol:
-            print("logm result may be inaccurate, approximate err =", errest)
+            message = f"logm result may be inaccurate, approximate err = {errest}"
+            warnings.warn(message, RuntimeWarning)
path: scipy/linalg/_matfuncs.py
text: @author:
```suggestion
            warnings.warn(message, RuntimeWarning, stacklevel=2)
```