如何在Python中计算矩阵集合的Commutant？

符号方法（SymPy实现） #

SciPy主打数值计算，符号运算得用SymPy来做。核心思路是把变量矩阵的元素设为符号，构造交换子等于零的线性方程组，解出元素间的依赖关系，得到Commutant的通式。

步骤与示例 #

假设我们要找和集合{B1, B2, ...}中所有矩阵可交换的矩阵M：

1. 导入SymPy，定义符号矩阵M的元素
2. 对每个Bi，生成M@Bi - Bi@M = 0的等式
3. 合并所有等式，求解线性方程组
4. 用解的结果构造M的通式

代码示例（以2阶矩阵集合为例）：

import sympy as sp

# 定义符号变量矩阵M
n = 2
M = sp.Matrix([[sp.symbols(f'm{i}{j}') for j in range(n)] for i in range(n)])

# 示例矩阵集合（这里用单个矩阵演示，多矩阵只需重复构造等式）
B = sp.Matrix([[1, 2], [3, 4]])

# 构造交换子等式
commutator = M @ B - B @ M
equations = [eq for row in commutator for eq in row]

# 解线性方程组
solutions = sp.solve(equations, sp.flatten(M))

# 构造Commutant的通式
M_comm = M.subs(solutions)
print("Commutant的通式：")
sp.pprint(M_comm)

运行后会得到M的元素满足的关系，比如对于上面的B，Commutant是所有形如[[a, b], [1.5b, a+2b]]的矩阵（a,b为任意常数）。

数值方法（NumPy+SciPy实现） #

当矩阵阶数较大或只有数值矩阵时，符号方法效率太低，可通过将矩阵元素向量化，把交换子等式转化为齐次线性方程组，求解零空间得到Commutant的基。

核心原理 #

n阶矩阵M的元素可以展开成n²维向量x = [m00, m01, ..., m(n-1)(n-1)]，对于数值矩阵B，M@B - B@M = 0等价于A @ x = 0，其中A是n²×n²的系数矩阵。将集合中所有矩阵对应的A堆叠，求解整个方程组的零空间，零空间的基向量还原成矩阵就是Commutant的一组基。

代码示例 #

import numpy as np
from scipy.linalg import null_space

def commutant_matrix(B):
    """生成对应矩阵B的系数矩阵A，使得vec(M)满足A@vec(M)=0"""
    n = B.shape[0]
    I = np.eye(n)
    A = np.kron(I, B) - np.kron(B.T, I)
    return A

# 示例矩阵集合（数值矩阵）
B1 = np.array([[1, 2], [3, 4]])
B2 = np.array([[0, 1], [1, 0]])  # 可添加更多矩阵

# 堆叠所有矩阵对应的系数矩阵
A_list = [commutant_matrix(B) for B in [B1, B2]]
A = np.vstack(A_list)

# 求解零空间
null_space_basis = null_space(A)
print("Commutant的基矩阵：")
for vec in null_space_basis.T:
    # 将向量还原为矩阵
    mat = vec.reshape(B1.shape)
    print(mat)

注意事项 #

• 数值方法得到的是近似解，受浮点数精度影响，可通过设置rcond参数控制零空间的阈值
• 若零空间维度为k，则Commutant是这k个基矩阵的线性组合空间

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