Lab U5: Hessian eigenvalues, least squares, and fixed nonlinear features

Unit: Unit 5, Eigenvalues and second-order optimization
Role: Required
Textbook sections: Hessians and the second derivative test; Least squares: three viewpoints; Least squares with fixed nonlinear features; Applications and computation recap
Unit 3 review: gradient descent updates, learning-rate diagnosis, debugging the descent sign, and a final-layer outer-product update
Unit 5 core path: Hessian eigenvalues and curvature; least squares through residual, gradient, and Hessian viewpoints; global minimality and uniqueness; fixed nonlinear features
Assessment note: This notebook is a supplement to the textbook, not a separate programming assignment. The main goal is to interpret computations, not to memorize code syntax.

No code submission is expected. Use the notebook for prediction, in-class discussion, and guided review.

Computational tools used in this lab

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

• Appendix A.4: Elementwise arithmetic versus linear algebra
• Appendix A.7: Reductions, axes, and row-wise computations
• Appendix A.8: Constructing common arrays and design matrices
• Appendix A.10: NumPy linear algebra commands
• Appendix A.12: Small Python patterns used in the labs

The goal is to interpret the mathematical computation, not to memorize every command.

Setup and reading habits

Before interpreting a computation:

1. identify the mathematical objects;
2. check their shapes;
3. identify the operation;
4. interpret the output.

The Unit 3 review uses gradient-descent updates and one outer product. The Unit 5 core uses Hessian eigenvalues, least-squares residuals and gradients, and fixed-feature design matrices.

Predict before running. Which objects below are vectors? Which object is a matrix? What shape should H_demo @ u have?

import numpy as np

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

x = np.array([4.0])
u = np.array([3.0, 4.0])

H_demo = np.array([
    [4.0, 0.0],
    [0.0, 1.0],
])

x.shape, u.shape, H_demo.shape, H_demo @ u

Run and compare. The objects x and u are vectors. The object H_demo is a matrix. The product H_demo @ u is a matrix-vector product, so its length equals the number of rows of H_demo.

Interpretation check. Shape checking should happen before mathematical interpretation.

Part A. Unit 3 review

The next four calculations review first-order optimization from Unit 3. They remain in this lab because Unit 5 will compare gradient-descent behavior with Hessian curvature.

A1. One-dimensional gradient descent

For (f(x)=x^2), [ f'(x)=2x. ] Gradient descent with learning rate (\alpha) uses [ x_{k+1}=x_k-\alpha f'(x_k). ]

Predict before running. Starting from (x_0=4), should the first few loss values (f(x_k)) decrease when (\alpha=0.2)?

def f1(x):
    return float(x @ x)

def grad_f1(x):
    return 2 * x

x = np.array([4.0])
alpha = 0.2
losses = []
iterates = []

for k in range(8):
    iterates.append(float(x[0]))
    losses.append(f1(x))
    x = x - alpha * grad_f1(x)

iterates, losses

Run and compare. The list iterates stores the successive values of (x_k). The list losses stores (f(x_k)).

Interpretation check. The line x = x - alpha * grad_f1(x) is the gradient descent update. The minus sign matters: the derivative points toward increasing (f), so the negative derivative points toward decreasing (f), at least locally.

Answer sketch. For this simple example, the loss decreases quickly for (\alpha=0.2). This is evidence of useful progress for this starting point and learning rate, but it is not a proof that every gradient descent run converges.

A2. Learning-rate comparison

Math reminder. The learning rate controls the step size. Too small can be slow. Too large can be unstable.

Predict before running. For the same function (f(x)=x^2), which learning rate should be slow, which should be useful, and which should be unstable: (0.05), (0.20), or (1.05)?

def run_gd(alpha, steps=4):
    x = np.array([4.0])
    losses = []
    for k in range(steps):
        losses.append(round(f1(x), 3))
        x = x - alpha * grad_f1(x)
    return losses

{alpha: run_gd(alpha) for alpha in [0.05, 0.20, 1.05]}

Run and compare. The learning rate (\alpha=0.05) makes slow but steady progress. The learning rate (\alpha=0.20) makes faster useful progress. The learning rate (\alpha=1.05) is unstable in this example.

Interpretation check. A loss table is a diagnostic. It tells what happened for the displayed iterates; it does not by itself prove that a global minimum has been found.

Common mistake. A larger learning rate is not always better.

