Optimization in several variables

Lesson 3 introduced critical points, the second derivative test and Lagrange multipliers. Optimization is so central to engineering that it deserves a closer look. Designers minimize weight, cost, energy use and error, and maximize strength, efficiency and output, almost always with several variables and several constraints. This lesson strengthens the theory and widens the toolkit. It explains the second derivative test through the second-order Taylor expansion and the Hessian matrix, which extends naturally to three or more variables; finds absolute extremes on closed regions; derives the method of least squares as an optimization problem; looks more deeply at Lagrange multipliers, including what the multiplier means and how to handle two constraints and inequality constraints; and finishes with the numerical methods, gradient descent and Newton's method, that solve the large problems found in practice.

The Hessian and the second derivative test

Near a point \(\mathbf{a}\), a smooth function of several variables is approximated by its second-order Taylor expansion. Writing \(\mathbf{h}\) for the small displacement from \(\mathbf{a}\):

Here H is the Hessian matrix of second partial derivatives, which is symmetric because mixed partial derivatives are equal:

At a critical point the gradient term vanishes, so the behavior is set by the quadratic form \(\mathbf{h}^TH\mathbf{h}\) (MATH-125). If it is positive for every direction, the point is a local minimum; if negative for every direction, a local maximum; if positive in some directions and negative in others, a saddle. Since a symmetric matrix can be diagonalized by rotating the axes, the test reduces to the signs of the eigenvalues of H:

• All eigenvalues positive (H positive definite): local minimum.
• All eigenvalues negative: local maximum.
• Eigenvalues of both signs: saddle point.
• Some eigenvalue zero: the test is inconclusive.

In two variables, the determinant \(D = f_{xx}f_{yy} - f_{xy}^2\) is the product of the two eigenvalues, which is why the test in Lesson 3 works: D > 0 means the eigenvalues share a sign, and \(f_{xx}\) tells which. The eigenvectors give the principal directions of curvature, the directions in which the surface curves most and least.

Absolute extremes on a closed region

A continuous function on a closed, bounded region always attains an absolute maximum and minimum. To find them: locate the critical points inside the region, find the extremes on the boundary (by parametrizing it, or with Lagrange multipliers), and compare all the candidate values.

Least squares as optimization

When a model \(y = mx + c\) is fitted to measured points \((x_i, y_i)\), the best fit in the least-squares sense minimizes the sum of squared errors:

Setting \(\partial E/\partial m = 0\) and \(\partial E/\partial c = 0\) gives the normal equations:

These are the same normal equations that MATH-125 derives geometrically by projection. The Hessian of E has positive eigenvalues whenever the \(x_i\) are not all equal, so the critical point is the minimum.

Lagrange multipliers revisited

What the multiplier means

In the method \(\nabla f = \lambda\nabla g\) with constraint \(g = c\), the multiplier is not just an auxiliary unknown. If \(f^*(c)\) is the optimal value, then

So \(\lambda\) is the rate at which the best achievable objective changes as the constraint is relaxed. Economists call it a shadow price; for an engineer it answers questions such as "how much more material does each extra liter of capacity cost?" without re-solving the problem.

Two constraints

With two constraints \(g = c_1\) and \(h = c_2\), the feasible set is usually a curve, and at the optimum the gradient of f lies in the plane spanned by the two constraint gradients:

Inequality constraints

Real constraints are usually inequalities: stress must not exceed the allowable value, a budget must not be overspent. At the optimum, each inequality is either inactive (not binding, so it can be ignored locally) or active (holding with equality, so it acts like an equality constraint with a multiplier of the correct sign). These are the Karush-Kuhn-Tucker (KKT) conditions, the basis of the optimization software used for structural design and operations planning. In practice, one tries the cases: solve with no constraints active, check feasibility, then make constraints active one at a time.

Numerical optimization

Most engineering optimization problems have too many variables, or functions too complex, for these hand methods. Numerical methods search iteratively. Gradient descent steps repeatedly in the direction of steepest descent, \(-\nabla f\) (Lesson 3), with a step size \(\alpha\):

It is simple and needs only the gradient, but it is slow when the bowl is long and narrow, that is, when the eigenvalues of the Hessian differ greatly: the steps zigzag across the valley. Newton's method uses the Hessian as well, jumping to the minimum of the local quadratic model, \(\mathbf{x}_{k+1} = \mathbf{x}_k - H^{-1}\nabla f\), and converges in very few steps near the optimum, at the cost of computing and solving with H. Practical software uses compromises (quasi-Newton methods such as BFGS) that build up an approximation of H from successive gradients.

import numpy as np

# Narrow valley: f = (x - 1)^2 + 10 (y - 2)^2, minimum at (1, 2)
def f(p):
    return (p[0] - 1)**2 + 10 * (p[1] - 2)**2

def grad(p):
    return np.array([2 * (p[0] - 1), 20 * (p[1] - 2)])

H = np.array([[2.0, 0.0], [0.0, 20.0]])
start = np.array([-2.0, 0.0])

for alpha in [0.01, 0.05, 0.09, 0.11]:
    p = start.copy()
    for k in range(1, 5001):
        p = p - alpha * grad(p)
        if f(p) < 1e-10 or f(p) > 1e12:
            break
    result = "converged" if f(p) < 1e-10 else "diverged"
    print(f"Gradient descent, step {alpha:.2f}: {result} after {k} iterations")

p = start - np.linalg.solve(H, grad(start))
print(f"Newton: one step to ({p[0]:.3f}, {p[1]:.3f}), f = {f(p):.1e}")

Gradient descent, step 0.01: converged after 625 iterations
Gradient descent, step 0.05: converged after 120 iterations
Gradient descent, step 0.09: converged after 64 iterations
Gradient descent, step 0.11: diverged after 66 iterations
Newton: one step to (1.000, 2.000), f = 0.0e+00

The output shows the trade-off clearly. The largest stable step for gradient descent is set by the steepest direction (here \(\alpha < 2/20 = 0.1\)), and progress along the shallow direction is then slow. Rescaling the variables so that the bowl is round, called preconditioning, is one of the most effective tricks in numerical optimization, and it is the same idea as choosing good units and dimensionless variables.