如何使用 scipy.optimize.minimize(...) 找到 z = f(x, y) （如椭圆）的最优参数？

Question

我正在尝试学习如何使用 python / numpy / scipy 优化更高维度的数据拟合。在成功地使用 scipy 拟合 y = f(x) 形式的直线后，我尝试扩展逻辑以拟合 z = f(x, y) 形式的椭圆；两者都如下所示。我希望可以推广这种方法以适应更高维度的形状（即球体）。

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

## LINE EXAMPLE
npoints = 30
x = np.linspace(-5, 5, npoints)
m, b = 3, -4 # slope, y-intercept
initial_parameter_guess = (2.5, -1)
y = m * x + b # exact line
noise = np.random.uniform(-1, 1, size=x.size)
yn = y + noise # line with noise

def get_residuals(prms, x, y):
    """ """
    return (y - (prms[0] * x + prms[1]))**2

def f_error(prms, x, y):
    """ """
    resid = get_residuals(prms, x, y)
    return np.sum(resid)

result = minimize(f_error, x0=initial_parameter_guess, args=(x, yn))
# print(result)
yf = result.x[0] * x + result.x[1]

fig, ax = plt.subplots()
ax.scatter(x, yn, color='b', marker='.')
ax.plot(x, yf, color='r', alpha=0.5)
ax.grid(color='k', alpha=0.3, linestyle=':')
plt.show()
plt.close(fig)

将此逻辑应用于椭圆的情况，

## ELLIPSE EXAMPLE
npoints = 75
theta = np.random.uniform(0, 2*np.pi, size=npoints)
a, b = (3, 5)
initial_parameter_guess = (2.5, 6)
xnoise = np.random.uniform(-1, 1, size=theta.size)
ynoise = np.random.uniform(-1, 1, size=theta.size)
x = a**2 * np.cos(theta)
xn = x + xnoise
y = b**2 * np.sin(theta)
yn = y + ynoise

def get_residuals(prms, x, y):
    """ """
    return 1 - ((x/prms[0])**2 + (y/prms[1])**2)

def f_error(prms, x, y):
    """ """
    resid = get_residuals(prms, x, y)
    return np.sum(resid)

result = minimize(f_error, x0=initial_parameter_guess, args=(xn, yn))
# print(result)

通过scipy的minimize算法没有找到最优参数；以下输出显示为 print(result):

      fun: -4243.611573066522
 hess_inv: array([[41.44400248, 39.68101343],
       [39.68101343, 37.99343048]])
      jac: array([-1496.81719971,  2166.68896484])
  message: 'Desired error not necessarily achieved due to precision loss.'
     nfev: 174
      nit: 1
     njev: 42
   status: 2
  success: False
        x: array([-1.51640333,  2.93879038])

我看过另一个，它使用了最小二乘矩阵公式，但是这个例子对我来说有点难理解。我见过的几乎所有方法都是基于已发布 link 中的方法。我更喜欢使用 minimize 除非线性代数方法由于我目前不知道的原因而更好。

无论如何，我的上述方法是否可以 adapted/tweaked 行得通并且可以推广到更高维度？

Answer 1

代码有两个问题。不是最小化残差和，而是最小化残差平方和。其次，椭圆方程应定义为x=a*np.cos(theta)和y=b*np.sin(theta)

npoints = 75
theta = np.random.uniform(0, 2*np.pi, size=npoints)
a, b = (3, 5)
initial_parameter_guess = (2.5, 6)
xnoise = np.random.uniform(-1, 1, size=theta.size)
ynoise = np.random.uniform(-1, 1, size=theta.size)
x = a * np.cos(theta)
xn = x + xnoise
y = b * np.sin(theta)
yn = y + ynoise
def get_residuals(prms, x, y):
    """ """
    return 1 - ((x/prms[0])**2 + (y/prms[1])**2)

def f_error(prms, x, y):
    """ """
    resid = get_residuals(prms, x, y)
    return np.sum(np.square(resid))
result = minimize(f_error, x0=initial_parameter_guess,args=(xn, yn))

result
  fun: 5.85099318913613
 hess_inv: array([[ 0.07025572, -0.02902779],
       [-0.02902779,  0.12040811]])
      jac: array([-5.96046448e-08,  1.19209290e-07])
  message: 'Optimization terminated successfully.'
     nfev: 48
      nit: 10
     njev: 12
   status: 0
  success: True
        x: array([3.35248219, 5.13728323])

如何使用 scipy.optimize.minimize(...) 找到 z = f(x, y) （如椭圆）的最优参数？

How to use scipy.optimize.minimize(...) to find the optimal parameters of z = f(x, y) (like an ellipse)?

python

optimization

numpy

scipy

n-dimensional