A3. Debugging the sign

Math reminder. Gradient descent uses a minus sign: [ x_{k+1}=x_k-\alpha\nabla f(x_k). ] Changing the minus sign to a plus sign gives a gradient ascent step.

Predict before running. Which loss list should decrease: the one labeled "minus sign" or the one labeled "plus sign"?

def run_gd_with_sign(alpha, sign=-1, steps=5):
    x = np.array([4.0])
    losses = []
    for k in range(steps):
        losses.append(round(f1(x), 3))
        x = x + sign * alpha * grad_f1(x)
    return losses

{
    "minus sign": run_gd_with_sign(0.2, sign=-1),
    "plus sign": run_gd_with_sign(0.2, sign=1),
}

Run and compare. The "minus sign" version decreases the loss. The "plus sign" version increases the loss.

Interpretation check. The expression grad_f1(x) points in the direction of steepest increase for this one-dimensional function. Subtracting it moves in the opposite direction.

Debugging habit. When a gradient descent loop makes the loss increase immediately, first check the sign and the learning rate.

A4. Final-layer outer-product update

Suppose a final linear layer has the form [ \boldsymbol{\ell}=W\mathbf h, ] and let [ \mathbf g=\nabla_{\boldsymbol{\ell}}\mathcal L ] be the gradient of the loss with respect to the output scores. By the Unit 3 chain-rule calculation, [ \nabla_W\mathcal L=\mathbf g\mathbf h^T. ]

If (\mathbf g\in\mathbb R^m) and (\mathbf h\in\mathbb R^d), then this gradient is an (m\times d) matrix. In NumPy, np.outer(g, h) forms the outer product when the vectors are stored as one-dimensional arrays.

Predict before running. What shape should G have? Why should its rank be at most one?

g = np.array([1.0, -2.0, 3.0])
h = np.array([4.0, 0.0])

G = np.outer(g, h)
W = np.ones((3, 2))
alpha = 0.1
W_new = W - alpha * G

G, G.shape, np.linalg.matrix_rank(G), W_new

Run and compare. The matrix G has shape (3, 2). It is the outer product [ \mathbf g\mathbf h^T. ] Every column is a scalar multiple of (\mathbf g), so its rank is at most one.

Interpretation check. The line

W_new = W - alpha * G

is a gradient-descent update of the final-layer matrix for this one example.

Common mistakes.

• For one-dimensional NumPy arrays, g @ h.T does not form the outer product.
• The update matrix has rank at most one. This does not imply that W or W_new has rank at most one.

Part B. Unit 5 core

The remaining calculations use Unit 5 tools: Hessian eigenvalues, quadratic curvature, least-squares structure, and fixed nonlinear features.

B1. A quadratic loss: gradient and Hessian

Consider [ f(\mathbf x) = \frac12\mathbf x^T H\mathbf x+\mathbf b^T\mathbf x, ] where (H) is symmetric. Then [ \nabla f(\mathbf x)=H\mathbf x+\mathbf b, \qquad H_f(\mathbf x)=H. ]

The gradient determines a first-order update. The Hessian is the matrix whose eigenvalues describe second-order curvature.

Predict before running.

1. What shape should grad_quad(x) have?
2. How can the critical point be found from (H\mathbf x+\mathbf b=\mathbf0)?
3. What should positive Hessian eigenvalues say about that critical point?

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

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

def f_quad(x):
    return float(0.5 * x @ H @ x + b @ x)

def grad_quad(x):
    return H @ x + b

x = np.array([2.0, -1.0])
alpha = 0.1

g = grad_quad(x)
x_next = x - alpha * g

x_star = np.linalg.solve(H, -b)
hessian_eigenvalues = np.linalg.eigvalsh(H)

(
    f_quad(x),
    g,
    x_next,
    f_quad(x_next),
    x_star,
    grad_quad(x_star),
    hessian_eigenvalues,
)

Run and compare.

• The gradient g has the same shape as x, so the update is defined.
• The displayed step lowers the loss from f_quad(x) to f_quad(x_next).
• The vector x_star solves [ H\mathbf x+\mathbf b=\mathbf0, ] so its gradient is numerically zero.
• Both Hessian eigenvalues are positive.

Interpretation check. The gradient supplies the update direction. The Hessian eigenvalues classify the critical point as a strict local minimum.

