how to implement least square polynomial with no built in methods using python?

Question

currently running into a problem solving this. The objective of the exercise given is to find a polynom of certian degree (the degree is given) from a dataset of points (that can be noist) and to best fit it using least sqaure method. I don't understand the steps that lead to solving the linear equations? what are the steps or should anyone provide such a python program that lead to the matrix that I put as an argument in my decomposition program? Note:I have a python program for cubic splines ,LU decomposition/Guassian decomposition. Thanks.

I tried to apply guassin / LU decomposition straight away on the dataset but I understand there are more steps to the solution... I donwt understand how cubic splines add to the mix either..

Edit: guassian elimintaion :

import numpy as np
import math

def swapRows(v,i,j):
    if len(v.shape) == 1:
        v[i],v[j] = v[j],v[i]
    else:
        v[[i,j],:] = v[[j,i],:]
def swapCols(v,i,j):
    v[:,[i,j]] = v[:,[j,i]]
def gaussPivot(a,b,tol=1.0e-12):
    n = len(b)
    # Set up scale factors
    s = np.zeros(n)
    for i in range(n):
        s[i] = max(np.abs(a[i,:]))
    for k in range(0,n-1):
        # Row interchange, if needed
        p = np.argmax(np.abs(a[k:n,k])/s[k:n]) + k
        if abs(a[p,k]) < tol: error.err('Matrix is singular')
        if p != k:
            swapRows(b,k,p)
            swapRows(s,k,p)
            swapRows(a,k,p)
        # Elimination
        for i in range(k+1,n):
            if a[i,k] != 0.0:
                lam = a[i,k]/a[k,k]
                a[i,k+1:n] = a[i,k+1:n] - lam*a[k,k+1:n]
                b[i] = b[i] - lam*b[k]
    if abs(a[n-1,n-1]) < tol: error.err('Matrix is singular')
    # Back substitution
    b[n-1] = b[n-1]/a[n-1,n-1]
    for k in range(n-2,-1,-1):
        b[k] = (b[k] - np.dot(a[k,k+1:n],b[k+1:n]))/a[k,k]
    return b
def polyFit(xData,yData,m):
    a = np.zeros((m+1,m+1))
    b = np.zeros(m+1)
    s = np.zeros(2*m+1)
    for i in range(len(xData)):
        temp = yData[i]
        for j in range(m+1):
            b[j] = b[j] + temp
            temp = temp*xData[i]
        temp = 1.0
        for j in range(2*m+1):
            s[j] = s[j] + temp
            temp = temp*xData[i]
    for i in range(m+1):
        for j in range(m+1):
            a[i,j] = s[i+j]
    return gaussPivot(a,b)
degree = 10 # can be any degree
polyFit(xData,yData,degree)    

I was under the impression the code above gets a dataset of points and a degree. The output should be coeefients of a polynom that fits those points but I have a grader that was provided by my proffesor , and after checking the grading the polynom that returns has a lrage error. After that I tried the following LU decomposition instead:

import numpy as np
def swapRows(v,i,j):
    if len(v.shape) == 1:
        v[i],v[j] = v[j],v[i]
    else:
        v[[i,j],:] = v[[j,i],:]
def swapCols(v,i,j):
    v[:,[i,j]] = v[:,[j,i]]
def LUdecomp(a,tol=1.0e-9):
    n = len(a)
    seq = np.array(range(n))
    # Set up scale factors
    s = np.zeros((n))
    for i in range(n):
        s[i] = max(abs(a[i,:]))
    for k in range(0,n-1):
        # Row interchange, if needed
        p = np.argmax(np.abs(a[k:n,k])/s[k:n]) + k
        if abs(a[p,k]) < tol: error.err('Matrix is singular')
        if p != k:
            swapRows(s,k,p)
            swapRows(a,k,p)
            swapRows(seq,k,p)
    # Elimination
        for i in range(k+1,n):
            if a[i,k] != 0.0:
                lam = a[i,k]/a[k,k]
                a[i,k+1:n] = a[i,k+1:n] - lam*a[k,k+1:n]
                a[i,k] = lam
        return a,seq
def LUsolve(a,b,seq):
    n = len(a)
    # Rearrange constant vector; store it in [x]
    x = b.copy()
    for i in range(n):
        x[i] = b[seq[i]]
    # Solution
    for k in range(1,n):
        x[k] = x[k] - np.dot(a[k,0:k],x[0:k])
    x[n-1] = x[n-1]/a[n-1,n-1]
    for k in range(n-2,-1,-1):
        x[k] = (x[k] - np.dot(a[k,k+1:n],x[k+1:n]))/a[k,k]
    return x

the results were a bit better but nowhere near what it should be Edit 2: I tried the chebyshev method suggested in the comments and came up with:

import numpy as np

def chebyshev_transform(x, n):
    """
    Transforms x-coordinates to Chebyshev coordinates
    """
    return np.cos(n * np.arccos(x))

def chebyshev_design_matrix(x, n):
    """
    Constructs the Chebyshev design matrix
    """
    x_cheb = chebyshev_transform(x, n)
    T = np.zeros((len(x), n+1))
    T[:,0] = 1
    T[:,1] = x_cheb
    for i in range(2, n+1):
        T[:,i] = 2 * x_cheb * T[:,i-1] - T[:,i-2]
    return T
degree =10
f = lambda x: np.cos(X)
xdata = np.linspace(-1,1,num=100)
ydata = np.array([f(i) for i in xdata])
M = chebyshev_design_matrix(xdata,degree)
D_x ,D_y = np.linalg.qr(M)
D_x, seq = LUdecomp(D_x)
A = LUsolve(D_x,D_y,seq)

• I can't use linalg.qr in my program , it was just for checking how it works.In addition , I didn't get the 'slow way' of the formula that were in the comment.

• The program cant get an x point that is not between -1 and 1 , is there any way around it , any normalizition?

Thanks a lot.

Answer 1

Hints:

You are probably asked for an unsophisticated method. If the degree of the polynomial remains low, you can use the straightforward approach below. For the sake of the explanation, I'll use a cubic model.

Assume that you want to fit your data to this polynomial, by observing that it seems to follow a cubic behavior:

ax³ + bx² + cx + d ~ y

[All x and y should be understood with an index i which is omitted for notational convenience.]

If there are more than four data points, you get an overdetermined system of equations, usually with no solution. The trick is to consider the error on the individual equations, e = ax³ + bx² + cx + d - y, and to minimize the total error. As the error is a signed number, negative errors would make minimization impossible. Instead, we minimize the sum of squared errors. (The sum of absolute errors is another option but it unfortunately leads to a much harder problem.)

Min(a, b, c, d) Σ(ax³ + bx² + cx + d - y)²

As the unknown parameters are unconstrained, it suffices to look for a stationary point, i.e. cancel the gradient of the total error. By differentiation on the unknowns a, b, c and d, we obtain

2Σ(ax³x³ + bx²x³ + cxx³ + dx³ - yx³) = 0
2Σ(ax³x² + bx²x² + cxx² + dx² - yx²) = 0
2Σ(ax³x  + bx²x  + cxx  + dx  - yx ) = 0
2Σ(ax³   + bx²   + cx   + d   - y  ) = 0

As you can recognize, this is a square linear system of equations.