BUG: signal: lsim returns incorrect results when the state space matrix A has integer type

[WarrenWeckesser]

Describe your issue.

This problem showed up in a stackoverflow question. The output of lsim did not match the output produced by another library, and in fact the result was radically different and obviously not correct.

Here is a simple example. The components of the solution should exhibit exponential decay, but instead they immediately jump to zero. And note that solution is an integer array:

In [24]: from scipy.signal import lsim

In [25]: A = np.array([[-1, 0], [0, -1]])

In [26]: B = np.array([[0], [1]])

In [27]: C = np.array([0, 0])

In [28]: lsim((A, B, C, None), U=None, T=np.linspace(0, 1, 5), X0=[1, 1])
Out[28]: 
(array([0.  , 0.25, 0.5 , 0.75, 1.  ]),
 array([0, 0, 0, 0, 0]),
 array([[1, 1],
        [0, 0],
        [0, 0],
        [0, 0],
        [0, 0]]))

If A is converted to float64, the result are as expected:

In [10]: lsim((A.astype(np.float64), B, C, None), U=None, T=np.linspace(0, 1, 5), X0=[1, 1])
Out[10]: 
(array([0.  , 0.25, 0.5 , 0.75, 1.  ]),
 array([0., 0., 0., 0., 0.]),
 array([[1.        , 1.        ],
        [0.77880078, 0.77880078],
        [0.60653066, 0.60653066],
        [0.47236655, 0.47236655],
        [0.36787944, 0.36787944]]))

The source of the problem is here:

scipy/scipy/signal/_ltisys.py

Line 2094 in 9e3f5fd

xout = np.empty((n_steps, n_states), sys.A.dtype)

The type of the xout array should not be copied directly from the type of A. When A is an integer type, out should be float64.

Reproducing Code Example

See above.

Error message

No error, just incorrect output.

SciPy/NumPy/Python version and system information

In [9]: import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info); scipy.show_config()
1.12.0.dev0+1399.9e3f5fd 1.25.1 sys.version_info(major=3, minor=11, micro=4, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /usr/include/x86_64-linux-gnu/openblas-pthread/
    lib directory: /usr/lib/x86_64-linux-gnu/openblas-pthread/
    name: openblas
    openblas configuration: USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER=1 NO_CBLAS=
      NO_LAPACK= NO_LAPACKE=1 NO_AFFINITY=1 USE_OPENMP=0 generic MAX_THREADS=64
    pc file directory: /usr/lib/x86_64-linux-gnu/pkgconfig
    version: 0.3.20
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /usr/include/x86_64-linux-gnu/openblas-pthread/
    lib directory: /usr/lib/x86_64-linux-gnu/openblas-pthread/
    name: openblas
    openblas configuration: USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER=1 NO_CBLAS=
      NO_LAPACK= NO_LAPACKE=1 NO_AFFINITY=1 USE_OPENMP=0 generic MAX_THREADS=64
    pc file directory: /usr/lib/x86_64-linux-gnu/pkgconfig
    version: 0.3.20
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.10.4
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 11.3.0
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 11.3.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 0.29.35
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 11.3.0
  pythran:
    include directory: ../../../../../py3.11.4/lib/python3.11/site-packages/pythran
    version: 0.13.1
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
Python Information:
  path: /home/warren/py3.11.4/bin/python3.11
  version: '3.11'

[roryyorke]

Related to gh-13125.

FWIW I had a go at this a while back #13235

[michaelsmatsko]

Why should the output be a float? It is perfectly reasonable for the output type to mimic the input type. If you want floats out, pass in floats. You easily imagine the alternative case where you want integer outputs.

[ilayn]

Because a few planets should align to keep things integer and still make the results correct.

It is not a both valid and common enough edge case to consider at the outset. And the consequence is just having float integers the user can truncate which is not terribly costly, unlike the checks to catch this edge case.