Repository: scipy/scipy
Pull request: BUG:linalg:Remove the 2x2 branch in expm
Closes #19584 

The 2x2 shortcut branch is not numerically stable for overflow/underflows. It might also be the numpy implementation of hyperbolic functions but there is not much performance gain to insist on that branch we can leave it for future. 

In the meantime, fixed a missing dtype branch separation. 

Problem statement:
BUG: MATLAB expm vs scipy.linalg.expm: overflow/invalid value warnings, resulting in nan
### Describe your issue.

I want to reproduce MATLAB expm-function in scipy

scipy.linalg.expm produces nan array

scipy.sparse.linalg.expm works, the result is exactly the same as in MATLAB.

All use the same algorithms, so i think there is some issue in scipy.linalg

### Reproducing Code Example

```python
import numpy as np
from scipy.linalg import expm
from scipy.sparse import csc_matrix
from scipy.sparse.linalg import expm as expm_sparse

def test():
    # this produces nan
    A = np.array([[-494.08845191, 0], [12566.3706, - 12566.3706]], dtype=np.float64)
    w = 2 * np.pi
    A = (np.pi / w) * A
    Aexp  = expm(A)
    print("scipy.linalg.expm: ",Aexp)
    # output is:
    # scipy.linalg.expm:  [[nan nan]
    #  [nan nan]]

    # this works
    A = csc_matrix([[-494.08845191, 0], [12566.3706, - 12566.3706]])
    w = 2* np.pi
    A = (np.pi / w) * A
    Aexp  = expm_sparse(A)
    print("scipy.sparse.linalg.expm: ",Aexp)
    # output is:
    # scipy.sparse.linalg.expm:    (1, 0)	5.339202199043829e-108
    #   (0, 0)	5.129273793067958e-108
    print("Aexp[0, 0] ", Aexp[0, 0])
    # Aexp[0, 0]  5.129273793067958e-108
    print("Aexp[0, 1] ", Aexp[0, 1])
    # Aexp[0, 1]  0.0
    print("Aexp[1, 0] ", Aexp[1, 0])
    # Aexp[1, 0]  5.339202199043829e-108
    print("Aexp[1, 1] ", Aexp[1, 1])
    # Aexp[1, 1]  0.0
```


### Error message

```shell
C:\Users\jne\.conda\envs\py310\lib\site-packages\scipy\linalg\_matfuncs.py:314: RuntimeWarning: overflow encountered in cosh
  coshMu = np.cosh(mu)
C:\Users\jne\.conda\envs\py310\lib\site-packages\scipy\linalg\_matfuncs.py:317: RuntimeWarning: overflow encountered in sinh
  sinchMu[mask] = np.sinh(mu[mask]) / mu[mask]
C:\Users\jne\.conda\envs\py310\lib\site-packages\scipy\linalg\_matfuncs.py:319: RuntimeWarning: invalid value encountered in multiply
  eA[..., [0], [0]] = eApD2 * (coshMu + AmD2*sinchMu)
C:\Users\jne\.conda\envs\py310\lib\site-packages\scipy\linalg\_matfuncs.py:320: RuntimeWarning: invalid value encountered in multiply
  eA[..., [0], [1]] = eApD2 * a2 * sinchMu
C:\Users\jne\.conda\envs\py310\lib\site-packages\scipy\linalg\_matfuncs.py:321: RuntimeWarning: invalid value encountered in multiply
  eA[..., [1], [0]] = eApD2 * a3 * sinchMu
C:\Users\jne\.conda\envs\py310\lib\site-packages\scipy\linalg\_matfuncs.py:322: RuntimeWarning: invalid value encountered in subtract
  eA[..., [1], [1]] = eApD2 * (coshMu - AmD2*sinchMu)

Process finished with exit code 0
```


### SciPy/NumPy/Python version and system information

```shell
C:\Users\jne\.conda\envs\py310\lib\site-packages\scipy\__config__.py:146: UserWarning: Install `pyyaml` for better output
  warnings.warn("Install `pyyaml` for better output", stacklevel=1)
{
  "Compilers": {
    "c": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "cc"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "0.29.36",
      "commands": "cython"
    },
    "c++": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "c++"
    },
    "fortran": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "gfortran"
    },
    "pythran": {
      "version": "0.14.0",
      "include directory": "C:\\Users\\runneradmin\\AppData\\Local\\Temp\\pip-build-env-qjiwhfk5\\overlay\\Lib\\site-packages/pythran"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "windows"
    },
    "build": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "windows"
    },
    "cross-compiled": false
  },
  "Build Dependencies": {
    "blas": {
      "name": "openblas",
    },
    "pybind11": {
      "name": "pybind11",
      "version": "2.11.0",
      "detection method": "config-tool",
      "include directory": "unknown"
    }
  },
  "Python Information": {
    "path": "C:\\Users\\runneradmin\\AppData\\Local\\Temp\\cibw-run-8hexkyrv\\cp310-win_amd64\\build\\venv\\Scripts\\python.exe",
    "version": "3.10"
  }
}
```

