Intrinsic calculation error exists in scipy.linalg.expm

Hi,

I am trying to use scipy.linalg.expm to obtain the exponential of a matrix. It was okay in the previous version (like 1.13.1 or 1.14.1). However, I upgrade the scipy to 1.15.1 recently and the result gose not correctly in some parameters region. I notice there is some rebuild in expm so I am douting that’s the reason for the inaccurate result.

Please see the Fig.

Thanks for looking.

Hi @zeling Indeed it looks like a bug. Can you please include the data as text so we can copy/paste and try ourselves?

I am sorry I forgot this. Here it is.

from scipy.linalg import expm
import numpy as np

A1=np.array([[ 0.        +0.j        ,  0.08660254-0.05j      ,
         0.        +0.j        ,  0.        +0.j        ,
         0.        +0.j        ],
       [-0.08660254-0.05j      ,  0.        +0.j        ,
         0.12247449-0.07071068j,  0.        +0.j        ,
         0.        +0.j        ],
       [ 0.        +0.j        , -0.12247449-0.07071068j,
         0.        +0.j        ,  0.15      -0.08660254j,
         0.        +0.j        ],
       [ 0.        +0.j        ,  0.        +0.j        ,
        -0.15      -0.08660254j,  0.        +0.j        ,
         0.17320508-0.1j       ],
       [ 0.        +0.j        ,  0.        +0.j        ,
         0.        +0.j        , -0.17320508-0.1j       ,
         0.        +0.j        ]])

A2=np.array([[ 0.   +0.j   ,  0.087-0.05j ,  0.   +0.j   ,  0.   +0.j   ,
         0.   +0.j   ],
       [-0.087-0.05j ,  0.   +0.j   ,  0.122-0.071j,  0.   +0.j   ,
         0.   +0.j   ],
       [ 0.   +0.j   , -0.122-0.071j,  0.   +0.j   ,  0.15 -0.087j,
         0.   +0.j   ],
       [ 0.   +0.j   ,  0.   +0.j   , -0.15 -0.087j,  0.   +0.j   ,
         0.173-0.1j  ],
       [ 0.   +0.j   ,  0.   +0.j   ,  0.   +0.j   , -0.173-0.1j  ,
         0.   +0.j   ]])
expm(A2)

I can’t replicate this in my system. This is what I am getting with 1.15.1


C:\Users\ilhan> ipython
Python 3.12.1 (tags/v3.12.1:2305ca5, Dec  7 2023, 22:03:25) [MSC v.1937 64 bit (AMD64)]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.25.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import numpy as np

In [2]: from scipy.linalg import expm

In [3]: import scipy as sp

In [4]: sp.__version__
Out[4]: '1.15.1'

In [5]: A1=np.array([[ 0.        +0.j        ,  0.08660254-0.05j,         0.        +0.j        ,  0.        +0.j        ,         0.        +0.j        ],
   ...:              [-0.08660254-0.05j      ,  0.        +0.j  ,         0.12247449-0.07071068j,  0.        +0.j        ,         0.        +0.j        ],
   ...:        [ 0.        +0.j        , -0.12247449-0.07071068j,         0.        +0.j        ,  0.15      -0.08660254j,         0.        +0.j        ],
   ...:        [ 0.        +0.j        ,  0.        +0.j        ,        -0.15      -0.08660254j,  0.        +0.j        ,         0.17320508-0.1j       ],
   ...:        [ 0.        +0.j        ,  0.        +0.j        ,         0.        +0.j        , -0.17320508-0.1j       ,         0.        +0.j        ]])
   ...:
   ...: A2=np.array([[ 0.   +0.j   ,  0.087-0.05j ,  0.   +0.j   ,  0.   +0.j   ,         0.   +0.j   ],
   ...:              [-0.087-0.05j ,  0.   +0.j   ,  0.122-0.071j,  0.   +0.j   ,         0.   +0.j   ],
   ...:              [ 0.   +0.j   , -0.122-0.071j,  0.   +0.j   ,  0.15 -0.087j,         0.   +0.j   ],
   ...:              [ 0.   +0.j   ,  0.   +0.j   , -0.15 -0.087j,  0.   +0.j   ,         0.173-0.1j  ],
   ...:              [ 0.   +0.j   ,  0.   +0.j   ,  0.   +0.j   , -0.173-0.1j  ,         0.   +0.j   ]])

In [6]: expm(A1)
Out[6]:
array([[ 9.95012479e-01-6.76162470e-20j,  8.61706080e-02-4.97506240e-02j,
         3.51790042e-03-6.09318231e-03j, -5.71230046e-13-4.06211183e-04j,
        -1.01722384e-05-1.76188336e-05j],
       [-8.61706080e-02-4.97506240e-02j,  9.85062354e-01-1.17664121e-18j,
         1.21254329e-01-7.00062198e-02j,  6.07283784e-03-1.05184637e-02j,
        -1.14246002e-12-8.12422367e-04j],
       [ 3.51790042e-03+6.09318231e-03j, -1.21254329e-01-7.00062198e-02j,
         9.75162063e-01-3.39199667e-18j,  1.47759357e-01-8.53089040e-02j,
         8.59548225e-03-1.48878120e-02j],
       [ 5.71230028e-13-4.06211183e-04j,  6.07283784e-03+1.05184637e-02j,
        -1.47759357e-01-8.53089040e-02j,  9.65228550e-01+1.25717314e-18j,
         1.71192277e-01-9.88379079e-02j],
       [-1.01722384e-05+1.76188336e-05j,  1.14245987e-12-8.12422367e-04j,
         8.59548225e-03+1.48878120e-02j, -1.71192277e-01-9.88379079e-02j,
         9.80116362e-01+4.30141027e-18j]])

I wasn’t able to reproduce the problem on Colab, either.

Same here: cannot repro on linux (bog-standard ubuntu, conda env from environment.yml.

Hi I have reproduced it in my colleague’s pc. We both use Anaconda in Windows system to run it. At first update to 1.15.1 by “conda update scipy”, the result is incorrect. Then we get back to 1.14.1 by “conda install scipy=1.14.1” and the result becomes correct.

I don’t see why it should be wrong unless the input array is different or MKL/OpenBLAS difference.

Can you logm the result and see what it gives back compared to the array you have in mind?

We both use Anaconda in Windows system

Do you get scipy from conda-forge or from the default channel?


It should be defaults channel.

@zeling Sorry for asking this again but screenshots are not valid ways of communication. For example we need to see what logm(A3) gives to find a hint of what the original array is and so on.

Please paste results like we have done above using a simple ipython console so that we can see what the data and the results are. Or run everything in a single cell and paste the code and the results separately. I don’t quite understand what I am seeing other than what you already pasted above.

Also I think it is better to open an issue on the GitHub repo instead of here.

So that we don’t ping everyone here.

1 Like

Discussion moved to BUG: linalg.expm: bug on Windows / conda · Issue #22558 · scipy/scipy · GitHub.