Common mistake. The expression H @ x + b is the gradient. The Hessian is H.

B2. Hessian eigenvalues and curvature

At a critical point, [ f(\mathbf a+\mathbf h)-f(\mathbf a) \approx \frac12\mathbf h^T H_f(\mathbf a)\mathbf h. ]

For a symmetric Hessian:

• all positive eigenvalues indicate a local minimum;
• mixed positive and negative eigenvalues indicate a saddle point;
• a zero eigenvalue makes the second derivative test inconclusive.

Predict before running. Which matrix below represents each of these three cases?

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

H_saddle = np.array([
    [4.0, 0.0],
    [0.0, -1.0],
])

H_flat = np.array([
    [4.0, 0.0],
    [0.0, 0.0],
])

(
    np.linalg.eigvalsh(H_min),
    np.linalg.eigvalsh(H_saddle),
    np.linalg.eigvalsh(H_flat),
)

Run and compare.

• H_min has two positive eigenvalues, so it is positive definite.
• H_saddle has one positive and one negative eigenvalue, so it is indefinite.
• H_flat has one positive and one zero eigenvalue.

Interpretation check. At a critical point, the first matrix gives a local minimum and the second gives a saddle point. The third is inconclusive because the quadratic approximation is flat in one direction.

Common mistake. A zero eigenvalue is not a classification.

B3. Curvature and stable step sizes

This calculation connects the Unit 3 update rule with a Unit 5 curvature diagnosis. For [ f(\mathbf x)=\frac12\mathbf x^T H\mathbf x, ] the gradient is (H\mathbf x). Large Hessian eigenvalues correspond to steep curvature directions, so a learning rate can be too large because of one direction.

Predict before running. Which learning rate should be slow, which should be useful, and which should be unstable: (0.05), (0.40), or (0.60)?

H = np.diag([4.0, 1.0])

def f_curv(x):
    return float(0.5 * (x @ H @ x))

def grad_curv(x):
    return H @ x

def run_quadratic_gd(alpha, steps=6):
    x = np.array([4.0, 4.0])
    losses = []
    for k in range(steps):
        losses.append(round(f_curv(x), 3))
        x = x - alpha * grad_curv(x)
    return losses

{alpha: run_quadratic_gd(alpha) for alpha in [0.05, 0.40, 0.60]}

Run and compare. The learning rate (\alpha=0.05) is slow. The learning rate (\alpha=0.40) makes useful progress. The learning rate (\alpha=0.60) is unstable for this quadratic loss.

Interpretation check. The largest eigenvalue of the Hessian corresponds to the steepest curvature direction. A step size can be too large because of one steep direction even when another direction is less steep.

B4. Least squares: three viewpoints

For [ L(\mathbf x)=|A\mathbf x-\mathbf b|^2, ] the same least-squares solution has three readings:

1. Projection [ A^T(\mathbf b-A\widehat{\mathbf x})=\mathbf0. ]

2. Gradient [ \nabla L(\widehat{\mathbf x}) = 2A^T(A\widehat{\mathbf x}-\mathbf b) = \mathbf0. ]

3. Curvature [ H_L=2A^TA. ]

If (\widehat{\mathbf x}) satisfies the normal equations, then [ L(\widehat{\mathbf x}+\mathbf z) = L(\widehat{\mathbf x})+|A\mathbf z|^2. ] This proves global minimality. Null-space directions determine whether the coefficient vector is unique.

Predict before running.

• Should A.T @ r and grad_at_xhat be close to zero?
• Should loss_gap equal direction_cost?
• What should happen after adding a null-space vector to a least-squares coefficient vector?

def squared_loss(A, b, x):
    return float(np.linalg.norm(A @ x - b)**2)

# Full-column-rank example
A = np.array([
    [1.0, 0.0],
    [1.0, 1.0],
    [1.0, 2.0],
])

b_ls = np.array([1.0, 2.0, 2.0])

xhat = np.linalg.lstsq(A, b_ls, rcond=None)[0]
r = b_ls - A @ xhat

grad_at_xhat = 2 * A.T @ (A @ xhat - b_ls)
H_lstsq = 2 * A.T @ A

z = np.array([1.5, -0.5])
loss_gap = squared_loss(A, b_ls, xhat + z) - squared_loss(A, b_ls, xhat)
direction_cost = float(np.linalg.norm(A @ z)**2)

