Lab U7: Polynomial Approximation

Unit: Unit 7, Abstract vector spaces and polynomial approximation
Role: Core computation companion
Textbook sections: 7.3 Polynomial spaces and course connections; 7.4 Applications and computation recap
Core path: polynomial coefficient vectors, feature matrices, Taylor versus least squares, residual orthogonality, shrinking intervals, Gram-Schmidt checks, sampled planes, and Hessian quadratic forms
Main habit: read the code by identifying the vector space, the approximation subspace, the feature columns, the coefficient vector, the fitted values, and the residual check.

This lab is about interpretation. You are not expected to write long programs. Most questions ask you to predict, interpret, explain, or debug short code snippets.

Computational tools used in this lab

Before starting, review these parts of Appendix B, NumPy and SymPy Quick Reference for the Labs:

• Appendix B.2: NumPy arrays, vectors, matrices, and shapes
• Appendix B.4: Elementwise arithmetic versus linear algebra
• Appendix B.8: Constructing common arrays and design matrices
• Appendix B.9: Numerical checks and roundoff
• Appendix B.10: NumPy linear algebra commands

The main commands are np.array, np.linspace, np.column_stack, matrix products with @, transposes such as .T, np.linalg.lstsq, and np.linalg.norm.

Part 0. Setup: polynomial values and feature columns

Math reminder. The polynomial [ p(x)=1+2x+x^2+x^3 ] has coefficient vector [ \begin{bmatrix} 1\\2\\1\\1 \end{bmatrix} ] in the ordered basis (1,x,x^2,x^3).

Predict before running.

1. What should the first column of a polynomial feature matrix be?
2. What should the columns represent if the degree is (3)?
3. What does the vector y5 store?

import numpy as np

np.set_printoptions(precision=4, suppress=True)

def p(x):
    return 1 + 2*x + x**2 + x**3

def dp(x):
    return 2 + 2*x + 3*x**2

def feature_matrix_1d(x_values, degree):
    return np.column_stack([x_values**k for k in range(degree + 1)])

xs5 = np.array([-1.0, -0.5, 0.0, 0.5, 1.0])
y5 = p(xs5)

coeffs_p = np.array([1.0, 2.0, 1.0, 1.0])
Phi3 = feature_matrix_1d(xs5, 3)

Phi3.shape, Phi3, Phi3 @ coeffs_p, y5

Run and compare. The matrix Phi3 has one row for each sample point and four columns for (1,x,x^2,x^3). The product Phi3 @ coeffs_p gives the sampled values of (p), so it agrees with y5.

Interpretation check. A polynomial can be read two ways: as a function (p(x)), or as a coordinate vector after a basis has been chosen.

Common mistake. The vector coeffs_p stores coefficients. The vector y5 stores sampled values. They are not the same kind of vector.

Part 1. Affine feature matrix

Math reminder. The affine-polynomial subspace is [ \mathcal P_{\le1}=\operatorname{span}{1,x}. ] A sampled affine fit uses a design matrix with two columns: the sampled constant feature and the sampled (x)-feature.

Predict before running.

1. How many rows should A1 have?
2. How many columns should A1 have?
3. Which polynomial subspace do the columns represent?

A1 = feature_matrix_1d(xs5, 1)

A1.shape, A1

Run and compare. The matrix A1 has shape (5, 2). The five rows correspond to the five sample points. The two columns correspond to (1) and (x).

Interpretation check. A coefficient vector [ \mathbf c= \begin{bmatrix} c_0\\c_1 \end{bmatrix} ] produces fitted values for [ q(x)=c_0+c_1x ] by the product A1 @ c.

Part 2. Taylor line versus sampled least-squares line

