WIP - BUG: scipy.signal.lsim overflow error

[zeehio]

The scipy.signal.lsim function raises an OverflowError on the new test committed here, whereas the same test using lsim2 does not. lsim2 returns a similar output to MATLAB. Given that #5676 aims to deprecate lsim2 I believe it would be good to fix this issue. Unfortunately I do not have the skills to properly address it.

This is the error log:

ERROR: test_ltisys.TestLsim.test_07
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/travis/virtualenv/python3.5.0/lib/python3.5/site-packages/nose/case.py", line 198, in runTest
    self.test(*self.arg)
  File "/home/travis/virtualenv/python3.5.0/lib/python3.5/site-packages/scipy/signal/tests/test_ltisys.py", line 405, in test_07
    (t2, y2, x2) = lsim((z, p, k), x1, t)
  File "/home/travis/virtualenv/python3.5.0/lib/python3.5/site-packages/scipy/signal/ltisys.py", line 1325, in lsim
    expMT = linalg.expm(transpose(M))
  File "/home/travis/virtualenv/python3.5.0/lib/python3.5/site-packages/scipy/linalg/matfuncs.py", line 261, in expm
    return scipy.sparse.linalg.expm(A)
  File "/home/travis/virtualenv/python3.5.0/lib/python3.5/site-packages/scipy/sparse/linalg/matfuncs.py", line 582, in expm
    return _expm(A, use_exact_onenorm='auto')
  File "/home/travis/virtualenv/python3.5.0/lib/python3.5/site-packages/scipy/sparse/linalg/matfuncs.py", line 643, in _expm
    s = s + _ell(2**-s * h.A, 13)
  File "/home/travis/virtualenv/python3.5.0/lib/python3.5/site-packages/scipy/sparse/linalg/matfuncs.py", line 832, in _ell
    value = int(np.ceil(log2_alpha_div_u / (2 * m)))
OverflowError: cannot convert float infinity to integer

[argriffing]

fwiw this is a bug or an unnecessary limitation in the expm implementation.

[argriffing]

In case anyone wants to work on this, the ndarray passed to expm has shape (24, 24) with nonzero entries only in the first column, last row, and superdiagonal. The largest entry magnitude is like 1e90 and the smallest nonzero magnitude is like 1e-6.

[argriffing]

Here's a reduced example that gives a similar error.

>>> import numpy as np
>>> from scipy.linalg import expm
>>> expm(np.eye(2) * 1e90)
Traceback (most recent call last):
  File ".../scipy/sparse/linalg/matfuncs.py", line 642, in _expm
    s = max(int(np.ceil(np.log2(eta_5 / theta_13))), 0)
ValueError: cannot convert float NaN to integer

[zeehio]

I know this is an old issue and I still don't have a solution but anyway... I don't think that the issue lies in expm, even though it does generate the error. I would not expect for expm to work with such large numbers... Am I wrong? After all MATLAB can't compute expm(np.eye(2) * 1e90) but its lsim implementation does not have this issue...

Wouldn't it make more sense that the lsim implementation did some scaling shenanigans to guarantee the numerical stability, and then let expm do the computation?

This is the test case given in the PR, in case someone wants to try to work on this:

import numpy as np
from scipy.signal import butter, lsim
from numpy.testing import assert_approx_equal
# Test chances of overflow in lsim
(z, p, k) = butter(22, 22029, analog=True, output='zpk')
t = np.linspace(0,10E-3,num=10000)
x1 = np.cos(2*np.pi*3500*t)
(t2, y2, x2) = lsim((z, p, k), x1, t)
assert_approx_equal(np.max(y2), 0.72, significant=2)

[ilayn]

This is not reproducible on the current 1.17.3 development version. It gives the following response which seems OK. Hence closing but feel free to report back if it persists.