full_rank_checks = {
    "xhat": xhat,
    "residual": r,
    "A.T @ r": A.T @ r,
    "gradient": grad_at_xhat,
    "Hessian eigenvalues": np.linalg.eigvalsh(H_lstsq),
    "loss gap": loss_gap,
    "||A z||^2": direction_cost,
    "global identity": np.isclose(loss_gap, direction_cost),
}

# Dependent-column example
A_dep = np.array([
    [1.0, 2.0],
    [2.0, 4.0],
    [3.0, 6.0],
])

b_dep = np.array([1.0, 2.0, 2.5])

xhat_dep = np.linalg.lstsq(A_dep, b_dep, rcond=None)[0]
z_null = np.array([-2.0, 1.0])

pred_dep = A_dep @ xhat_dep
pred_shifted = A_dep @ (xhat_dep + z_null)

dependent_checks = {
    "xhat": xhat_dep,
    "null direction": A_dep @ z_null,
    "same predictions": np.allclose(pred_dep, pred_shifted),
    "same loss": np.isclose(
        squared_loss(A_dep, b_dep, xhat_dep),
        squared_loss(A_dep, b_dep, xhat_dep + z_null),
    ),
    "Hessian eigenvalues": np.linalg.eigvalsh(2 * A_dep.T @ A_dep),
}

full_rank_checks, dependent_checks

Run and compare.

For the full-column-rank example:

• A.T @ r is numerically zero: the residual is orthogonal to the columns of (A).
• grad_at_xhat is numerically zero: the same vector is a critical point of the squared-residual loss.
• the Hessian eigenvalues are positive;
• loss_gap equals direction_cost, verifying [ L(\widehat{\mathbf x}+\mathbf z)-L(\widehat{\mathbf x}) = |A\mathbf z|^2. ]

For the dependent-column example:

• z_null is a nonzero null-space direction;
• adding it does not change the prediction vector;
• adding it does not change the loss;
• the Hessian has a zero eigenvalue.

Interpretation check. The fitted vector is still a global closest vector, but the coefficient vector is not unique when (A) has a nontrivial null space.

Common mistakes.

• A.T @ r being zero does not mean that r is zero.
• A local Hessian test alone does not prove the global least-squares result; the exact loss identity does.
• Existence of a least-squares fit does not require independent columns. Independence gives uniqueness of the coefficient vector.

B5. Fixed nonlinear features

Let [ \boldsymbol{\phi}(t) = \begin{bmatrix} 1\\ \tanh(t)\\ \tanh(t-1)\\ \tanh(t+1) \end{bmatrix}, \qquad N_{\mathbf c}(t) = \boldsymbol{\phi}(t)^T\mathbf c. ]

After evaluating these fixed features at data inputs (t_1,\ldots,t_m), their values form a design matrix (A). The prediction vector is [ A\mathbf c, ] and the squared-error loss is [ L(\mathbf c)=|A\mathbf c-\mathbf y|^2. ]

The model is nonlinear in (t), but it is linear in the trained coefficient vector (\mathbf c).

Predict before running.

1. Which objects are fixed?
2. Which vector is fitted?
3. Should A.T @ residual be close to zero?

t = np.array([-2.0, -1.0, 0.0, 1.0, 2.0])
y = np.array([4.0, 1.0, 0.0, 1.0, 4.0])

A = np.column_stack([
    np.ones_like(t),
    np.tanh(t),
    np.tanh(t - 1),
    np.tanh(t + 1),
])

c = np.linalg.lstsq(A, y, rcond=None)[0]
yhat = A @ c
residual = y - yhat

loss = float(np.linalg.norm(residual)**2)
normal_check = A.T @ residual

fit_table = np.column_stack([
    t,
    y,
    yhat,
    residual,
])

A.shape, c, fit_table, loss, normal_check

Run and compare.

• The four columns of A are the fixed feature vectors.
• The vector c is fitted by least squares.
• The vector yhat contains the fitted values.
• The variable loss stores [ |A\mathbf c-\mathbf y|^2. ]
• The vector normal_check is numerically zero.

The columns of fit_table are:

1. input (t_i);
2. target (y_i);
3. fitted value (N_{\mathbf c}(t_i));
4. residual (y_i-N_{\mathbf c}(t_i)).

Interpretation check. The model is nonlinear in the input (t), but linear in (\mathbf c). That is why fitting (\mathbf c) with squared error is an ordinary least-squares problem.

