`scipy.optimize.minimize` 中的 Jacobian 和 Hessian 输入

Posted 2023-03-05

【中文标题】`scipy.optimize.minimize` 中的 Jacobian 和 Hessian 输入

【英文标题】：Jacobian and Hessian inputs in `scipy.optimize.minimize`

【发布时间】：2017-04-29 10:59:37

【问题描述】：

我试图了解 Python 的 scipy.optimize.minimize 函数中的“dogleg”方法是如何工作的。我正在修改帮助页面底部的示例。

根据注释，dogleg 方法需要 Jacobian 和 Hessian 参数。为此，我使用numdifftools 包：

import numpy as np
from scipy.optimize import minimize
from numdifftools import Jacobian, Hessian

def fun(x,a):
    return (x[0] - 1)**2 + (x[1] - a)**2

x0 = np.array([2,0]) # initial guess
a = 2.5

res = minimize(fun, x0, args=(a), method='dogleg',
               jac=Jacobian(fun)([2,0]), hess=Hessian(fun)([2,0]))

print(res)

编辑：

如果我按照以下帖子的建议进行更改，

res = minimize(fun, x0, args=a, method='dogleg',
               jac=Jacobian(lambda x: fun(x,a)),
               hess=Hessian(lambda x: fun(x,a)))

我收到一个错误TypeError: <lambda>() takes 1 positional argument but 2 were given。我做错了什么？

在初始猜测x0 处计算雅可比和黑森矩阵也是正确的吗？

【参考方案1】：

我知道这是一个玩具示例，但我想指出，使用像 Jacobian 或 Hessian 这样的工具来计算导数而不是推导函数本身是相当昂贵的。例如用你的方法：

x0 = np.array([2, 0])
a = 2.5
%timeit minimize(fun, x0, args=(a,), method='dogleg', jac=fun_der, hess=fun_hess)
100 loops, best of 3: 13.6 ms per loop

但是你可以这样计算导数函数：

def fun_der(x, a):
    dx = 2 * (x[0] - 1)
    dy = 2 * (x[1] - a)
    return np.array([dx, dy]

def fun_hess(x, a):
    dx = 2
    dy = 2
    return np.diag([dx, dy])

%timeit minimize(fun, x0, args=(a,), method='dogleg', jac=fun_der, hess=fun_hess)
1000 loops, best of 3: 279 µs per loop

如您所见，这几乎快了 50 倍。它真的开始增加复杂的功能。因此，我总是尝试自己明确地导出函数，无论这可能有多困难。一个有趣的例子是基于内核的Inductive Matrix Completion 实现。

argmin --> sum((A - gamma_u(X) Z gamma_v(Y))**2 - lambda * ||Z||**2)
where gamma_u = (1/sqrt(m_x)) * [cos(UX), sin(UX)] and
gamma_v = (1/sqrt(m_y)) * [cos(VY), sin(VY)]
X.shape = n_x, p; Y.shape = n_y, q; U.shape = m_x, p; V.shape = m_y, q; Z.shape = 2m_x, 2m_y

与显式推导和利用这些函数相比，从这个方程计算梯度和粗麻布是极其不合理的。正如@bnaul 指出的那样，如果您的函数确实具有封闭形式的导数，那么您确实想要计算和使用它们。

【讨论】：

不知道numdifftools，但你不能计算一次jacf = Jacobian(lambda x: fun(x, 2.5))，然后minimize( ... jac=jacf, hess=hessf )？

@denis 我认为您的意思是将其指向 bnaul 我不使用 numdifftools 或这些方法。正如我所提到的，如果可能的话，我喜欢直接派生函数，以避免在每次迭代时重复计算它们的开销。

【参考方案2】：

该错误来自对Jacobian 和Hessian 的调用，而不是来自minimize。用Jacobian(lambda x: fun(x, a)) 替换Jacobian(fun) 和类似的Hessian 应该可以解决问题（因为现在被微分的函数只有一个向量参数）。

另外一件事：(a) 就是 a，如果您希望它是一个元组，请使用 (a,)。

import numpy as np
from scipy.optimize import minimize
from numdifftools import Jacobian, Hessian

def fun(x, a):
    return (x[0] - 1) **2 + (x[1] - a) **2

def fun_der(x, a):
    return Jacobian(lambda x: fun(x, a))(x).ravel()

def fun_hess(x, a):
    return Hessian(lambda x: fun(x, a))(x)

x0 = np.array([2, 0]) # initial guess
a = 2.5

res = minimize(fun, x0, args=(a,), method='dogleg', jac=fun_der, hess=fun_hess)
print(res)

【讨论】：

谢谢。如果我对Jacobian 和Hessian 进行这些更改，我会收到错误ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()。此外，如果我想指定超过 1 个参数，我是否需要使用 args=(a,b,c) 之类的东西？

在这种特殊情况下，您实际上想要传递一个可调用的雅可比函数，因此您应该省略([2, 0])（我编辑了我的回复以反映这一点）。

好的...如果我这样做，我会得到TypeError: <lambda>() takes 1 positional argument but 2 were given。然后，如果我尝试使用lambda x,a: fun(x, a) 作为Jacobian 和Hessian 输入，我会得到ValueError: incompatible dimensions.

再想一想：无论如何，你真的没有理由想要这样做。如果您没有封闭式导数，则应该使用不需要的方法；像这样用数值计算导数是非常低效的。

公平点，但我特别需要使用 Python 中的“dogleg”算法（需要 Jacobian 和 Hessian）来解决我的问题。我还不关心代码的效率，我只想知道scipy.optimize.minimize(method='dogleg') 函数是如何工作的。

【参考方案3】：

你可以使用 autograd 代替

import numpy as np
from scipy.optimize import minimize
from autograd import jacobian, hessian

def fun(x, a):
    return (x[0] - 1) **2 + (x[1] - a) **2

def fun_der(x, a):
    return jacobian(lambda x: fun(x, a))(x).ravel()

def fun_hess(x, a):
    return hessian(lambda x: fun(x, a))(x)

x0 = np.array([2, 0]) # initial guess
a = 2.5

res = minimize(fun, x0, args=(a,), method='dogleg', jac=fun_der, hess=fun_hess)
print(res)