fprime的Python fsolve形状错误？答案

【问题标题】：Python fsolve shape of fprime is wrong?fprime的Python fsolve形状错误？
【发布时间】：2018-04-11 00:18:01
【问题描述】：

以下内容再简单不过了，但我找不到如何使它工作......

from scipy.optimize import fsolve

def a(x):
    return x*x-x

def ap(x):
    return 2*x-1

#works
print(fsolve(a, 0.3))

#works
print(fsolve(a, 0.3, fprime=ap))

#works
print(fsolve(a, [0.3], fprime=ap))

#works
print(fsolve(a, [0.3, 0.7]))

#crashes
print(fsolve(a, [0.3, 0.7], fprime=ap))

当它崩溃时给出错误

TypeError: fsolve: there is a mismatch between the input and output shape of the 'fprime' argument 'ap'.Shape should be (2, 2) but it is (2,).

ap 的输出维度看起来肯定应该和输入一样。这怎么可能出错（以及如何解决）？

我认为一些反对者忽略了这个问题的微妙之处，所以这里更深入地解释了我为什么感到困惑：

似乎 scipy 将 a 解释为一个变量的函数，[0.3,0.7] 开始估计 fsolve(a, [0.3, 0.7]) 中的两个根，但在 fsolve(a, [0.3, 0.7], fprime=ap) 中，它将 a 解释为带有 [.3,.7] 的两个变量的函数是单根的估计。根据文档（https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html），第二个参数

x0 : ndarray
    The starting estimate for the roots of func(x) = 0.

说它正在寻找根的估计（复数）。但是它在给定fprime 的情况下的行为听起来像是x0 被解释为对单个根的估计。

【问题讨论】：

你认为你会在下面的代码中得到什么输出：ap([0.3, 0.7])
scipy 要求 ap(x) 的 2x2 输出，因为这是 jacobian matrix 的形式。
@eyllanesc 我不希望 scipy 传入一个列表，而是一个数组，我希望 ap(np.array([0.3, 0.7])) == np.array([ap(0.3), ap(0.7)]) 确实如此。如果要允许ap 接收列表，我会将其定义为return 2*np.array(x)-1，但这并不能解决形状不匹配问题。
能否解释一下为什么这是一个糟糕的问题？对我来说，倒数第二个 fsolve(a, [0.3, 0.7]) 有效但 fsolve(a, [0.3, 0.7], fprime=ap) 无效的原因似乎很重要。
另外@sascha我意识到两个变量的函数的雅可比是一个2x2矩阵，但是一个变量的函数的雅可比也是一个变量的函数，所以我的问题是为什么scipy将a 解释为一个变量的函数，[.3,.7] 在fsolve(a, [0.3, 0.7]) 中开始估计两个根，但在fsolve(a, [0.3, 0.7], fprime=ap) 中它将a 解释为两个变量的函数，[.3,.7] 是单个根的估计?

标签： python scipy

【解决方案1】：

我怀疑您希望fsolve(a, [0.3, 0.7], fprime=ap) 返回标量问题的两个解决方案。 fsolve 不能那样工作。

当您调用fsolve(a, x0, fprime=ap) 时，fsolve 函数会根据x0 的形状推断问题的维度。如果x0 是一个标量，它期望a 接受一个标量，而fprime 必须接受一个标量并返回一个标量（或1x1 数组）。如果 x0 是长度为 2 的序列（如在您的示例中不起作用），fsolve 期望 a 接受长度为 2 的数组作为其 x 参数并返回长度为 2 的序列。它还期望fprime 接受一个长度为 2 的数组并返回一个形状为 2x2 的数组（雅可比矩阵）。

【讨论】：

确实，scipy.optimize.fsolve 的文档说第二个参数x0 是The starting estimate for the roots of func(x) = 0. 对我来说，这听起来像是它想要对所有根进行估计。事实上，fsolve(a, np.array([0.3, 0.7])) 返回[ -2.76444148e-16 1.00000000e+00] ，所以这似乎是正确的。为什么当我给它 fprime 时行为会改变？
如果你调用a(np.array([0.3, 0.7]))，函数返回array([-0.21, -0.21])，所以它是一个完美的两个值返回两个值的函数。这种函数的雅可比行列式是一个形状为 2x2 的二维矩阵（或者，在 numpy 术语中，一个形状为 (2, 2) 的二维数组）。这是ap 必须返回的矩阵，它才能与fsolve 保持一致。使用您的函数，ap(np.array([0.3, 0.7]) 返回array([-0.4, 0.4])。这是a(x) 的雅可比矩阵的错误形状。（如果你不使用参数fprime，fsolve 以数字方式计算雅可比行列式。）

【解决方案2】：

感谢所有解释我为什么会遇到我现在的行为的人。但实际上没有人告诉我如何解决这个问题，所以我想出了这个答案，所以我正在写这个答案。

混淆源于 scipy 文档介于错误和误导之间的事实。

文档说：

在给定初始估计的情况下，返回由 func(x) = 0 定义的（非线性）方程的根。

建议对于标量函数，它可能会返回一个根数组。实际行为更接近

在给定初始估计的情况下，返回由 func(x) = 0 定义的（非线性）方程的根。

这会使做什么变得显而易见。为了找到等式的多个根，我们可以让 scipy（一如既往）将a(x) 解释为n 变量的函数，其中n 是x 的长度。要找到多个根，让 a 向量化，即

a([x1, x2, ..., xn]) = [a(x1), a(x2), ..., a(xn)].

这个函数的雅可比是一个对角矩阵

diag([ap(x1), ap(x2), ..., ap(xn)]).

因此我的问题的答案是：

改变

print(fsolve(a, [0.3, 0.7], fprime=ap))

到

print(fsolve(a, [0.3, 0.7], fprime=lambda x: np.diag(ap(x))))
# correctly outputs
# [ -3.69321143e-16   1.00000000e+00]

或者简单地循环 [0.3, 0.7] 调用 fsolve 每次初始猜测一次。

【讨论】：