BUG: linalg.solve_lyapunov/solve_sylvester: incorrect solution for real `a`/`b`, complex `q`
Author: BasilLiekensCreated Sep 13, 2026Updated Sep 17, 2026
Labelsdefectscipy.linalg
Describe your issue.
Found when working on #26103. When solve_discrete_lyapunov(a, q, method) is called with a a real array, q a complex array, method="bilinear". The returned solution is (completely) wrong. Remarkably, this issue even exists when taking q a real array of which the dtype was set to be complex and vanishes when a is also taken to be complex, see also the accompanying snippet.
Reproducing Code Example
>>> import numpy as np
>>> import scipy as sp
>>> sp.__version__
'1.18.1'
>>>
>>>
>>> rng = np.random.default_rng(seed=12345)
>>> n = 5
>>>
>>> a = rng.normal(size=(n, n))
>>> q = rng.normal(size=(n, n))
>>>
>>> x_direct = sp.linalg.solve_discrete_lyapunov(a, q, method="direct")
>>> x_bilinear = sp.linalg.solve_discrete_lyapunov(a, q, method="bilinear")
>>> np.linalg.matrix_norm(a @ x_direct @ np.conj(a.T) - x_direct + q)
np.float64(4.359459152030078e-15)
>>> np.linalg.matrix_norm(a @ x_bilinear @ np.conj(a.T) - x_bilinear + q)
np.float64(2.7917387325281098e-14)
>>>
>>>
>>> q = q.astype(np.complex128)
>>> x_direct = sp.linalg.solve_discrete_lyapunov(a, q, method="direct")
>>> x_bilinear = sp.linalg.solve_discrete_lyapunov(a, q, method="bilinear")
>>> np.linalg.matrix_norm(a @ x_direct @ np.conj(a.T) - x_direct + q)
np.float64(5.771808487210764e-15)
>>> np.linalg.matrix_norm(a @ x_bilinear @ np.conj(a.T) - x_bilinear + q)
np.float64(7.2563938234006615)
>>>
>>>
>>> a = a.astype(np.complex128)
>>> x_direct = sp.linalg.solve_discrete_lyapunov(a, q, method="direct")
>>> x_bilinear = sp.linalg.solve_discrete_lyapunov(a, q, method="bilinear")
>>> np.linalg.matrix_norm(a @ x_direct @ np.conj(a.T) - x_direct + q)
np.float64(5.771808487210764e-15)
>>> np.linalg.matrix_norm(a @ x_bilinear @ np.conj(a.T) - x_bilinear + q)
np.float64(3.610440832212805e-14)Error message
N/ASciPy/NumPy/Python version and system information
1.18.1 2.5.2 sys.version_info(major=3, minor=14, micro=0, releaselevel='final', serial=0)Source: scipy/scipy