如何将SymPy多元多项式转换为矩阵表示并导出为NumPy数组？

假设我们有一个多元多项式，通过SymPy定义如下：

import sympy

x_1, x_2, x_3, x_4 = sympy.symbols('x_1 x_2 x_3 x_4')

expressions = [
           -x_1*x_1-x_2*x_2+x_1*x_2,
           x_2*x_2-x_1*x_2,
           x_1*x_1-x_1*x_2,
           x_1*x_2,
           x_1*x_1-x_1*x_3,
           x_1*x_3,
           x_3*x_3-x_2*x_3,
           -x_2*x_2-x_4*x_4+x_2*x_4,
           x_2*x_2-x_2*x_3,
           x_2*x_3,
           -x_3*x_3-x_4*x_4+x_3*x_4,
           x_3*x_4,
          ]

model = sympy.Poly(sympy.Add(*expressions))
# 输出结果：Poly(x_1**2 - x_2*x_3 + x_2*x_4 + 2*x_3*x_4 - 2*x_4**2, x_1, x_2, x_3, x_4, domain='ZZ')

对于变量[x_1, x_2, x_3, x_4]，我们可以将多项式系数表示为4×4方阵：

• 对角线元素对应平方项x_i*x_i的系数
• 非对角线元素对应x_i*x_j的系数

示例矩阵如下：

[[1, 0,  0,  0],
 [0, 0, -1,  1],
 [0, 0,  0,  2],
 [0, 0, 0, -2]
]

需求是：针对任意含变量[x_1, x_2, ..., x_N]的SymPy多项式，提取系数并构建上述格式的SymPy矩阵，最终转换为NumPy数组用于后续计算。

实现步骤 #

我们可以编写一个通用函数完成该任务，核心步骤如下：

1. 获取多项式的变量列表，确定矩阵维度N
2. 初始化N×N的零矩阵
3. 遍历多项式的每一项，根据项的类型（平方项/交叉项）填充矩阵对应位置的系数
4. 将SymPy矩阵转换为NumPy数组

代码实现 #

import sympy
import numpy as np

def poly_to_coeff_matrix(poly):
    # 获取多项式的变量列表
    vars_list = poly.gens
    N = len(vars_list)
    # 初始化N×N的零矩阵
    coeff_matrix = sympy.zeros(N, N)
    
    # 遍历多项式的每一项及其系数
    for term, coeff in poly.terms():
        exponents = list(term)
        # 筛选出指数不为0的变量索引
        non_zero_indices = [i for i, exp in enumerate(exponents) if exp > 0]
        
        if len(non_zero_indices) == 1:
            # 处理平方项x_i²
            idx = non_zero_indices[0]
            coeff_matrix[idx, idx] = coeff
        elif len(non_zero_indices) == 2:
            # 处理交叉项x_i*x_j
            i, j = non_zero_indices
            coeff_matrix[i, j] = coeff
    
    return coeff_matrix

# 测试示例
if __name__ == "__main__":
    x_1, x_2, x_3, x_4 = sympy.symbols('x_1 x_2 x_3 x_4')
    expressions = [
               -x_1*x_1-x_2*x_2+x_1*x_2,
               x_2*x_2-x_1*x_2,
               x_1*x_1-x_1*x_2,
               x_1*x_2,
               x_1*x_1-x_1*x_3,
               x_1*x_3,
               x_3*x_3-x_2*x_3,
               -x_2*x_2-x_4*x_4+x_2*x_4,
               x_2*x_2-x_2*x_3,
               x_2*x_3,
               -x_3*x_3-x_4*x_4+x_3*x_4,
               x_3*x_4,
              ]
    model = sympy.Poly(sympy.Add(*expressions))
    
    # 构建SymPy系数矩阵
    sympy_matrix = poly_to_coeff_matrix(model)
    print("SymPy系数矩阵：")
    print(sympy_matrix)
    
    # 转换为NumPy数组
    numpy_array = np.array(sympy_matrix).astype(np.int64)
    print("\nNumPy数组：")
    print(numpy_array)

输出结果 #

运行代码后会得到如下输出：

SymPy系数矩阵：
Matrix([
[1, 0,  0,  0],
[0, 0, -1,  1],
[0, 0,  0,  2],
[0, 0,  0, -2]])

NumPy数组：
[[ 1  0  0  0]
 [ 0  0 -1  1]
 [ 0  0  0  2]
 [ 0  0  0 -2]]

补充说明 #

• 该函数支持任意数量变量的多元多项式，输入需为SymPy的Poly对象
• 若需要对称矩阵（即交叉项x_i*x_j的系数同时填充到(i,j)和(j,i)位置），可在代码中添加coeff_matrix[j, i] = coeff
• 转换为NumPy数组时，可根据系数类型选择对应的数据类型（如float64用于含小数的系数）

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