Hints from issue discussion:
I forgot to post the Matlab code. Here it is:

>> A = [-494.08845191 0; 12566.3706 -12566.3706;]

A =

   1.0e+04 *

   -0.0494         0
    1.2566   -1.2566

>> w = 2*pi

w =

    6.2832

>> expm((pi/w)*A)

ans =

  1.0e-107 *

    0.5129         0
    0.5339         0
Might be a different issue but xref gh-19020 
Probably this is because of the 2x2 special case formula saturating because I get 

```python
A = np.array([[-494.08845191, 0], [12566.3706, - 12566.3706]], dtype=np.float64)
w = 2 * np.pi
A = (np.pi / w) * A
Aexp  = la.expm(A)
C:\Users\ilhan\AppData\Local\Programs\Python\Python310\lib\site-packages\scipy\linalg\_matfuncs.py:314: RuntimeWarning: overflow encountered in cosh
  coshMu = np.cosh(mu)
C:\Users\ilhan\AppData\Local\Programs\Python\Python310\lib\site-packages\scipy\linalg\_matfuncs.py:317: RuntimeWarning: overflow encountered in sinh
  sinchMu[mask] = np.sinh(mu[mask]) / mu[mask]
C:\Users\ilhan\AppData\Local\Programs\Python\Python310\lib\site-packages\scipy\linalg\_matfuncs.py:319: RuntimeWarning: invalid value encountered in multiply
  eA[..., [0], [0]] = eApD2 * (coshMu + AmD2*sinchMu)
C:\Users\ilhan\AppData\Local\Programs\Python\Python310\lib\site-packages\scipy\linalg\_matfuncs.py:320: RuntimeWarning: invalid value encountered in multiply
  eA[..., [0], [1]] = eApD2 * a2 * sinchMu
C:\Users\ilhan\AppData\Local\Programs\Python\Python310\lib\site-packages\scipy\linalg\_matfuncs.py:321: RuntimeWarning: invalid value encountered in multiply
  eA[..., [1], [0]] = eApD2 * a3 * sinchMu
C:\Users\ilhan\AppData\Local\Programs\Python\Python310\lib\site-packages\scipy\linalg\_matfuncs.py:322: RuntimeWarning: invalid value encountered in subtract
  eA[..., [1], [1]] = eApD2 * (coshMu - AmD2*sinchMu)
```

But if I extend the array with a dummy 0

```python
AA = la.block_diag(A, [[0]])

AA
Out[6]: 
array([[ -247.04422596,     0.        ,     0.        ],
       [ 6283.1853    , -6283.1853    ,     0.        ],
       [    0.        ,     0.        ,     0.        ]])

la.expm(AA)
Out[7]: 
array([[5.12927379e-108, 0.00000000e+000, 0.00000000e+000],
       [5.33920220e-108, 0.00000000e+000, 0.00000000e+000],
       [0.00000000e+000, 0.00000000e+000, 1.00000000e+000]])
```

So much so for the fast lane I guess. We have to scale the arrays before we go fast apparently.

Commit message: BUG:linalg:Remove the 2x2 branch in expm

Patch to review:
diff --git a/scipy/linalg/_matfuncs.py b/scipy/linalg/_matfuncs.py
index 6364c72b62ff..5089c11a7f55 100644
--- a/scipy/linalg/_matfuncs.py
+++ b/scipy/linalg/_matfuncs.py
@@ -299,35 +299,14 @@ def expm(A):
         return np.exp(a)
 
     if not np.issubdtype(a.dtype, np.inexact):
-        a = a.astype(float)
+        a = a.astype(np.float64)
     elif a.dtype == np.float16:
         a = a.astype(np.float32)
 
-    # Explicit formula for 2x2 case, formula (2.2) in [1]
-    # without Kahan's method numerical instabilities can occur.
-    if a.shape[-2:] == (2, 2):
-        a1, a2, a3, a4 = (a[..., [0], [0]],
-                          a[..., [0], [1]],
-                          a[..., [1], [0]],
-                          a[..., [1], [1]])
-        mu = csqrt((a1-a4)**2 + 4*a2*a3)/2.  # csqrt slow but handles neg.vals
-
-        eApD2 = np.exp((a1+a4)/2.)
-        AmD2 = (a1 - a4)/2.
-        coshMu = np.cosh(mu)
-        sinchMu = np.ones_like(coshMu)
-        mask = mu != 0
-        sinchMu[mask] = np.sinh(mu[mask]) / mu[mask]
-        eA = np.empty((a.shape), dtype=mu.dtype)
-        eA[..., [0], [0]] = eApD2 * (coshMu + AmD2*sinchMu)
-        eA[..., [0], [1]] = eApD2 * a2 * sinchMu
-        eA[..., [1], [0]] = eApD2 * a3 * sinchMu
-        eA[..., [1], [1]] = eApD2 * (coshMu - AmD2*sinchMu)
-        if np.isrealobj(a):
-            return eA.real
-        return eA
-
-    # larger problem with unspecified stacked dimensions.
+    # Explicit formula for 2x2 case exists; formula (2.2) in [1].
+    # Howver without Kahan's method numerical instabilities can occur (See gh-19584).
+    # Hence removed here until we have a more stable implementation.
+
     n = a.shape[-1]
     eA = np.empty(a.shape, dtype=a.dtype)
     # working memory to hold intermediate arrays
