Repository: scipy/scipy
Pull request: BUG: linalg: fix `expm` for large arrays
#### Reference issue
Supersedes gh-20234 and closes gh-18086.

#### What does this implement/fix?
The first commit is gh-19020. The second commit is from gh-20234. The third commit fixes gh-18086 in accordance with [this comment](https://github.com/scipy/scipy/pull/20234#issuecomment-2001625323). The fourth commit adds a CI test for this bug.


Problem statement:
BUG: `scipy.linalg.expm` generates inconsistent results between runs
### Describe your issue.

For a custom `(400, 400)` matrix , `scipy.linalg.expm` generates inconsistent results between runs. I understand that the documentation mentions for `n>=400`, some approximation will be applied. However, there are just too much differences between runs, potentially causing backward incompatibility. Is there a way to make it deterministic like `scipy.sparse.linalg.expm` and matlab?

The matrix I used is uploaded here: *removed for security reasons*. 


### Reproducing Code Example

```python
import numpy as np
import scipy
data = np.load("laplacian.npy")
scipy.linalg.expm(data)
```


### Error message

```shell
N/A
```


### SciPy/NumPy/Python version and system information

```shell
1.10.0 1.23.4 sys.version_info(major=3, minor=10, micro=9, releaselevel='final', serial=0)
lapack_mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/path/to/miniconda3/envs/dlrl/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/path/to/miniconda3/envs/dlrl/include']
lapack_opt_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/path/to/miniconda3/envs/dlrl/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/path/to/miniconda3/envs/dlrl/include']
blas_mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/path/to/miniconda3/envs/dlrl/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/path/to/miniconda3/envs/dlrl/include']
blas_opt_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/path/to/miniconda3/envs/dlrl/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/path/to/miniconda3/envs/dlrl/include']
Supported SIMD extensions in this NumPy install:
    baseline = SSE,SSE2,SSE3
    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
    not found = AVX512F,AVX512CD,AVX512_KNL,AVX512_KNM,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL
```

Hints from issue discussion:
What is the inconsistency nature you mention? The entries or the eigenvalues or something else? 
The results, they are just ... different between runs.
```
In [5]: scipy.linalg.expm(data)
Out[5]:
array([[ 2.88426183e+00, -6.02503658e-01, -5.65874903e-01, ...,
        -7.52170398e-09,  2.10243969e-07,  1.34680865e-07],
       [-6.02889369e-01,  2.87222335e+00,  4.77662584e-03, ...,
         3.45144086e-07,  1.24199908e-06, -1.15446044e-06],
       [-5.88734086e-01,  4.96640387e-03,  2.85936730e+00, ...,
         1.36811044e-08,  3.25666085e-07,  2.81800070e-08],
       ...,
       [-2.31330659e-09,  1.06081452e-07,  4.04426276e-09, ...,
         2.78491643e+00, -4.07905987e-02,  9.98847016e-05],
       [ 1.01611502e-07,  5.99877620e-07,  1.51284038e-07, ...,
        -6.41006517e-02,  2.86275668e+00, -1.68356888e-01],
       [ 7.23144718e-08, -6.19470231e-07,  1.45432568e-08, ...,
         1.74381875e-04, -1.87038447e-01,  2.81744986e+00]])

In [6]: scipy.linalg.expm(data)
Out[6]:
array([[ 2.88426183e+00, -6.02503658e-01, -5.65874903e-01, ...,
        -7.52170398e-09,  2.10243969e-07,  1.34680865e-07],
       [-6.02889369e-01,  2.87222335e+00,  4.77662584e-03, ...,
         3.45144086e-07,  1.24199908e-06, -1.15446044e-06],
       [-5.88734086e-01,  4.96640387e-03,  2.85936730e+00, ...,
         1.36811044e-08,  3.25666085e-07,  2.81800070e-08],
       ...,
       [-2.31330659e-09,  1.06081452e-07,  4.04426276e-09, ...,
         2.78491643e+00, -4.07905987e-02,  9.98847016e-05],
       [ 1.01611502e-07,  5.99877620e-07,  1.51284038e-07, ...,
        -6.41006517e-02,  2.86275668e+00, -1.68356888e-01],
       [ 7.23144718e-08, -6.19470231e-07,  1.45432568e-08, ...,
         1.74381875e-04, -1.87038447e-01,  2.81744986e+00]])

In [7]: scipy.linalg.expm(data)
Out[7]:
array([[ 1.04446589e+00, -1.01868285e-02, -9.71354146e-03, ...,
        -9.04431893e-12, -1.30537946e-11,  8.54133833e-11],
       [-1.01933499e-02,  1.04446107e+00, -1.16936358e-03, ...,
        -3.18725179e-11,  9.56010610e-11,  5.98025893e-10],
       [-1.01059314e-02, -1.21582305e-03,  1.04445220e+00, ...,
        -5.66295403e-12, -5.48993744e-12,  6.31820620e-11],
       ...,
       [-2.78158814e-12, -9.79614926e-12, -1.67402233e-12, ...,
         1.04440582e+00, -1.41848165e-03, -6.35730662e-04],
       [-6.30893564e-12,  4.61747016e-11, -2.55028062e-12, ...,
        -2.22908222e-03,  1.04445711e+00, -3.90461631e-03],
       [ 4.58611825e-11,  3.20893835e-10,  3.26072648e-11, ...,
        -1.10987872e-03, -4.33788830e-03,  1.04442715e+00]])

In [8]: scipy.linalg.expm(data)
Out[8]:
array([[ 1.03355230e+00, -7.66605929e-03, -7.31059618e-03, ...,
        -3.12051959e-12, -3.46496003e-12,  3.31767941e-11],
       [-7.67096695e-03,  1.03354919e+00, -8.86366846e-04, ...,
        -8.40867011e-12,  4.09066913e-11,  2.21392513e-10],
       [-7.60591633e-03, -9.21582700e-04,  1.03354374e+00, ...,
        -1.85122411e-12, -6.10891246e-13,  2.40090584e-11],
       ...,
       [-9.59718509e-13, -2.58443929e-12, -5.47239209e-13, ...,
         1.03351499e+00, -1.07125330e-03, -4.81828271e-04],
       [-1.67462493e-12,  1.97576705e-11, -2.83781759e-13, ...,
        -1.68342797e-03,  1.03354665e+00, -2.94353145e-03],
       [ 1.78136839e-11,  1.18796683e-10,  1.23906960e-11, ...,
        -8.41191052e-04, -3.27015758e-03,  1.03352816e+00]])

In [9]: scipy.linalg.expm(data)
Out[9]:
array([[ 1.02988489e+00, -6.82196969e-03, -6.50585710e-03, ...,
        -1.78284695e-12, -1.70706759e-12,  1.86487577e-11],
       [-6.82633698e-03,  1.02988226e+00, -7.90677206e-04, ...,
        -4.15389270e-12,  2.41679821e-11,  1.22468984e-10],
       [-6.76866887e-03, -8.22091257e-04,  1.02987778e+00, ...,
        -1.03059709e-12, -9.22630450e-15,  1.33836715e-11],
       ...,
       [-5.48316127e-13, -1.27671598e-12, -3.04654165e-13, ...,
         1.02985405e+00, -9.54432013e-04, -4.29796999e-04],
       [-8.25030567e-13,  1.16729809e-11, -4.28596242e-15, ...,
        -1.49984840e-03,  1.02988014e+00, -2.62096077e-03],
       [ 1.00131156e-11,  6.57154520e-11,  6.90710158e-12, ...,
        -7.50353209e-04, -2.91179316e-03,  1.02986490e+00]])

In [10]: scipy.linalg.expm(data)
Out[10]:
array([[ 1.02804559e+00, -6.39920509e-03, -6.10278178e-03, ...,
        -1.22569479e-12, -1.05161603e-12,  1.20881721e-11],
       [-6.40330174e-03,  1.02804320e+00, -7.42575075e-04, ...,
        -2.57710961e-12,  1.63735118e-11,  7.86866116e-11],
       [-6.34931084e-03, -7.72078001e-04,  1.02803917e+00, ...,
        -6.96378267e-13,  1.46775083e-13,  8.62632498e-12],
       ...,
       [-3.76963497e-13, -7.92085221e-13, -2.05855947e-13, ...,
         1.02801777e+00, -8.95817416e-04, -4.03642865e-04],
       [-5.08248988e-13,  7.90830156e-12,  6.81824984e-14, ...,
        -1.40773811e-03,  1.02804127e+00, -2.45925877e-03],
       [ 6.49052693e-12,  4.22223333e-11,  4.45191015e-12, ...,
        -7.04692494e-04, -2.73214805e-03,  1.02802754e+00]])

In [11]: scipy.linalg.expm(data)
\Out[11]:
array([[ 1.02694021e+00, -6.14531586e-03, -5.86071056e-03, ...,
        -9.25986177e-13, -7.23439130e-13,  8.39895388e-12],
       [-6.14924997e-03,  1.02693795e+00, -7.13631095e-04, ...,
        -1.79089083e-12,  1.19207786e-11,  5.43268027e-11],
       [-6.09746087e-03, -7.41984060e-04,  1.02693418e+00, ...,
        -5.19006776e-13,  2.00543538e-13,  5.96306748e-12],
       ...,
       [-2.84787853e-13, -5.50437649e-13, -1.53423271e-13, ...,
         1.02691414e+00, -8.60583168e-04, -3.87905867e-04],
       [-3.49640167e-13,  5.75765988e-12,  9.31599502e-14, ...,
        -1.35236902e-03,  1.02693614e+00, -2.36210392e-03],
       [ 4.50966747e-12,  2.91511392e-11,  3.07744499e-12, ...,
        -6.77218345e-04, -2.62421251e-03,  1.02692329e+00]])

In [12]: scipy.linalg.expm(data)
Out[12]:
array([[ 1.02620254e+00, -5.97596032e-03, -5.69923555e-03, ...,
        -7.40305308e-13, -5.30238431e-13,  6.04766973e-12],
       [-5.97978601e-03,  1.02620037e+00, -6.94300565e-04, ...,
        -1.32944336e-12,  9.05508076e-12,  3.89043610e-11],
       [-5.92946289e-03, -7.21885518e-04,  1.02619677e+00, ...,
        -4.10120022e-13,  2.21346405e-13,  4.27039267e-12],
       ...,
       [-2.27681540e-13, -4.08609876e-13, -1.21235325e-13, ...,
         1.02617760e+00, -8.37066331e-04, -3.77395941e-04],
       [-2.56265725e-13,  4.37354615e-12,  1.02823657e-13, ...,
        -1.31541333e-03,  1.02619863e+00, -2.29727826e-03],
       [ 3.24718767e-12,  2.08756339e-11,  2.20388224e-12, ...,
        -6.58869783e-04, -2.55219353e-03,  1.02618635e+00]])
In [13]: scipy.linalg.expm(data.copy())
Out[13]:
array([[ 2.88426183e+00, -6.02503658e-01, -5.65874903e-01, ...,
        -7.52170398e-09,  2.10243969e-07,  1.34680865e-07],
       [-6.02889369e-01,  2.87222335e+00,  4.77662584e-03, ...,
         3.45144086e-07,  1.24199908e-06, -1.15446044e-06],
       [-5.88734086e-01,  4.96640387e-03,  2.85936730e+00, ...,
         1.36811044e-08,  3.25666085e-07,  2.81800070e-08],
       ...,
       [-2.31330659e-09,  1.06081452e-07,  4.04426276e-09, ...,
         2.78491643e+00, -4.07905987e-02,  9.98847016e-05],
       [ 1.01611502e-07,  5.99877620e-07,  1.51284038e-07, ...,
        -6.41006517e-02,  2.86275668e+00, -1.68356888e-01],
       [ 7.23144718e-08, -6.19470231e-07,  1.45432568e-08, ...,
         1.74381875e-04, -1.87038447e-01,  2.81744986e+00]])

In [14]: scipy.linalg.expm(data.copy())
Out[14]:
array([[ 1.04446589e+00, -1.01868285e-02, -9.71354146e-03, ...,
        -9.04431893e-12, -1.30537946e-11,  8.54133833e-11],
       [-1.01933499e-02,  1.04446107e+00, -1.16936358e-03, ...,
        -3.18725179e-11,  9.56010610e-11,  5.98025893e-10],
       [-1.01059314e-02, -1.21582305e-03,  1.04445220e+00, ...,
        -5.66295403e-12, -5.48993744e-12,  6.31820620e-11],
       ...,
       [-2.78158814e-12, -9.79614926e-12, -1.67402233e-12, ...,
         1.04440582e+00, -1.41848165e-03, -6.35730662e-04],
       [-6.30893564e-12,  4.61747016e-11, -2.55028062e-12, ...,
        -2.22908222e-03,  1.04445711e+00, -3.90461631e-03],
       [ 4.58611825e-11,  3.20893835e-10,  3.26072648e-11, ...,
        -1.10987872e-03, -4.33788830e-03,  1.04442715e+00]])
```
Ouch doesn't look fun indeed. Very strange that it oscillates between two approximate solutions apparently. I'll have a look. Do you mind pasting the array somewhere like pastebin.com or Github gists so that we don't have to download the zip file for security reasons? 
Here https://gist.github.com/liuzhenqi77/e30005e4d8accbd390510d25145cd035
@liuzhenqi77 is that array correct that it has the shape 91*400 or was there a copy paste issue ? 
I think it's gist not rendering the full array -- you could see the full array using the "Raw" button at upper right corner of the gist. Since it's a pure text file there should be no security issues.
Thanks for the tip. I've checked the array with my own Win10 box and Ubuntu VM but I can't seem to replicate the oscillation that you are showing here. I'll try to find a system that uses MKL instead of OpenBLAS when I have the chance. In the meantime could you test running the same code with a fixed seed so that the randomization does not play a role? 
For the interested I've used 

```python
import numpy as np
import scipy.linalg as la
A = np.genfromtxt('laplacian.csv', dtype=float)
la.expm(A)
```
Hi, could you maybe try running `expm` a few times by hand in the REPL like below? 
If I run it as a script it (seems) won't change.

```
Python 3.10.9 | packaged by conda-forge | (main, Feb  2 2023, 20:20:04) [GCC 11.3.0]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.10.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import numpy as np
   ...: import scipy.linalg as la
   ...: A = np.genfromtxt('laplacian.csv', dtype=float)

In [2]: la.expm(A)
Out[2]: 
array([[ 2.88426183e+00, -6.02503658e-01, -5.65874903e-01, ...,
        -7.52170398e-09,  2.10243969e-07,  1.34680865e-07],
       [-6.02889369e-01,  2.87222335e+00,  4.77662584e-03, ...,
         3.45144086e-07,  1.24199908e-06, -1.15446044e-06],
       [-5.88734086e-01,  4.96640387e-03,  2.85936730e+00, ...,
         1.36811044e-08,  3.25666085e-07,  2.81800070e-08],
       ...,
       [-2.31330659e-09,  1.06081452e-07,  4.04426276e-09, ...,
         2.78491643e+00, -4.07905987e-02,  9.98847016e-05],
       [ 1.01611502e-07,  5.99877620e-07,  1.51284038e-07, ...,
        -6.41006517e-02,  2.86275668e+00, -1.68356888e-01],
       [ 7.23144718e-08, -6.19470231e-07,  1.45432568e-08, ...,
         1.74381875e-04, -1.87038447e-01,  2.81744986e+00]])

In [3]: la.expm(A)
Out[3]: 
array([[ 2.88426183e+00, -6.02503658e-01, -5.65874903e-01, ...,
        -7.52170398e-09,  2.10243969e-07,  1.34680865e-07],
       [-6.02889369e-01,  2.87222335e+00,  4.77662584e-03, ...,
         3.45144086e-07,  1.24199908e-06, -1.15446044e-06],
       [-5.88734086e-01,  4.96640387e-03,  2.85936730e+00, ...,
         1.36811044e-08,  3.25666085e-07,  2.81800070e-08],
       ...,
       [-2.31330659e-09,  1.06081452e-07,  4.04426276e-09, ...,
         2.78491643e+00, -4.07905987e-02,  9.98847016e-05],
       [ 1.01611502e-07,  5.99877620e-07,  1.51284038e-07, ...,
        -6.41006517e-02,  2.86275668e+00, -1.68356888e-01],
       [ 7.23144718e-08, -6.19470231e-07,  1.45432568e-08, ...,
         1.74381875e-04, -1.87038447e-01,  2.81744986e+00]])

In [4]: la.expm(A)
Out[4]: 
array([[ 1.04446589e+00, -1.01868285e-02, -9.71354146e-03, ...,
        -9.04431893e-12, -1.30537946e-11,  8.54133833e-11],
       [-1.01933499e-02,  1.04446107e+00, -1.16936358e-03, ...,
        -3.18725179e-11,  9.56010610e-11,  5.98025893e-10],
       [-1.01059314e-02, -1.21582305e-03,  1.04445220e+00, ...,
        -5.66295403e-12, -5.48993744e-12,  6.31820620e-11],
       ...,
       [-2.78158814e-12, -9.79614926e-12, -1.67402233e-12, ...,
         1.04440582e+00, -1.41848165e-03, -6.35730662e-04],
       [-6.30893564e-12,  4.61747016e-11, -2.55028062e-12, ...,
        -2.22908222e-03,  1.04445711e+00, -3.90461631e-03],
       [ 4.58611825e-11,  3.20893835e-10,  3.26072648e-11, ...,
        -1.10987872e-03, -4.33788830e-03,  1.04442715e+00]])

In [5]: la.expm(A)
Out[5]: 
array([[ 1.03355230e+00, -7.66605929e-03, -7.31059618e-03, ...,
        -3.12051959e-12, -3.46496003e-12,  3.31767941e-11],
       [-7.67096695e-03,  1.03354919e+00, -8.86366846e-04, ...,
        -8.40867011e-12,  4.09066913e-11,  2.21392513e-10],
       [-7.60591633e-03, -9.21582700e-04,  1.03354374e+00, ...,
        -1.85122411e-12, -6.10891246e-13,  2.40090584e-11],
       ...,
       [-9.59718509e-13, -2.58443929e-12, -5.47239209e-13, ...,
         1.03351499e+00, -1.07125330e-03, -4.81828271e-04],
       [-1.67462493e-12,  1.97576705e-11, -2.83781759e-13, ...,
        -1.68342797e-03,  1.03354665e+00, -2.94353145e-03],
       [ 1.78136839e-11,  1.18796683e-10,  1.23906960e-11, ...,
        -8.41191052e-04, -3.27015758e-03,  1.03352816e+00]])

In [6]: import scipy; scipy.__version__
Out[6]: '1.10.0'
```

Setting the random seed still gives changing results
```
Python 3.10.9 | packaged by conda-forge | (main, Feb  2 2023, 20:20:04) [GCC 11.3.0]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.10.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import numpy as np
   ...: import scipy.linalg as la
   ...: A = np.genfromtxt('laplacian.csv', dtype=float)

In [2]: np.random.seed(0); np.random.rand(4)
Out[2]: array([0.5488135 , 0.71518937, 0.60276338, 0.54488318])

In [3]: np.random.seed(0); np.random.rand(4)
Out[3]: array([0.5488135 , 0.71518937, 0.60276338, 0.54488318])

In [4]: np.random.seed(0); la.expm(A)
Out[4]: 
array([[ 2.88426183e+00, -6.02503658e-01, -5.65874903e-01, ...,
        -7.52170398e-09,  2.10243969e-07,  1.34680865e-07],
       [-6.02889369e-01,  2.87222335e+00,  4.77662584e-03, ...,
         3.45144086e-07,  1.24199908e-06, -1.15446044e-06],
       [-5.88734086e-01,  4.96640387e-03,  2.85936730e+00, ...,
         1.36811044e-08,  3.25666085e-07,  2.81800070e-08],
       ...,
       [-2.31330659e-09,  1.06081452e-07,  4.04426276e-09, ...,
         2.78491643e+00, -4.07905987e-02,  9.98847016e-05],
       [ 1.01611502e-07,  5.99877620e-07,  1.51284038e-07, ...,
        -6.41006517e-02,  2.86275668e+00, -1.68356888e-01],
       [ 7.23144718e-08, -6.19470231e-07,  1.45432568e-08, ...,
         1.74381875e-04, -1.87038447e-01,  2.81744986e+00]])

In [5]: np.random.seed(0); la.expm(A)
Out[5]: 
array([[ 2.88426183e+00, -6.02503658e-01, -5.65874903e-01, ...,
        -7.52170398e-09,  2.10243969e-07,  1.34680865e-07],
       [-6.02889369e-01,  2.87222335e+00,  4.77662584e-03, ...,
         3.45144086e-07,  1.24199908e-06, -1.15446044e-06],
       [-5.88734086e-01,  4.96640387e-03,  2.85936730e+00, ...,
         1.36811044e-08,  3.25666085e-07,  2.81800070e-08],
       ...,
       [-2.31330659e-09,  1.06081452e-07,  4.04426276e-09, ...,
         2.78491643e+00, -4.07905987e-02,  9.98847016e-05],
       [ 1.01611502e-07,  5.99877620e-07,  1.51284038e-07, ...,
        -6.41006517e-02,  2.86275668e+00, -1.68356888e-01],
       [ 7.23144718e-08, -6.19470231e-07,  1.45432568e-08, ...,
         1.74381875e-04, -1.87038447e-01,  2.81744986e+00]])

In [6]: np.random.seed(0); la.expm(A)
Out[6]: 
array([[ 1.04446589e+00, -1.01868285e-02, -9.71354146e-03, ...,
        -9.04431893e-12, -1.30537946e-11,  8.54133833e-11],
       [-1.01933499e-02,  1.04446107e+00, -1.16936358e-03, ...,
        -3.18725179e-11,  9.56010610e-11,  5.98025893e-10],
       [-1.01059314e-02, -1.21582305e-03,  1.04445220e+00, ...,
        -5.66295403e-12, -5.48993744e-12,  6.31820620e-11],
       ...,
       [-2.78158814e-12, -9.79614926e-12, -1.67402233e-12, ...,
         1.04440582e+00, -1.41848165e-03, -6.35730662e-04],
       [-6.30893564e-12,  4.61747016e-11, -2.55028062e-12, ...,
        -2.22908222e-03,  1.04445711e+00, -3.90461631e-03],
       [ 4.58611825e-11,  3.20893835e-10,  3.26072648e-11, ...,
        -1.10987872e-03, -4.33788830e-03,  1.04442715e+00]])

In [7]: np.random.seed(0); la.expm(A)
Out[7]: 
array([[ 1.03355230e+00, -7.66605929e-03, -7.31059618e-03, ...,
        -3.12051959e-12, -3.46496003e-12,  3.31767941e-11],
       [-7.67096695e-03,  1.03354919e+00, -8.86366846e-04, ...,
        -8.40867011e-12,  4.09066913e-11,  2.21392513e-10],
       [-7.60591633e-03, -9.21582700e-04,  1.03354374e+00, ...,
        -1.85122411e-12, -6.10891246e-13,  2.40090584e-11],
       ...,
       [-9.59718509e-13, -2.58443929e-12, -5.47239209e-13, ...,
         1.03351499e+00, -1.07125330e-03, -4.81828271e-04],
       [-1.67462493e-12,  1.97576705e-11, -2.83781759e-13, ...,
        -1.68342797e-03,  1.03354665e+00, -2.94353145e-03],
       [ 1.78136839e-11,  1.18796683e-10,  1.23906960e-11, ...,
        -8.41191052e-04, -3.27015758e-03,  1.03352816e+00]])
```

Let me know if there are other things I could test or if I did something wrong. Thanks!
> Hi, could you maybe try running expm a few times by hand in the REPL like below? 

I went about 40 times however the results is all the same unfortunately. 

> Setting the random seed still gives changing results

That basically removes all the statistical stuff happening behind the scenes. Hence it makes it a code error which I don't quite understand how it might be possible.

Commit message: BUG: expm must compute A^8 when needed

Patch to review:
diff --git a/scipy/linalg/_matfuncs_expm.pyx.in b/scipy/linalg/_matfuncs_expm.pyx.in
index 67006c042fc0..eddd77fae65d 100644
--- a/scipy/linalg/_matfuncs_expm.pyx.in
+++ b/scipy/linalg/_matfuncs_expm.pyx.in
@@ -1,5 +1,10 @@
 # cython: language_level=3
+# cython: boundscheck=False
+# cython: initializedcheck=False
+# cython: wraparound=False
+# cython: cdivision=True
 # cython: cpow=True
+
 cimport cython
 import math
 import random
@@ -18,6 +23,16 @@ from scipy.linalg.cython_blas cimport (sgemm, saxpy, sscal, scopy, sgemv,
 from ._cythonized_array_utils cimport (lapack_t, lapack_cz_t,
         lapack_sd_t, swap_c_and_f_layout)
 
+ctypedef fused numpy_lapack_t:
+    cnp.float32_t
+    cnp.float64_t
+    cnp.complex64_t
+    cnp.complex128_t
+
+ctypedef fused numpy_lapack_sc_t:
+    cnp.float32_t
+    cnp.complex64_t
+
 __all__ = ['pick_pade_structure', 'pade_UV_calc']
 
 # See GH-14813
@@ -73,50 +88,35 @@ cdef double norm1(lapack_t[:, ::1] A):
         return temp
     finally:
         free(work)
-# ============================================================================
 
-# ========================= kth power norm : d ===============================
-@cython.wraparound(False)
-@cython.boundscheck(False)
-@cython.initializedcheck(False)
-cdef double kth_power_norm_d(double* A, double* v1, double* v2, Py_ssize_t n, int spins) noexcept:
-    cdef int k, int_one = 1, intn = <int>n
-    cdef double one =1., zero = 0.
-    for k in range(spins):
-        dgemv(<char*>'C', &intn, &intn, &one, &A[0], &intn, &v1[0], &int_one, &zero, &v2[0], &int_one)
-        dcopy(&intn, &v2[0], &int_one, &v1[0], &int_one)
-    k = idamax(&intn, &v1[0], &int_one)
-    return v1[k-1]
-# ============================================================================
 
 # ====================== pick_pade_structure : s, d, c, z ====================
-{{for tname, pfx, dchar, dt in zip(typenames, prefixes, ['f', 'd', 'F', 'D'],
-                                   ['float32', 'float64', 'complex64', 'complex128'],
-                                  )}}
-@cython.cdivision(True)
-@cython.wraparound(False)
-@cython.boundscheck(False)
-@cython.initializedcheck(False)
-cdef pick_pade_structure_{{ pfx }}({{ tname }}[:, :, ::1] Am):
-    cdef Py_ssize_t n = Am.shape[1], i, j, k
-    cdef int lm = 0, s = 0, intn = <int>n, n2= intn*intn, int_one = 1
-    cdef {{ tname }}[:, :, ::1] Amv = Am
+def pick_pade_structure(cnp.ndarray[ndim=3, dtype=numpy_lapack_t] Am):
     cdef:
-        {{ tname }} one = 1.0, zero = 0.0
+        cnp.npy_intp n = Am.shape[1]
+        cnp.npy_intp i
+        cnp.npy_intp [2] dims
+        int lm = 0
+        int s = 0
         double normA
         double d4, d6, d8, d10
         double eta0, eta1, eta2, eta3, eta4
-        double u = (2.)**({{if pfx in ['s', 'c']}}-24.{{else}}-53{{endif}})
+        double u = (2.)**(-53)
         double two_pow_s
         double temp
-    cdef double *absA = <double*>malloc(n*n*sizeof(double))
-    cdef double *work = <double*>malloc(n*sizeof(double))
-    cdef double *work2 = <double*>malloc(n*sizeof(double))
-    if not work or not work2 or not absA:
-        raise MemoryError('Internal function "pick_pade_structure" failed to allocate memory.')
-    cdef double [5] theta
-    cdef double [5] coeff
-    cdef char* cN = 'N'
+        double [5] theta
+        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)
 
     theta[0] = 1.495585217958292e-002
     theta[1] = 2.539398330063230e-001
@@ -128,85 +128,103 @@ cdef pick_pade_structure_{{ pfx }}({{ tname }}[:, :, ::1] Am):
     coeff[2] = u*4487938430976000.
     coeff[3] = u*5914384781877411840000.
     coeff[4] = u*113250775606021113483283660800000000.
-    try:
-        for j in range(n):
-            work[j] = 1.
 
-        # {{ pfx }}copy(&n2, &A[0, 0], &int_one, &Amv[0, 0, 0], &int_one)
-        for i in range(n):
-            for j in range(n):
-                absA[i*n + j] = abs(Am[0, i, j])
-
-        # First spin = normest(|A|, 1), increase m when spun more
-        normA = kth_power_norm_d(absA, work, work2, n, 1)
-        {{ pfx }}gemm(cN, cN, &intn, &intn, &intn, &one, &Amv[0, 0, 0], &intn, &Amv[0, 0, 0], &intn, &zero, &Amv[1, 0, 0], &intn)
-        {{ pfx }}gemm(cN, cN, &intn, &intn, &intn, &one, &Amv[1, 0, 0], &intn, &Amv[1, 0, 0], &intn, &zero, &Amv[2, 0, 0], &intn)
-        {{ pfx }}gemm(cN, cN, &intn, &intn, &intn, &one, &Amv[2, 0, 0], &intn, &Amv[1, 0, 0], &intn, &zero, &Amv[3, 0, 0], &intn)
-        d4 = norm1(Amv[2]) ** (1./4.)
-        d6 = norm1(Amv[3]) ** (1./6.)
-        eta0 = max(d4, d6)
-        eta1 = eta0
-
-        # m = 3
-        temp = kth_power_norm_d(absA, work, work2, n, 6)
-        lm = max(<int>ceil(log2(temp/normA/coeff[0])/6), 0)
-        if eta0 < theta[0] and lm == 0:
-            return 3, s
-
-        # m = 5
-        temp = kth_power_norm_d(absA, work, work2, n, 4)
-        lm = max(<int>ceil(log2(temp/normA/coeff[1])/10), 0)
-        if eta1 < theta[1] and lm == 0:
-            return 5, s
-
-        # m = 7
-        if n < 400:
-            {{ pfx }}gemm(cN, cN, &intn, &intn, &intn, &one, &Amv[2, 0, 0], &intn, &Amv[2, 0, 0], &intn, &zero, &Amv[4, 0, 0], &intn)
-            d8 = norm1(Amv[4, :, :]) ** (1./8.)
-        else:
-            d8  = _norm1est(np.asarray(Am[0]), m=8) ** (1./8.)
-
-        eta2 = max(d6, d8)
-        temp = kth_power_norm_d(absA, work, work2, n, 4)
-        lm = max(<int>ceil(log2(temp/normA/coeff[2])/14), 0)
-        if eta2 < theta[2] and lm == 0:
-            return 7, s
-
-        # m = 9
-        temp = kth_power_norm_d(absA, work, work2, n, 4)
-        lm = max(<int>ceil(log2(temp/normA/coeff[3])/18), 0)
-        if eta2 < theta[3] and lm == 0:
-            return 9, s
-
-        # m = 13
-        # Scale-square
-        if n < 400:
-            {{ pfx }}gemm(cN, cN, &intn, &intn, &intn, &one, &Amv[3, 0, 0], &intn, &Amv[2, 0, 0], &intn, &zero, &Amv[4, 0, 0], &intn)
-            d10 = norm1(Amv[4, :, :]) ** (1./10.)
-        else:
-            d10  = _norm1est(np.asarray(Am[0]), m=10) ** (1./10.)
-
-        eta3 = max(d8, d10)
-        eta4 = min(eta2, eta3)
-        s = max(<int>ceil(log2(eta4/theta[4])), 0)
-        if s != 0:
-            two_pow_s = 2.** (-s)
-            dscal(&n2, &two_pow_s, absA, &int_one)
-            # kth_power_norm has spun 19 times already
-            two_pow_s = 2.** ((-s)*19.)
-            for i in range(n):
-                work[i] *= two_pow_s
-            normA *= 2.**(-s)
-        temp = kth_power_norm_d(absA, work, work2, n, 8)
-        s += max(<int>ceil(log2(temp/normA/coeff[4])/26), 0)
-        return 13, s
-    finally:
-        free(work)
-        free(work2)
-        free(absA)
 
+    absA = np.abs(Am[0, :, :])
+
+    # First spin = normest(|A|, m=1), increase m when spun more
+    np.matmul(absA, work_arr, out=work_arr)
+    normA = np.max(work_arr)
+
+    np.matmul(Am[0, :, :], Am[0, :, :], out=Am[1, :, :])
+    np.matmul(Am[1, :, :], Am[1, :, :], out=Am[2, :, :])
+    np.matmul(Am[2, :, :], Am[1, :, :], out=Am[3, :, :])
+    d4 = norm1(Amv[2]) ** (1./4.)
+    d6 = norm1(Amv[3]) ** (1./6.)
+    eta0 = max(d4, d6)
+    eta1 = eta0
+
+    # m = 3
+    # -------
+    # 1-norm of A ** 7
+    for i in range(6):
+        np.matmul(absA, work_arr, out=work_arr)
+    temp = np.max(work_arr)
+
+    lm = max(<int>ceil(log2(temp/normA/coeff[0])/6), 0)
+    if eta0 < theta[0] and lm == 0:
+        return 3, s
+
+    # m = 5
+    # -------
+    # 1-norm of A ** 11
+    for i in range(4):
+        np.matmul(absA, work_arr, out=work_arr)
+    temp = np.max(work_arr)
+
+    lm = max(<int>ceil(log2(temp/normA/coeff[1])/10), 0)
+    if eta1 < theta[1] and lm == 0:
+        return 5, s
+
+    # m = 7
+    # -------
+    if n < 400:
+        np.matmul(Am[2, :, :], Am[2, :, :], out=Am[4, :, :])
+        d8 = norm1(Amv[4, :, :]) ** (1./8.)
+    else:
+        d8  = _norm1est(Am[0], m=8) ** (1./8.)
+
+    eta2 = max(d6, d8)
+
+    # 1-norm of A ** 15
+    for i in range(4):
+        np.matmul(absA, work_arr, out=work_arr)
+    temp = np.max(work_arr)
+
+    lm = max(<int>ceil(log2(temp/normA/coeff[2])/14), 0)
+    if eta2 < theta[2] and lm == 0:
+        return 7, s
+
+    # m = 9
+    # -------
+    # 1-norm of A ** 19
+    for i in range(4):
+        np.matmul(absA, work_arr, out=work_arr)
+    temp = np.max(work_arr)
+
+    lm = max(<int>ceil(log2(temp/normA/coeff[3])/18), 0)
+    if eta2 < theta[3] and lm == 0:
+        return 9, s
+
+    # m = 13
+    # -------
+    # Scale-square
+    if n < 400:
+        np.matmul(Am[3, :, :], Am[2, :, :], out=Am[4, :, :])
+        d10 = norm1(Amv[4, :, :]) ** (1./10.)
+    else:
+        d10  = _norm1est(Am[0], m=10) ** (1./10.)
+
+    eta3 = max(d8, d10)
+    eta4 = min(eta2, eta3)
+    s = max(<int>ceil(log2(eta4/theta[4])), 0)
+    if s != 0:
+        two_pow_s = 2.** (-s)
+        absA *= two_pow_s
+
+        # work_arr has spun 19 times already
+        two_pow_s = 2.** ((-s)*19.)
+        work_arr *= two_pow_s
+        normA *= 2.**(-s)
+
+    # 1-norm of A ** 27
+    for i in range(8):
+        np.matmul(absA, work_arr, out=work_arr)
+    temp = np.max(work_arr)
+
+    s += max(<int>ceil(log2(temp/normA/coeff[4])/26), 0)
+    return 13, s
 
-{{endfor}}
 # ============================================================================
 
 # ====================== pade_m_UV_calc : s, d, c, z =========================
@@ -324,6 +342,7 @@ cdef void pade_9_UV_calc_{{ pfx }}({{ tname }}[:, :, ::]Am, int n) noexcept nogi
     cdef int i, info, n2 = n*n, int_one = 1
     cdef {{ tname }} one = 1.0, zero = 0.0, neg_one = -1.0
     cdef {{if pfx in ['s', 'c']}}float{{else}}double{{endif}} two = 2.0
+    cdef char* cN = 'N'
 
     if not (work and ipiv):
         raise MemoryError('Internal function "pade_9_UV_calc" failed to allocate memory.')
@@ -339,6 +358,11 @@ cdef void pade_9_UV_calc_{{ pfx }}({{ tname }}[:, :, ::]Am, int n) noexcept nogi
         b[7] = 3960.
         b[8] = 90.
 
+        # Compute this in case we skipped it before
+        if n >= 400:
+            {{ pfx }}gemm(cN, cN, &n, &n, &n, &one, &Am[2, 0, 0], &n,
+                &Am[2, 0, 0], &n, &zero, &Am[4, 0, 0], &n)
+
         # Utilize the unused powers of Am as scratch memory
         # U = Am[0] @ (b[9]*Am[4] + b[7]*Am[3] + b[5]*Am[2] + b[3]*Am[1] + b[1]*I_n)
         {{ pfx }}copy(&n2, &Am[1, 0, 0], &int_one, &work[0], &int_one)
@@ -638,18 +662,3 @@ def pade_UV_calc(lapack_t[:, :, ::1]Am, int n, int m):
             pade_13_UV_calc_z(Am, n)
     else:
         raise ValueError('Internal function "pade_UV_calc" received an unsupported dtype')
-
-
-@cython.initializedcheck(False)
-def pick_pade_structure(lapack_t[:, :, ::1]A):
-    """Helper functions for expm to choose Pade approximation order"""
-    if lapack_t is float:
-        return pick_pade_structure_s(A)
-    elif lapack_t is double:
-        return pick_pade_structure_d(A)
-    elif lapack_t is floatcomplex:
-        return pick_pade_structure_c(A)
-    elif lapack_t is doublecomplex:
-        return pick_pade_structure_z(A)
-    else:
-        raise ValueError('Internal function "pick_pade_structure" received an unsupported dtype.')


Code review comments:
diff_hunk: @@ -728,6 +729,24 @@ def test_readonly(self):
         a.flags.writeable = False
         expm(a)
 
+    def test_gh18086(self):
+        A = np.zeros((400, 400), dtype=float)
+        i = np.random.default_rng(0).integers(0, 399, 500)
+        j = np.random.default_rng(100).integers(0, 399, 500)
+        A[i, j] = np.random.default_rng(1000).random(500)
path: scipy/linalg/tests/test_matfuncs.py
text: @user1:
This pattern of multiple default_rngs with varying seeds is a bit unusual. Is there a reason to not use a single `rng=np.random.default_rng(some seed) ` and drawing from it throughout the test?