sympy.linsolve求解耦合线性代数方程组返回错误结果的问题

问题背景 #

本次涉及的耦合线性代数方程组为欠定系统，存在无穷多组解。使用Sympy的linsolve方法求解该系统时，返回的解空间不符合预期。

复现代码 #

import numpy as np
import sympy

A = np.array([[1,1,0,1,1,-3,-1,1,],[0,0,0,-1,-1,0,-1,-2,],[0,0,0,0,0,1,1,1,]])
b = np.array([[0],[0],[0]])

variable_symbols = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h']

# 转换为Sympy矩阵格式
A_sympy = sympy.Matrix(A)
b_sympy = sympy.Matrix(b)

# 生成符号变量
syms = sympy.symbols(variable_symbols)
syms = list(syms)

# 调用线性求解器
sol_sympy = sympy.linsolve((A_sympy,b_sympy),syms)
print(sol_sympy)

输出对比 #

异常输出（Sympy 1.6.1版本） #

FiniteSet((-b - g - 2*h, b, c, -c - g - 2*h, e, -g - h, g, h))

正确预期输出 #

FiniteSet((-b - g - 2*h, b, c, -e - g - 2*h, e, -g - h, g, h))

错误验证 #

两种输出并非是自由参数选择不同导致的等价表达，可通过代入验证真伪：
取正确解的自由参数b=0、c=1、e=0、g=0、h=0，可得到合法解向量(a=0, b=0, c=1, d=0, e=0, f=0, g=0, h=0)，该解无法通过Sympy 1.6.1返回的错误表达式构造得到。

修复方案 #

该问题为Sympy 1.6.1版本的已知缺陷，将Sympy升级至1.7.1及以上版本即可解决。

内容的提问来源于stack exchange，提问作者chr218