diff --git a/scipy/linalg/_matfuncs_expm.pyx.in b/scipy/linalg/_matfuncs_expm.pyx.in
index eddd77fae65d..c5c577a7e343 100644
--- a/scipy/linalg/_matfuncs_expm.pyx.in
+++ b/scipy/linalg/_matfuncs_expm.pyx.in
@@ -105,18 +105,17 @@ def pick_pade_structure(cnp.ndarray[ndim=3, dtype=numpy_lapack_t] Am):
         double two_pow_s
         double temp
         double [5] theta
-        double [5] coeff        
+        double [5] coeff
         numpy_lapack_t [:, :, ::1] Amv = Am
-        cnp.ndarray[ndim=2, dtype=cnp.float64_t] absA
-        cnp.ndarray[ndim=1, dtype=cnp.float64_t] work_arr
-
-    if numpy_lapack_t in numpy_lapack_sc_t:
-        u = (2.)**(-24)
 
     dims[0] = n
     dims[1] = n
 
-    work_arr = cnp.PyArray_ZEROS(1, dims, cnp.NPY_FLOAT64, 0)
+    if numpy_lapack_t in numpy_lapack_sc_t:
+        u = (2.)**(-24)
+        work_arr = cnp.PyArray_ZEROS(1, dims, cnp.NPY_FLOAT32, 0)
+    else:
+        work_arr = cnp.PyArray_ZEROS(1, dims, cnp.NPY_FLOAT64, 0)
 
     theta[0] = 1.495585217958292e-002
     theta[1] = 2.539398330063230e-001
@@ -637,28 +636,28 @@ def pade_UV_calc(lapack_t[:, :, ::1]Am, int n, int m):
             pade_357_UV_calc_s(Am, n, m)
         elif m == 9:
             pade_9_UV_calc_s(Am, n)
-        else: 
+        else:
             pade_13_UV_calc_s(Am, n)
     elif lapack_t is double:
         if m in [3, 5, 7]:
             pade_357_UV_calc_d(Am, n, m)
         elif m == 9:
             pade_9_UV_calc_d(Am, n)
-        else: 
+        else:
             pade_13_UV_calc_d(Am, n)
     elif lapack_t is floatcomplex:
         if m in [3, 5, 7]:
             pade_357_UV_calc_c(Am, n, m)
         elif m == 9:
             pade_9_UV_calc_c(Am, n)
-        else: 
+        else:
             pade_13_UV_calc_c(Am, n)
     elif lapack_t is doublecomplex:
         if m in [3, 5, 7]:
             pade_357_UV_calc_z(Am, n, m)
         elif m == 9:
             pade_9_UV_calc_z(Am, n)
-        else: 
+        else:
             pade_13_UV_calc_z(Am, n)
     else:
         raise ValueError('Internal function "pade_UV_calc" received an unsupported dtype')


Code review comments:
diff_hunk: @@ -299,35 +299,14 @@ def expm(A):
         return np.exp(a)
 
     if not np.issubdtype(a.dtype, np.inexact):
-        a = a.astype(float)
+        a = a.astype(np.float64)
     elif a.dtype == np.float16:
         a = a.astype(np.float32)
 
-    # Explicit formula for 2x2 case, formula (2.2) in [1]
-    # without Kahan's method numerical instabilities can occur.
-    if a.shape[-2:] == (2, 2):
-        a1, a2, a3, a4 = (a[..., [0], [0]],
-                          a[..., [0], [1]],
-                          a[..., [1], [0]],
-                          a[..., [1], [1]])
-        mu = csqrt((a1-a4)**2 + 4*a2*a3)/2.  # csqrt slow but handles neg.vals
-
-        eApD2 = np.exp((a1+a4)/2.)
-        AmD2 = (a1 - a4)/2.
-        coshMu = np.cosh(mu)
-        sinchMu = np.ones_like(coshMu)
-        mask = mu != 0
-        sinchMu[mask] = np.sinh(mu[mask]) / mu[mask]
-        eA = np.empty((a.shape), dtype=mu.dtype)
-        eA[..., [0], [0]] = eApD2 * (coshMu + AmD2*sinchMu)
-        eA[..., [0], [1]] = eApD2 * a2 * sinchMu
-        eA[..., [1], [0]] = eApD2 * a3 * sinchMu
-        eA[..., [1], [1]] = eApD2 * (coshMu - AmD2*sinchMu)
-        if np.isrealobj(a):
-            return eA.real
-        return eA
-
-    # larger problem with unspecified stacked dimensions.
+    # Explicit formula for 2x2 case exists; formula (2.2) in [1].
+    # Howver without Kahan's method numerical instabilities can occur (See gh-19584).
path: scipy/linalg/_matfuncs.py
text: @user1:
```suggestion
    # However, without Kahan's method, numerical instabilities can occur (See gh-19584).
```