Curve Fitting with Least Squares

Some time ago, I created a script to fit a 2D function. At the time I was dealing with noisy functions and was a bit less knowledgeable than today.

Engineering Way

I know that I wanted to create a fit function without actually worrying about the method. I created a github script around 8 years ago. The code with some minimal modifications is as follows.

import numpy as np
from scipy.optimize import curve_fit
np.random.seed(0)

# some function
def func(x, a, b, c):
    return a*np.sin(x[0])+b*np.cos(x[1])+c

# creating fake data
limits = [0, 2*np.pi, 0, 2*np.pi]  # [x1_min, x1_max, x2_min, x2_max]
side_x = np.linspace(limits[0], limits[1], 100)
side_y = np.linspace(limits[2], limits[3], 100)
X1, X2 = np.meshgrid(side_x, side_y)
size = X1.shape
x1_1d = X1.reshape((1, np.prod(size)))
x2_1d = X2.reshape((1, np.prod(size)))

xdata = np.vstack((x1_1d, x2_1d))
# original value or the result we expect to obtain
original = (3, 1, 0.5)
z = func(xdata, *original)
Z = z.reshape(size)
z_noise = z + .2*np.random.randn(len(z))
Z_noise = z_noise.reshape(size)
ydata = z_noise

# fit the noisy data to the curve
popt, pcov = curve_fit(func, xdata, ydata)
print("original: {}\nfitted: {}".format(original, popt))
# result
# original: (3, 1, 0.5)
# fitted: [3.00411704 0.99672656 0.49634599]

Thus, we created fake values that we want to fit in a function. In this example we use the same function as the generating function and the function we want to fit, but this is not necessary, we can have data coming from somewhere else.

Also, as you can see, the function curve_fit gives two values, popt and pcov. The first one is the response and the second one is the uncertainty quantification (or a covariance matrix on the error)

z_fit = func(xdata, *popt)
Z_fit = z_fit.reshape(size)
import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title("Real Function")
plt.pcolormesh(X1, X2, Z)
plt.axis(limits)
plt.colorbar()
plt.subplot(1, 3, 2)
plt.title("Function w/ Noise")
plt.pcolormesh(X1, X2, Z_noise)
plt.axis(limits)
plt.colorbar()
plt.subplot(1, 3, 3)
plt.title("Fitted Function from Noisy One")
plt.pcolormesh(X1, X2, Z_fit)
plt.axis(limits)
plt.colorbar()
plt.show()

The result is going to be the plot you see at the beginning (actually with other colors because jet is not the standard palette in matplotlib anymore).

More In Depth

The function curve_fit will fit the value of data to a given function. The documentation says they use non-linear least squares. We can think of this as solving for the Root Mean Squared Error.

$$ \begin{equation} \text{RME}=\sum_i (y_i - f(x_i))^2 \end{equation} $$

In our case, we are solving for an equation where we want to find the values of polynomials. Of course, the data of $X$ seems tricky because it goes through $\sin$ and $\cos$ but we can actually define.

$$ \bar{X} = \begin{bmatrix} | & | & | \\ \sin X_0 & \cos X_1 & \ \textbf{1} \ \\ | & | & | \\ \end{bmatrix} $$

Since we have a linear function, we can rewrite equation (1) as follows. $$ \text{RME}= |Y - X\beta|^2 = (Y - X\beta)^\mathsf{T}(Y - X\beta) $$

We then take the derivative of this equation w.r.t. $X$, set it to zero and solve for $\beta$. The Wikipedia Article will point us to the following solution, which is the standard Linear Least Squares.

$$ \begin{equation} \beta=\left (X^\mathsf{T} X \right )^{-1} X^\mathsf{T} Y \end{equation} $$

We have to replace for our variable $\bar{X}$ instead of $X$. Now, we can code the small bit of code given the information we have.

X_bar = np.c_[np.sin(xdata[0]), np.cos(xdata[1]), np.ones(xdata.shape[1])]
# writting as in the equation
result = np.linalg.inv(X_bar.T@X_bar)@X_bar.T@ydata
# rewrite the function in order not to use inverse (because it's a costly and unstable function)
result = np.linalg.solve(X_bar.T@X_bar, X_bar.T@ydata)
# rewrite the function because it exists actually a least squares function in numpy
result = np.linalg.lstsq(X_bar,ydata)[0]
# all results give [3.00411704 0.99672656 0.49634599] as expected (the same value given by curve_fit)

I wrote 3 different times the same thing just in a different form. The first time, I respected equation (1). However this way is computationally more expensive and less reliable that the two other ways, in which the last is the fastest, since it explicitly solves the linear least squares problem.

As an interesting note, you will find that numpy and scipy have two really different meanings for the function least_squares, the first one is to find the values given the variables (as we saw in our example), the second is to fit a given equation or function.

Further reading

• Curve Fitting Least Squares
• Curve Fitting with Python
• Linear Least Squares
• Gilbert Strang’s book Introduction to Linear Algebra Section 4.3 Least Squares Approximations