Polyfit - Executable Engineering: Interactive Numerical Methods in Python

import numpy as np
import plotly.graph_objects as go
import plotly.colors as colors

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

Introduction to Polynomial Fitting¶

Let’s return to our polynomial fitting armed with our new tool, and use it to separate the order of the polynomial from the number of data points.

An $n$ th degree polynomial,

y(x) = a_n x^n + a_{n-1} x^{n-1} \dots a_2 x^2 + a_1 x +a_0

(1)

can be applied to $m$ data points,

y(x_i) = a_n x_i^n + a_{n-1} x_i^{n-1} \dots a_2 x_i^2 + a_1 x_i +a_0 = y_i

(2)

to generate an $m \times n$ matrix, multiplied by an $n$ vector of polynomial coefficients to equal an $m$ vector of data:

\begin{aligned} \begin{bmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^n \\ 1 & x_2 & x_2^2 & \cdots & x_2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_m & x_m^2 & \cdots & x_m^n \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_m \end{bmatrix} \end{aligned}

(3)

Example: polynomial fitting¶

Let’s generate some noisy data from a known 4th order polynomial. We’ll add normally distributed noise to the true function values.

# Generate x values
x = np.linspace(-5, 5, 100)

# Define the true function
def true_function(x):
  return 3 * x**4 - 2 * x**2 + x - 9

# Calculate the true y values
y_true = true_function(x)

# Generate noisy data
np.random.seed(0)  # For reproducibility
noise = np.random.normal(0, 100, 100)
x_data = np.linspace(-5, 5, 100)
y_data = true_function(x_data) + noise

# Plot the true curve and the data
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y_true, mode='lines', line=dict(dash='dash', color='blue'), name='True Function'))
fig.add_trace(go.Scatter(x=x_data, y=y_data, mode='markers', marker=dict(size=5, color='red'), name='Noisy Data'))

fig.update_layout(title='True Function and Noisy Data', xaxis_title='x', yaxis_title='y')
fig.show()

Fitting a Cubic Polynomial¶

We can solve for the polynomial coefficients directly using the pseudoinverse of the Vandermonde matrix, or using numpy.polyfit. Let’s compare them and plot the result for a 3rd degree (cubic) polynomial.

# Generate the Vandermonde matrix for a cubic polynomial
n = 3  # Degree of the polynomial
X = np.vander(x_data, n + 1)

# Calculate the coefficients using the pseudoinverse
coefficients = np.linalg.pinv(X) @ y_data

coeffs = np.polyfit(x_data, y_data, n)
print('Coefficients calculated manually:', coefficients)
print('Coefficients calculated with polyfit:', coeffs)

# Generate y values for the fitted polynomial
fitted_polynomial = np.poly1d(coeffs)
y_fitted = fitted_polynomial(x)

# Plot the fitted polynomial along with the data
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y_fitted, mode='lines', line=dict(color='green'), name='Fitted Cubic Polynomial'))
fig.add_trace(go.Scatter(x=x_data, y=y_data, mode='markers', marker=dict(size=5, color='red'), name='Noisy Data'))

fig.update_layout(title='Fitted Cubic Polynomial and Noisy Data', xaxis_title='x', yaxis_title='y')
fig.show()

Fitting Polynomials of Varying Degrees¶

Let’s fit polynomials from degree 1 to 10 and calculate the residual norm for each.

Even though the true underlying function is 4th order, adding higher-order terms will always reduce the residual error on the training data. This makes it non-trivial to determine the true order just by looking at the error! Let’s plot all of them:

fig = go.Figure()
fig.add_trace(go.Scatter(x=x_data, y=y_data, mode='markers', marker=dict(size=5, color='red'), name='Noisy Data'))

palette = colors.qualitative.Plotly
for degree in range(1, 11):
    coeffs = np.polyfit(x_data, y_data, degree)
    fitted_polynomial = np.poly1d(coeffs)
    y_fitted = fitted_polynomial(x)
    
    fig.add_trace(go.Scatter(x=x, y=y_fitted, mode='lines', 
                             line=dict(color=palette[degree % len(palette)]), 
                             name=f'Degree {degree}', 
                             visible='legendonly' if degree > 4 else True))
                             
    residuals = y_data - fitted_polynomial(x_data)
    residual_norm = np.linalg.norm(residuals)
    print(f"Degree {degree}: Residual Norm = {residual_norm:.2f}")

fig.update_layout(title='Fitted Polynomials of Varying Degrees (Toggle in Legend)', xaxis_title='x', yaxis_title='y')
fig.show()