Sympy的diff函数为何返回4D数组？预期为(3,1)的2D数组

复现代码 #

import sympy as sym
from sympy import symbols

x1,x2,x3 = symbols('x1,x2,x3')
X = sym.Matrix([[x1],[x2],[x3]])
print(X)
A = sym.Matrix([[1,0,1], [0,1,1],[1,1,0]])
print(A)
print(X.T*A*X)
D = sym.diff(X.T*A*X , X)
print(D.shape)
print(D)
print(D[0])
print(D[0][0])
print(D[0][0][0])
print(D[0][0][0][0])

代码输出结果 #

Matrix([[x1], [x2], [x3]])
Matrix([[1, 0, 1], [0, 1, 1], [1, 1, 0]])
Matrix([[x1*(x1 + x3) + x2*(x2 + x3) + x3*(x1 + x2)]])
(3, 1, 1, 1)
[[[[2*x1 + 2*x3]]], [[[2*x2 + 2*x3]]], [[[2*x1 + 2*x2]]]]
[[[2*x1 + 2*x3]]]
[[2*x1 + 2*x3]]
[2*x1 + 2*x3]
2*x1 + 2*x3

预期结果 #

[[2*x1 + 2*x3], [2*x2 + 2*x3], [2*x1 + 2*x2]]

原因说明 #

SymPy的diff函数处理矩阵对矩阵求导时，会严格保留原矩阵的维度嵌套关系：

• 首先，X^T A X的计算结果是一个**(1,1)的标量矩阵**，而非单纯的标量
• 当对列向量X（形状(3,1)）求导时，函数会对这个(1,1)矩阵里的唯一元素，逐个对X中的每个元素求导
• 每个求导结果都会保留原(1,1)矩阵的维度结构，所以每个导数都是(1,1)的矩阵
• 最终3个(1,1)矩阵嵌套起来，就形成了形状为(3,1,1,1)的4维数组

如果要得到预期的(3,1)形状矩阵，可以提取内层元素重新构造：

# 提取核心导数结果并生成目标矩阵
target_matrix = sym.Matrix([[D[i][0][0][0]] for i in range(3)])
print(target_matrix)

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