如何计算 Sympy 表达式的梯度和 Hessian 矩阵

Question

假设我们有一个 Sympy 表达式：

import sympy as sp 

a = sp.Symbol('a', real=True)
b = sp.Symbol('b', real=True)
c = sp.Symbol('c', real=True)

f = 0.2*a**3 - 0.1*a**2*b - 0.9*a**2*c + 0.5*a**2 - 0.1*a*b**2 - 0.4*a*b*c + 0.5*a*b - 
0.6*a*c**2 - 0.5*a*c - 0.6*a - 0.1*b**3 + 0.4*b**2*c - 0.4*b**2 + 0.7*b*c**2 - 0.4*b*c + 
0.9*b + 0.09*c**3 - 0.6*c**2 + 0.4*c + 0.22

如何获取可在任意点求值的梯度表达式列表以及可再次求值为任意点的常量矩阵的 Hessian 矩阵表达式列表列表。

Answer 1

我确信有一些更聪明、更优雅的解决方案，但为了以防万一这可以节省任何人的时间，下面是我手动执行的操作。我为渐变定义了一个类：

class GradF:
    def __init__(self, function_expression, variables: list):
        self.function_expression = function_expression
        self.variables = variables
        self.exprs = []

        for i in range(len(self.variables)):
            self.exprs.append(diff(self.function_expression, self.variables[i]))

        self.lambdas = [sp.lambdify(self.variables, self.exprs[i], "numpy") for i in 
                       range(len(self.variables))]

    def __call__(self,x):
        assert len(x) == len(self.variables), "number of variables must match number of 
                                               arguments"
        values_list = [self.lambdas[i](*x) for i in range(len(self.lambdas))]

        return np.array(values_list).astype(np.float64)

    def latex(self):
        return sp.print_latex(self.exprs)

还有黑森人的课程：

class HessF:

    def __init__(self, gradf: GradF):
        self.gradf = gradf
        self.variables = gradf.variables

       self.hess = [[0 for _ in range(len(self.variables))] for _ in range(len(self.variables))]

       col = 0
       row = 0
       while col < len(self.variables):
           while row < len(self.variables):
               self.hess[row][col] = self.gradf.exprs[col].diff(self.variables[row])
               row += 1

           col += 1
           row = col

        self.hess = self.hess + np.transpose(self.hess) - np.diag(np.diag(self.hess))


    def __call__(self, x: np.ndarray):

        vals = np.zeros((len(self.variables), len(self.variables)))
        col = 0
        row = 0
        while col < len(self.variables):
            while row < len(self.variables):
                vals[row][col] = self.hess[row][col].subs(zip(self.variables, x))
                row += 1
            col += 1
            row = col

        vals = vals + np.transpose(vals) - np.diag(np.diag(vals))

        return vals


    def latex(self):
        return sp.print_latex(self.hess)

然后对于任何符号列表和函数，我们都可以获得 Gradient 和 Hessian 对象：

variables = [a,b,c]
function = f

gradient = GradF(function, variables)
hessian = HessF(gradient)

对于表示为 Numpy 数组的任何点，例如

x = np.array([0.1,3,-4])

我们可以得到这些值：

f_prime = gradient(x)
f_hessian = hessian(x)