Math reminder. At (0), the Taylor line for [ p(x)=1+2x+x^2+x^3 ] is [ T(x)=p(0)+p'(0)x=1+2x. ] The sampled least-squares line minimizes squared error at the sample points.

Predict before running.

1. Which coefficient vector should represent the Taylor line?
2. Which coefficient vector is computed by np.linalg.lstsq?
3. What should A1.T @ r_sampled check?

c_taylor = np.array([p(0.0), dp(0.0)])

c_sampled = np.linalg.lstsq(A1, y5, rcond=None)[0]
fit_sampled = A1 @ c_sampled
r_sampled = y5 - fit_sampled

c_taylor, c_sampled, A1.T @ r_sampled, np.linalg.norm(r_sampled)

Run and compare. The Taylor coefficient vector is approximately [ \begin{bmatrix} 1\\2 \end{bmatrix}. ] The sampled least-squares coefficient vector is approximately [ \begin{bmatrix} 1.5\\2.85 \end{bmatrix}

\begin{bmatrix} \frac32\\ \frac{57}{20} \end{bmatrix}. ]

Interpretation check. The expression A1.T @ r_sampled checks that the residual vector is orthogonal to the sampled constant feature and the sampled (x)-feature. It should be close to zero.

Common mistake. A least-squares residual usually is not the zero vector. Orthogonal to the feature columns does not mean equal to zero.

Part 3. Continuous (L^2) coefficients versus sampled coefficients

Math reminder. For the same polynomial, the continuous (L^2[-1,1])-least-squares line is [ q_{L^2}(x)=\frac43+\frac{13}{5}x. ] The sampled line from the five sample points is [ q_{\mathrm{sample}}(x)=\frac32+\frac{57}{20}x. ]

Predict before running.

1. Which coefficient vector is local?
2. Which coefficient vector comes from an integral over ([-1,1])?
3. Which coefficient vector comes from the five sample points?

c_continuous = np.array([4/3, 13/5], dtype=float)

grid = np.linspace(-1, 1, 1001)
A_grid = feature_matrix_1d(grid, 1)

def approx_squared_error(c):
    residual = p(grid) - A_grid @ c
    return np.trapz(residual**2, grid)

errors = np.array([
    approx_squared_error(c_taylor),
    approx_squared_error(c_continuous),
    approx_squared_error(c_sampled),
])

c_taylor, c_continuous, c_sampled, errors

Run and compare. The three coefficient vectors represent three different lines in the same subspace (\mathcal P_{\le1}). The approximate squared errors are computed over a dense grid on ([-1,1]), so the continuous (L^2) line should have the smallest value among the three for this grid-based check.

Interpretation check. The Taylor line is local. The continuous line depends on the interval. The sampled line depends on the sample points.

Review check. If the sample points changed, which coefficient vector would change?

Part 4. Debugging a residual check

Math reminder. For sampled least squares with design matrix (A), the residual condition is [ A^T\mathbf r=\mathbf0. ] The transpose matters.

Predict before running.

1. What is the shape of A1?
2. What is the shape of r_sampled?
3. Why is A1.T @ r_sampled the correct check?
4. What is wrong with trying to use A1 @ r_sampled?

correct_check = A1.T @ r_sampled
shape_information = {
    "A1 shape": A1.shape,
    "A1.T shape": A1.T.shape,
    "r_sampled shape": r_sampled.shape,
    "A1.T @ r_sampled shape": correct_check.shape,
}

correct_check, shape_information

Run and compare. The product A1.T @ r_sampled has one entry for each feature column. Here there are two features, (1) and (x), so the check has two entries.

Debug check. The product A1 @ r_sampled is not the residual-orthogonality check. It either has incompatible dimensions or represents the wrong idea. Least-squares residuals are tested against the columns of (A), so the transpose appears.

Part 5. Shrinking interval coefficients

Math reminder. The (L^2[-t,t])-projection of (x^2) onto (\operatorname{span}{1,x}) is [ P_t(x^2)=\frac{t^2}{3}. ] The (L^2[-t,t])-projection of (x^3) onto (\operatorname{span}{1,x}) is [ P_t(x^3)=\frac{3t^2}{5}x. ]

Predict before running.

1. What should happen to both coefficient arrays as (t\to0)?
2. How does this connect to the Taylor line at (0)?
3. Are these fixed-interval least-squares approximations exactly Taylor approximations?

ts = np.array([1.0, 0.5, 0.25, 0.125, 0.0625])

x2_const = ts**2 / 3
x3_slope = 3 * ts**2 / 5

np.column_stack([ts, x2_const, x3_slope])

Run and compare. The second column gives the constant term in (P_t(x^2)). The third column gives the slope in (P_t(x^3)). Both tend to zero as (t) shrinks.

Interpretation check. For fixed (t), these are least-squares approximations on an interval. In the limit (t\to0), they recover the Taylor linear approximation at (0).

Common mistake. The inner products (\int_{-t}^{t}p(x)q(x),dx) collapse as (t\to0). The projections can still have a useful limit because their coefficients are ratios of quantities that vanish at related rates.

Part 6. Quadratic feature matrix

Math reminder. Polynomial regression with quadratic features uses [ \operatorname{span}{1,x,x^2}. ] The model [ y\approx c_0+c_1x+c_2x^2 ] is nonlinear in (x), but linear in the coefficient vector.

Predict before running.

1. How many columns should A2 have?
2. What does each column represent?
3. What should A2.T @ r2 check?

A2 = feature_matrix_1d(xs5, 2)

c2 = np.linalg.lstsq(A2, y5, rcond=None)[0]
fit2 = A2 @ c2
r2 = y5 - fit2

A2.shape, c2, A2.T @ r2, np.linalg.norm(r2)

Run and compare. The matrix A2 has three columns, corresponding to (1,x,x^2). The vector c2 stores the coefficients of the sampled least-squares quadratic.

Interpretation check. The product A2.T @ r2 has three entries because the residual is checked against three feature columns: [ \sum_i r_i,\qquad \sum_i x_ir_i,\qquad \sum_i x_i^2r_i. ]

Review check. Why can a quadratic fit have a smaller residual norm than an affine fit?

Part 7. Gram-Schmidt check for polynomial features

Math reminder. With [ \langle f,g\rangle=\int_{-1}^{1}f(x)g(x),dx, ] Gram-Schmidt applied to (1,x,x^2,x^3) gives the orthogonal basis [ 1,\qquad x,\qquad x^2-\frac13,\qquad x^3-\frac35x. ]

The code below approximates the inner products using a dense grid.

Predict before running.

1. Which entries of the Gram matrix should be close to zero?
2. Why should the matrix be close to diagonal?
3. Does orthogonal mean normalized?

xg = np.linspace(-1, 1, 10001)

phi0 = np.ones_like(xg)
phi1 = xg
phi2 = xg**2 - 1/3
phi3 = xg**3 - (3/5)*xg

Phi_orth = np.column_stack([phi0, phi1, phi2, phi3])

def approx_ip(u, v):
    return np.trapz(u * v, xg)

G_orth = np.array([
    [approx_ip(Phi_orth[:, i], Phi_orth[:, j]) for j in range(4)]
    for i in range(4)
])

G_orth

Run and compare. The off-diagonal entries should be close to zero. Small nonzero values can appear because the integral is approximated numerically.

Interpretation check. In an orthogonal basis, the Gram matrix is diagonal. Projection computations decouple because each basis direction can be handled separately.

Common mistake. The diagonal entries are not all (1), so this basis is orthogonal but not orthonormal.

Part 8. Two-variable sampled best-fit plane

Math reminder. A sampled best-fit plane uses the approximation subspace [ \operatorname{span}{1,x,y}. ] For sample points ((x_i,y_i)), the design matrix has columns (1,x,y).

We use [ f(x,y)=1+2x-y+x^2+xy+2y^2+x^3. ]

Predict before running.

1. How many points are in a (5\times5) grid?
2. What should the three columns of A_plane represent?
3. What should A_plane.T @ r_plane check?

grid5 = np.linspace(-1, 1, 5)
X, Y = np.meshgrid(grid5, grid5)

xx = X.ravel()
yy = Y.ravel()

z = 1 + 2*xx - yy + xx**2 + xx*yy + 2*yy**2 + xx**3

A_plane = np.column_stack([np.ones_like(xx), xx, yy])
c_plane = np.linalg.lstsq(A_plane, z, rcond=None)[0]
r_plane = z - A_plane @ c_plane

A_plane.shape, c_plane, A_plane.T @ r_plane, np.linalg.norm(r_plane)

Run and compare. The grid has (25) sample points. The coefficient vector should be approximately [ \begin{bmatrix} 2.5\\2.85\\-1 \end{bmatrix}

\begin{bmatrix} \frac52\\\frac{57}{20}\\-1 \end{bmatrix}. ]

Interpretation check. The fitted plane is [ q_{\mathrm{sample}}(x,y)=\frac52+\frac{57}{20}x-y. ] The residual check says [ \sum_i r_i=0,\qquad \sum_i x_ir_i=0,\qquad \sum_i y_ir_i=0. ]

Review check. How is this the two-variable version of sampled line fitting?

Part 9. Hessian quadratic form

Math reminder. For the same function, [ H_f(\mathbf0)= \begin{bmatrix} 2&1\\ 1&4 \end{bmatrix}. ] The quadratic part of the second-order Taylor polynomial is [ \frac12\mathbf h^TH_f(\mathbf0)\mathbf h. ]

Predict before running.

1. What shape should H have?
2. What does each row of directions represent?
3. What does the expression 0.5 * h @ H @ h compute for one direction (\mathbf h)?

H = np.array([
    [2.0, 1.0],
    [1.0, 4.0],
])

directions = np.array([
    [1.0, 0.0],
    [0.0, 1.0],
    [1.0, -1.0],
    [2.0, 1.0],
])

quadratic_values = 0.5 * np.sum((directions @ H) * directions, axis=1)

directions, quadratic_values

Run and compare. Each row of directions is one vector (\mathbf h). The array quadratic_values stores [ \frac12\mathbf h^TH\mathbf h ] for each direction.

Interpretation check. These values evaluate the quadratic part of the local Taylor model in several directions.

Common mistake. The Hessian does not store the constant or linear term. It stores the quadratic part.

Part 10. Review code-reading bank

For each snippet, answer the interpretation question without changing the code.

Snippet A

A = np.column_stack([np.ones_like(xs), xs])
c = np.linalg.lstsq(A, y, rcond=None)[0]
r = y - A @ c
A.T @ r

1. What are the feature columns?
2. What does (c) store?
3. What does (A^T\mathbf r) check?

Snippet B

c_taylor = np.array([1.0, 2.0])
c_continuous = np.array([4/3, 13/5])
c_sampled = np.array([3/2, 57/20])

1. Which approximation is local?
2. Which approximation depends on the interval ([-1,1])?
3. Which approximation depends on sample points?

Snippet C

A = np.column_stack([np.ones_like(xx), xx, yy])
c = np.linalg.lstsq(A, z, rcond=None)[0]
r = z - A @ c
A.T @ r

1. What plane is being fitted?
2. What are the three residual equations?
3. Why are there three entries in the check?

Snippet D

H = np.array([[2.0, 1.0],
              [1.0, 4.0]])
h = np.array([1.0, -1.0])
0.5 * h @ H @ h

1. What mathematical expression is computed?
2. Which part of the Taylor polynomial does it evaluate?
3. Why does the factor (1/2) appear?

Lab summary

In this lab, the same linear algebra pattern appeared several times:

[ \text{features}\quad\longrightarrow\quad \text{design matrix}\quad\longrightarrow\quad \text{coefficients}\quad\longrightarrow\quad \text{fitted values}\quad\longrightarrow\quad \text{residual check}. ]

Taylor approximation is local. Continuous least squares depends on an interval or region. Sampled least squares depends on sample points. Gram-Schmidt changes coordinates to make projection easier. Hessians encode the quadratic part of a local polynomial model.