SciPy Lyapunov 方程只把 q 转成复数，解却明显偏离

矩阵内容没有变化，只把右端项 q 的 dtype 改成 complex128，solve_discrete_lyapunov(..., method="bilinear") 的结果却不再满足原方程。SciPy 1.18.1 存在这一混合类型路径：系数矩阵保持实数，右端项为复数时，内部 Schur 分解与所选复数求解例程的输入形式可能不匹配。

同类路径还涉及连续 Lyapunov 与 Sylvester 方程。计算可能正常返回数组，因此先检查方程残差，不能仅以没有 LinAlgError 判断求解成功。

对应好方程，再比较残差

离散方程的约定是 A X Aᴴ - X + Q = 0；连续方程则是 A X + X Aᴴ = Q。两者的符号与乘法次序不同。用错残差公式，会把正确结果也判成错误。

对离散接口，可以在原求解位置保留输入，并独立计算：

residual = a @ x @ a.conj().T - x + q
residual_norm = np.linalg.norm(residual)

这里假设 a、q、x 是同一次求解对应的 NumPy 数组。大矩阵还应按矩阵和解的规模使用项目定义的相对残差与容限，不把一个固定绝对阈值用于所有数据。奇异或病态问题也可能造成大残差，不能只凭 dtype 混用就确定根因。

在求解前统一到所需复数精度

在已经确认该类型组合的情况下，可把参与运算的系数矩阵也转换成复数，使 Schur 路径使用相应形式。下面是保留原数组的离散接口调用片段：

import numpy as np
from scipy.linalg import solve_discrete_lyapunov

a_work = np.asarray(a, dtype=np.complex128)
q_work = np.asarray(q, dtype=np.complex128)
x = solve_discrete_lyapunov(a_work, q_work, method="bilinear")
residual = a_work @ x @ a_work.conj().T - x + q_work

a 和 q 需要换成原求解任务的同形方阵。complex128 是这里明确选定的工作精度；大型矩阵需要评估存储开销，特殊高精度需求也不能靠固定转成 complex128 自动满足。

不能只在求解后把 x 转成复数。错误结果的存储类型改变了，之前的分解路径不会重新计算。也不要取 q.real 来让所有输入变回实数：若虚部携带实际信息，那是在改方程。

对于 solve_sylvester(a, b, q)，若采用同类处理，两个系数矩阵和右端项都要保持一致的工作类型，并按 A X + X B = Q 检查结果。不要把离散 Lyapunov 的残差公式复制过去。

方法切换需要考虑规模与数值条件

离散接口的 method="direct" 可以作为小规模对照，但直接法会建立与矩阵维数平方相关的线性系统，不适合不看资源预算就替换大矩阵的默认路径。另一种方法算出了接近的数值，也不能取代原方程残差与业务约束的核对。

正式 1.18.1 的相关源码仍保留按实系数做 Schur 分解、随后可能选择复数 TRSYL 的组合。不要据主分支提案、问题关闭或本地补丁，就宣称某个正式版本已经消除了该路径；更不应直接编辑环境里的私有 SciPy 函数作为长期方案。

保留原依赖、方法参数和输入快照，分别覆盖全实、实系数配复数右端项、全复数及项目使用的精度。验收既比较残差，也检查解对后续计算的影响。已保存的错误解需要按任务确定重算范围；回退工作类型调整不会自动修复历史结果。