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

Question

I am trying to learn how to optimize data fits in higher dimensions using python / numpy / scipy. Having successfully used scipy to fit a line of the form y = f(x), I tried extending the logic to fit an ellipse of the form z = f(x, y); both are shown below. I'm hoping this approach can be generalized to fit shapes in higher dimensions (ie, sphere).

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)

Applying this logic to the case of an ellipse,

## 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)

The minimize algorithm via scipy does not find the optimal parameters; the following output is shown with 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])

I have seen another solution to this problem, which use the matrix formulation of least squares, but the example is a bit difficult for me to follow. Almost all approaches I've seen are based on the approach in the posted link. I'd prefer using minimize unless the linear algebra approach was better for a reason that I am currently unaware of.

Anyways, can my approach above be adapted/tweaked in a way that will work and can be generalized for higher dimensions?

Answer 1

Two issues with the code . Instead of minimizing the sum of residuals , minimize the sum of squares of residuals. Secondly , elliptic equation should be defined as x=a*np.cos(theta) and 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])