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?
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])