Fit quadratic function to data using MSE

Question

So my idea is (borrowed from neural network folks) that if I have data set D, I can fit a quadratic curve to it by first calculating the derivative of the error with respect to the parameters (a, b, and c), and then do a small update that minimizes the error. My problem is, the following code doesn't actually manage fit the curve. For linear stuff, a similar approach works, but quadratic seems to fail for some reason. Can you see what I have done wrong (either assumption or just implementation error)

EDIT: The question was not specific enough: This following code will not deal with bias in data very well. For some reason it updates a and b parameters somehow that c gets left behind. This method is similar to robotics (find path using Jacobians) and neural networks (find parameters based on error) so it is not unreasonable algorithm, now the question is, why this specific implementation does not produce results that I would expect.

In the following Python code, I use math as m, and MSE is a function that calculates Mean Squared Error between two arrays. Other than that, the code is self contained

Code:

def quadraticRegression(data, dErr):
    a = 1 #Starting values
    b = 1
    c = 1

    a_momentum = 0 #Momentum to counter steady state error
    b_momentum = 0
    c_momentum = 0

    estimate = [a*x**2 + b*x + c for x in range(len(data))] #Estimate curve
    error = MSE(data, estimate) #Get errors 'n stuff
    errorOld = 0

    lr = 0.0000000001 #learning rate
    while abs(error - errorOld) > dErr:
        #Fit a (dE/da)
        deda = sum([ 2*x**2 * (a*x**2 + b*x + c - data[x]) for x in range(len(data)) ])/len(data)
        correction = deda*lr
        a_momentum = (a_momentum)*0.99 + correction*0.1 #0.99 is to slow down momentum when correction speed changes
        a = a - correction - a_momentum
    
        #fit b (dE/db)
        dedb = sum([ 2*x*(a*x**2 + b*x + c - data[x]) for x in range(len(data))])/len(data)
        correction = dedb*lr
        b_momentum = (b_momentum)*0.99 + correction*0.1
        b = b - correction - b_momentum

        #fit c (dE/dc)
        dedc = sum([ 2*(a*x**2 + b*x + c - data[x]) for x in range(len(data))])/len(data)
        correction = dedc*lr
        c_momentum = (c_momentum)*0.99 + correction*0.1
        c = c - correction - c_momentum

        #Update model and find errors
        estimate = [a*x**2 +b*x + c for x in range(len(data))]
        errorOld = error
        print(error)
        error = MSE(data, estimate)

    return a, b, c, error

Answer 1

To me it looks like your code works totally correct! At least algorithm is correct. I've changed your code to use numpy for fast computation instead of pure Python. Also I've configured some params a bit, like changed momentum and learning rate, also implemented MSE.

Then I used matplotlib to draw plot animation. Finally on animation it looks like that your regression actually tries to fit the curve to data. Although at the last iteration for fitting sin(x) for x in [0; 2 * pi] it looks like linear approximation, but still close to data points as much as possible for quadratic curve to be. But for sin(x) for x in [0; pi] it looks like ideal approximation (it starts fitting from around 12-th iteration).

i-th frame of animation just does regression with dErr = 0.7 ** (i + 15).

My animation is a bit slow to just run script, but if you add param save like this python script.py save it will render/save to line.gif plot drawing animation. If you run the script without params it will animation on your PC screen plotting/fitting in real time.

Full code goes after graphics, code needs installing some python modules by running once python -m pip install numpy matplotlib.

Next is sin(x) for x in (0, pi):

Next is sin(x) for x in (0, 2 * pi):

Next is abs(x) for x in (-1, 1):

# Needs: python -m pip install numpy matplotlib
import math, sys
import numpy as np, matplotlib.pyplot as plt, matplotlib.animation as animation
from matplotlib.animation import FuncAnimation


x_range = (0., math.pi, 0.1) # (xmin, xmax, xstep)
y_range = (-0.2, 1.2) # (ymin, ymax)
num_iterations = 50

def f(x):
    return np.sin(x)

def derr(iteration):
    return 0.7 ** (iteration + 15)
    
    
def MSE(a, b):
    return (np.abs(np.array(a) - np.array(b)) ** 2).mean()

def quadraticRegression(*, x, data, dErr):
    x, data = np.array(x), np.array(data)
    assert x.size == data.size, (x.size, data.size)

    a = 1 #Starting values
    b = 1
    c = 1

    a_momentum = 0.1 #Momentum to counter steady state error
    b_momentum = 0.1
    c_momentum = 0.1

    estimate = a*x**2 + b*x + c #Estimate curve
    error = MSE(data, estimate) #Get errors 'n stuff
    errorOld = 0.

    lr = 10. ** -4 #learning rate
    while abs(error - errorOld) > dErr:
        #Fit a (dE/da)
        deda = np.sum(2*x**2 * (a*x**2 + b*x + c - data))/len(data)
        correction = deda*lr
        a_momentum = (a_momentum)*0.99 + correction*0.1 #0.99 is to slow down momentum when correction speed changes
        a = a - correction - a_momentum
    
        #fit b (dE/db)
        dedb = np.sum(2*x*(a*x**2 + b*x + c - data))/len(data)
        correction = dedb*lr
        b_momentum = (b_momentum)*0.99 + correction*0.1
        b = b - correction - b_momentum

        #fit c (dE/dc)
        dedc = np.sum(2*(a*x**2 + b*x + c - data))/len(data)
        correction = dedc*lr
        c_momentum = (c_momentum)*0.99 + correction*0.1
        c = c - correction - c_momentum

        #Update model and find errors
        estimate = a*x**2 +b*x + c
        errorOld = error
        #print(error)
        error = MSE(data, estimate)

    return a, b, c, error


        
fig, ax = plt.subplots()
fig.set_tight_layout(True)

x = np.arange(x_range[0], x_range[1], x_range[2])
#ax.scatter(x, x + np.random.normal(0, 3.0, len(x)))
line0, line1 = None, None

do_save = len(sys.argv) > 1 and sys.argv[1] == 'save'

def g(x, derr):
    a, b, c, error = quadraticRegression(x = x, data = f(x), dErr = derr)
    return a * x ** 2 + b * x + c
    
def dummy(x):
    return np.ones_like(x, dtype = np.float64) * 100.

def update(i):
    global line0, line1
    
    de = derr(i)
    
    if line0 is None:
        assert line1 is None
        line0, = ax.plot(x, f(x), 'r-', linewidth=2)
        line1, = ax.plot(x, g(x, de), 'r-', linewidth=2, color = 'blue')
        ax.set_ylim(y_range[0], y_range[1])
        
    if do_save:
        sys.stdout.write(str(i) + ' ')
        sys.stdout.flush()

    label = 'iter {0} derr {1}'.format(i, round(de, math.ceil(-math.log(de) / math.log(10)) + 2))
    line1.set_ydata(g(x, de))
    ax.set_xlabel(label)
    return line1, ax

if __name__ == '__main__':
    anim = FuncAnimation(fig, update, frames = np.arange(0, num_iterations), interval = 200)
    if do_save:
        anim.save('line.gif', dpi = 200, writer = 'imagemagick')
    else:
        plt.show()