如何为Sympy添加支持求导的自定义矩阵函数？

针对你提出的三个问题，以下是具体解答：

1. 子类化MatrixExpr/Expr时，_eval_derivative方法的实现逻辑 #

_eval_derivative的核心作用是接收一个求导变量（Symbol或MatrixSymbol元素），返回自定义函数对该变量的导数表达式。具体逻辑需贴合函数的数学定义：

• 若为标量输出的矩阵函数（如谱范数）：需根据其数学导数规则推导表达式，比如谱范数（最大奇异值）的导数依赖于对应的左/右奇异向量，可结合SVD分解构建导数公式。
• 若为矩阵输出的函数（如特征向量）：要返回与原函数同形状的MatrixExpr，需基于特征值扰动理论推导，比如非简并特征向量的导数涉及特征值差与投影矩阵的组合。
• 方法内部可调用Sympy内置的微分工具（如diff）、矩阵分解（如SingularValueDecomposition）或符号特征值/向量函数，构建符号化的导数结果。
• 注意返回值类型：标量函数返回Expr，矩阵函数返回MatrixExpr。

2. 还需定义或配置的内容 #

除了_eval_derivative，还有几个关键部分需要实现：

• __new__方法：用于输入验证（确保输入是MatrixExpr）和实例缓存（Sympy的标准优化方式，避免重复创建相同参数的对象）。
• 属性封装：比如定义arg属性存储输入矩阵，方便在导数方法中快速引用。
• 形状定义（仅MatrixExpr子类）：实现shape()方法，返回输出矩阵的维度（如特征向量是列向量，形状为(n,1)）。
• LaTeX输出（可选）：实现_repr_latex_方法，让函数在Jupyter或LaTeX环境中显示标准数学符号（比如谱范数的双竖线）。
• 多值函数的歧义处理：比如特征向量是多值的，需在类中明确规则（如取最大特征值对应的归一化向量），避免导数结果出现歧义。
• 简化逻辑（可选）：实现_simplify方法，或注册自定义简化规则，让导数结果自动简化。

3. 可参考的实现示例 #

以下是两个简化的自定义矩阵函数示例，涵盖标量输出（谱范数）和矩阵输出（主特征向量）的场景：

示例1：自定义谱范数（标量输出，继承Expr） #

from sympy import Expr, MatrixSymbol, diff, SingularValueDecomposition
from sympy.matrices.expressions import MatrixExpr

class SpectralNorm(Expr):
    def __new__(cls, mat):
        # 验证输入类型
        if not isinstance(mat, MatrixExpr):
            raise TypeError("Input must be a MatrixExpr")
        # 缓存实例，Sympy标准做法
        return super().__new__(cls, mat)
    
    @property
    def arg(self):
        return self.args[0]
    
    def _eval_derivative(self, var):
        mat = self.arg
        # 用SVD分解推导谱范数导数（假设最大奇异值非简并）
        svd = SingularValueDecomposition(mat)
        sigma_max = svd.S[0, 0]
        # 返回最大奇异值对变量的导数，等价于谱范数的导数
        return diff(sigma_max, var)
    
    def _repr_latex_(self):
        return r"\|%s\|_2" % self.arg._repr_latex_()

# 测试
A = MatrixSymbol('A', 2, 2)
norm_A = SpectralNorm(A)
print(diff(norm_A, A[0,0]))

示例2：自定义主特征向量（矩阵输出，继承MatrixExpr） #

from sympy import MatrixSymbol, diff, eigenvects
from sympy.matrices.expressions import MatrixExpr

class LeadingEigenvector(MatrixExpr):
    def __new__(cls, mat):
        if not isinstance(mat, MatrixExpr):
            raise TypeError("Input must be a MatrixExpr")
        return super().__new__(cls, mat)
    
    @property
    def arg(self):
        return self.args[0]
    
    def shape(self):
        # 返回列向量形状
        return (self.arg.shape[0], 1)
    
    def _eval_derivative(self, var):
        mat = self.arg
        # 获取最大特征值对应的特征向量（简化版，假设非简并）
        ev_pairs = eigenvects(mat)
        max_ev_pair = max(ev_pairs, key=lambda x: x[0])
        eigen_vec = max_ev_pair[2][0]
        # 返回特征向量对变量的导数
        return diff(eigen_vec, var)
    
    def _repr_latex_(self):
        return r"\mathbf{v}_{\text{max}}(%s)" % self.arg._repr_latex_()

# 测试
A = MatrixSymbol('A', 2, 2)
vec_A = LeadingEigenvector(A)
print(diff(vec_A, A[1,0]))

注：实际生产环境中，需补充非简并性判断、归一化处理等逻辑，避免符号计算中的异常。

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