scipy最小化没有找到解决方案

Question

我正在尝试使用 scipy.minimize 解决一组方程，但是我没有得到令人满意的结果，所以也许我弄错了。 我想求解以下方程组。

12.25 * (x + y * 2.2 + z * 4.84) - 8.17437483750257 = 0
12.25 * (x + y * 3.1 + z * 9.61) - 21.9317236606432 = 0
12.25 * (x + y * 4   + z * 16)   - 107.574834524443 = 0

使用 Wolfram Alpha 我得到答案

x=22.626570068753, y=-17.950683342597, z=3.6223614029055

这确实解决了方程组，给出的残差为

9.407585821463726e-12

现在我使用 scipy.minimize：

import numpy as np
from scipy.optimize import fsolve
from scipy.optimize import minimize

def my_func(p):
    points = [8.17437483750257, 21.9317236606432, 107.574834524443]
    h1 = abs(12.25 * (p[0] + p[1] * 2.2 + p[2] * 4.84) - points[0])
    h2 = abs(12.25 * (p[0] + p[1] * 3.1 + p[2] * 9.61) - points[1])
    h3 = abs(12.25 * (p[0] + p[1] * 4   + p[2] * 16)   - points[2])
    return h1 + h2 + h3

ini = np.array([22, -15, 5])   # Initial points close to solution
res = minimize(my_func, ini)
print(res)



  fun: 1.4196640741924451
  hess_inv: array([[ 20.79329103, -14.63447889,   2.36145776],
   [-14.63447889,  10.30037625,  -1.66214485],
   [  2.36145776,  -1.66214485,   0.26822135]])
  jac: array([ 12.25      ,  60.02499545, 254.43249989])
  message: 'Desired error not necessarily achieved due to precision loss.'
  nfev: 261
  nit: 8
  njev: 64
  status: 2
  success: False
    x: array([ 21.39197235, -17.08623345,   3.48344393])

首先，它说success=False，其次它找到不是最优的解决方案。

为什么初始值接近最优解却无法找到这些解。

优化器的定义有问题吗？

尝试运行它并给出 [0,0,0] 的初始值，但结果很糟糕

ini = np.array([0, 0, 0])   # Initial points close to solution
res = minimize(my_func, ini)
print(res)

fun: 73.66496363902732
 hess_inv: array([[ 0.98461683, -0.04223651, -0.1207056 ],
       [-0.04223651,  0.88596592, -0.31885642],
       [-0.1207056 , -0.31885642,  0.13448927]])
      jac: array([ 12.25      ,  15.92499924, -18.98750019])
  message: 'Desired error not necessarily achieved due to precision loss.'
     nfev: 164
      nit: 1
     njev: 40
   status: 2
  success: False
        x: array([0.02901304, 0.08994042, 0.29448233])

注意：我不想使用fsolve 来查找解决方案，而是使用minimize。 原因是我的真正问题涉及到的方程多于未知数，所以最后我想要一个能够最大限度地减少所有这些方程的误差的解决方案。 然而，由于它没有给出好的结果，我想首先测试一个存在精确解决方案的简单问题。但即使在这种情况下，它也不起作用。 一旦我使它适用于这个问题，我将扩展它添加更多方程。

Answer 1

...我真正的问题是方程比未知数多，所以最后我想要一个解决方案，使所有这些方程的误差最小化

这听起来很像在广义矩量法 (GMM) 中解决的问题，其中方程比未知数还多。

这类问题通常使用最小二乘法来解决。假设您的整个系统如下所示：

h1(x, y, z) = 0
h2(x, y, z) = 0
h3(x, y, z) = 0
h4(x, y, z) = 0

它有 3 个未知数和 4 个方程。那么你的目标函数将是：

F(x, y, z) = H(x, y, z)' * W * H(x, y, z)

• H(x, y, z)是上面所有hj(x, y, z)的向量
• H(x, y, z)' 是它的转置
• W是权重矩阵

如果W 是单位矩阵，则得到最小二乘目标函数。那么F(x, y, z)是一个二次型（基本上是多维抛物线），应该很容易优化，因为它是凸的和光滑的。

您的代码使用像 h1 = abs(12.25 * (p[0] + p[1] * 2.2 + p[2] * 4.84) - points[0]) 这样的绝对值，但在原点附近可能难以区分 abs，但这正是您的最佳选择所在，因为您本质上希望 h1 等于零。

您可以通过平方误差来近似绝对值函数：

h1 =(12.25 * (p[0] + p[1] * 2.2 + p[2] * 4.84) - points[0])**2

这产生了与 GMM（或最小二乘法）基本相同的方法，并为您提供了一个易于优化的函数，因为正方形在原点附近是平滑的。

Answer 2

优化问题（和求解器）通常受益于表现良好（平滑）的“优化表面”。当您使用 abs 函数时，它会创建一个“不稳定”的曲面，其中的点是导数不连续的。

如果不使用abs 而使用二次函数（具有相同的效果），您将获得接近预期的解决方案。只需将my_func 更改为：

def my_func(p):
    points = [8.17437483750257, 21.9317236606432, 107.574834524443]
    h1 = (12.25 * (p[0] + p[1] * 2.2 + p[2] * 4.84) - points[0])**2
    h2 = (12.25 * (p[0] + p[1] * 3.1 + p[2] * 9.61) - points[1])**2
    h3 = (12.25 * (p[0] + p[1] * 4   + p[2] * 16)   - points[2])**2
    return h1 + h2 + h3

我得到的是：

      fun: 8.437863292878727e-10
 hess_inv: array([[ 0.64753863, -0.43474506,  0.06909179],
       [-0.43474506,  0.29487798, -0.04722923],
       [ 0.06909179, -0.04722923,  0.00761762]])
      jac: array([-6.26698693e-12,  6.22490927e-10, -5.11716516e-10])
  message: 'Optimization terminated successfully.'
     nfev: 55
      nit: 7
     njev: 11
   status: 0
  success: True
        x: array([ 22.62653789, -17.95066124,   3.62235782])