Warning. If shifts or other parameters inside the feature functions are also trained, then the entries of (A) depend on trainable parameters. One fixed least-squares solve no longer trains the whole model.

Part C. Review bank

These are short discussion or self-review questions. They do not require writing a new program.

Unit 3 review

1. Update rule. What mathematical update is represented by
   x = x - alpha * grad_f(x)?

2. Learning rate. What role does alpha play?

3. Shape check. Why must grad_f(x) have the same shape as x?

4. Sign check. What usually happens if the minus sign is changed to a plus sign?

5. Loss history. What does losses.append(f(x)) record? Does a decreasing finite list prove that a global minimum has been found?

6. Final-layer gradient. If [ \nabla_W\mathcal L=\mathbf g\mathbf h^T, ] what is its shape when (\mathbf g\in\mathbb R^m) and (\mathbf h\in\mathbb R^d)?

7. Rank-one warning. Why does (\mathbf g\mathbf h^T) have rank at most one, and what does this not imply about (W_{\mathrm{new}})?

Unit 5 core

8. Gradient versus Hessian. For [ f(\mathbf x)=\frac12\mathbf x^TH\mathbf x+\mathbf b^T\mathbf x, ] what are the gradient and Hessian?

9. Positive eigenvalues. What do all-positive Hessian eigenvalues imply at a critical point?

10. Mixed eigenvalues. What do positive and negative Hessian eigenvalues imply at a critical point?

11. Zero eigenvalue. Why does a zero Hessian eigenvalue make the test inconclusive?

12. Curvature and step size. Why can one large Hessian eigenvalue make a learning rate unstable?

13. Least-squares solve. What does
    np.linalg.lstsq(A, b, rcond=None)[0] return?

14. Residual. In r = b - A @ xhat, what is r?

15. Two readings. How do projection and gradient interpret A.T @ r?

16. Global identity. What does [ L(\widehat{\mathbf x}+\mathbf z) = L(\widehat{\mathbf x})+|A\mathbf z|^2 ] prove?

17. Uniqueness. What condition on (A) makes the least-squares coefficient vector unique?

18. Null direction. What happens to predictions and loss if [ \mathbf z\in\operatorname{null}(A)? ]

19. Fixed features. What is stored in the columns of the fixed-feature design matrix?

20. Trained vector. Which vector is fitted in [ \mathbf y_{\mathrm{hat}}=A\mathbf c? ]

21. Squared-error loss. What loss is minimized in the fixed-feature problem, and is the model linear in (t) or in (\mathbf c)?

22. Trainable features. Why does one fixed least-squares solve no longer suffice when parameters inside the feature functions are trained?

Answer sketches

1. (\mathbf x_{k+1}=\mathbf x_k-\alpha\nabla f(\mathbf x_k)).
2. It is the learning rate or step size.
3. The update subtracts one vector from another.
4. The step points toward steepest increase rather than steepest decrease.
5. It records the current loss. No finite decreasing list proves a global minimum or convergence.
6. The gradient has shape (m\times d).
7. Every column is a scalar multiple of (\mathbf g). This does not imply that (W) or (W_{\mathrm{new}}) has rank one.
8. (\nabla f(\mathbf x)=H\mathbf x+\mathbf b), and the Hessian is (H).
9. The Hessian is positive definite, so the critical point is a strict local minimum.
10. The Hessian is indefinite, so the critical point is a saddle point.
11. The quadratic approximation is flat in one direction, so higher-order terms may matter.
12. The corresponding direction has steep curvature, so a large step can overshoot there.
13. A least-squares coefficient vector.
14. The residual (\mathbf b-A\widehat{\mathbf x}).
15. Projection reads it as residual orthogonality; the gradient calculation reads it as a zero-gradient condition.
16. It proves that (\widehat{\mathbf x}) is a global minimizer.
17. The columns of (A) must be linearly independent, equivalently (\operatorname{null}(A)={\mathbf0}).
18. The prediction vector and the loss do not change.
19. Values of the fixed feature functions at the data inputs.
20. The coefficient vector (\mathbf c).
21. The loss is (|A\mathbf c-\mathbf y|^2). The model may be nonlinear in (t), but it is linear in (\mathbf c).
22. The design matrix then depends on trainable parameters, so the loss is not one fixed quadratic function of all trained parameters.