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

...my real problem involves having more equations than unknowns, so at the end I want a solution that minimizes the errors of all this equations

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

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

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

它有 3 个未知数和 4 个方程。那么您的 objective 函数将是：

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是单位矩阵，你得到最小二乘objective函数。然后，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])

scipy 尽量减少找不到解决方案

scipy minimize not finding solution

python

scipy

scipy-optimize-minimize