Optimization with Linear Algebra

Gradient Descent

In the previous lesson we established that the gradient ∇f(x) points in the direction of steepest ascent of f at x, and that −∇f(x) therefore points in the direction of steepest descent. Gradient descent is the iterative algorithm that exploits this fact: starting from an initial guess x₀, it repeatedly takes steps in the downhill direction:

The algorithm is deceptively simple — one line of math — yet it underlies virtually every trained machine learning model in existence. The scalar α > 0 is the step size (also called the learning rate), and it controls how far we move in the downhill direction at each iteration. Choosing α wisely is crucial: too large and the iterates may overshoot and diverge; too small and convergence is painfully slow.

The convergence behavior of gradient descent is governed by the eigenvalues of the Hessian. For a strictly convex quadratic f(x) = ½xᵀAx − bᵀx with A positive definite, gradient descent converges geometrically: the error at step k satisfies ‖x_k − x*‖ ≤ ((λ_max − λ_min)/(λ_max + λ_min))^k ‖x_0 − x*‖. The ratio κ = λ_max/λ_min is the condition number of A — a large condition number means a poorly conditioned problem and slow convergence, no matter how well you tune α. This is a fundamental reason why preconditioning (transforming the problem to reduce κ) matters in practice.

For non-quadratic f, the step size can be adapted at each iteration using a line search: rather than fixing α, find the α_k that minimizes f(x_k − α∇f(x_k)) along the gradient direction. In deep learning, exact line searches are too expensive; instead, practitioners use fixed learning rates with schedules, momentum (adding a fraction of the previous step to the current), or adaptive methods like Adam that estimate per-parameter step sizes.

The Loss Landscape

To understand why optimization is hard in general, it helps to visualize the loss landscape — the graph of the objective function f over its domain. For a convex function, the landscape has a single bowl-shaped valley: any local minimum is automatically the global minimum, and gradient descent is guaranteed to find it. For a non-convex function (like a neural network loss), the landscape is rugged, riddled with local minima, saddle points, and flat plateaus.

Three types of critical points (where ∇f = 0) are possible:

• Local minimum: Hessian is positive definite (all eigenvalues positive). f increases in every direction from this point.
• Local maximum: Hessian is negative definite (all eigenvalues negative). f decreases in every direction — rare as optimization targets.
• Saddle point: Hessian is indefinite (mixed-sign eigenvalues). f increases in some directions and decreases in others. These are exponentially more common than local minima in high-dimensional problems and are a primary obstacle for first-order methods.

In high-dimensional optimization (parameters numbering in millions), the landscape geometry changes character. Research suggests that for over-parameterized models, most local minima have similar function values to the global minimum — a reassuring result. However, saddle points can cause gradient descent to stall: at a saddle point, the gradient is zero, so gradient descent halts. The slightest perturbation (noise from stochastic batches) usually escapes saddle points in practice, which is one reason why stochastic gradient descent (SGD) — using noisy gradient estimates from random minibatches — often outperforms full-batch gradient descent on non-convex problems.

Newton's Method

Gradient descent uses only first-order information (the gradient). Newton's method uses second-order information — the Hessian — to make much larger, more targeted steps. The idea: at the current iterate x_k, form the second-order Taylor approximation of f:

q(δ) = f(x_k) + ∇f(x_k)ᵀδ + ½δᵀH(x_k)δ

This quadratic q(δ) is minimized exactly by setting its gradient to zero: ∇_δ q = ∇f(x_k) + H(x_k)δ = 0, giving δ = −H(x_k)⁻¹∇f(x_k). The Newton update moves to x_k + δ:

Newton's method has quadratic convergence near a strict local minimum: once the iterates are close to x*, the number of correct decimal places roughly doubles with each step. Gradient descent has only linear convergence (constant fraction of error removed per step). In practice, however, Newton's method has serious drawbacks:

• Cost: Computing H(x_k) requires O(n²) storage and O(n²) gradient evaluations. Solving the linear system H δ = −∇f costs O(n³) (Cholesky or LU). For n = 10⁶ (typical neural network), this is completely infeasible.
• Non-convex case: If H is indefinite (not positive definite), the Newton direction may point uphill rather than downhill. Safeguards like Hessian modification (adding λI to force positive definiteness) are required.
• Quasi-Newton methods (L-BFGS, BFGS) approximate H⁻¹ using gradient differences between iterates, achieving superlinear convergence without ever forming the full Hessian. They are the method of choice for moderate-dimensional smooth optimization.

Despite these practical limitations, Newton's method is the gold standard conceptually. Every advanced optimization algorithm can be understood as an approximation or modification of Newton's method.

Convexity and Positive Definite Hessians

A function f : ℝⁿ → ℝ is convex if its graph lies below or on any chord connecting two points on the graph. Formally, for all x, y in the domain and all λ ∈ [0, 1]:

For twice-differentiable functions, convexity has an elegant characterization in terms of the Hessian: f is convex if and only if H(x) is positive semidefinite (PSD) for all x in the domain. If H(x) is positive definite (PD) for all x, then f is strictly convex, which guarantees a unique global minimum.

Why does this matter so much for optimization? For convex f, every local minimum is a global minimum. Gradient descent (with appropriate step size) converges to the global minimum. There are no saddle points to get stuck in (all critical points are global minima). Duality theory (Lagrange multipliers, KKT conditions) gives tight bounds and certificates of optimality. In summary: convexity transforms optimization from a search problem into a computation.

• H(x) ≻ 0 (PD) everywhere: f is strictly convex → unique global minimum exists → gradient descent and Newton converge globally.
• H(x) ⪰ 0 (PSD) everywhere: f is convex → global minimum exists but may not be unique → still tractable.
• H(x) indefinite somewhere: f is non-convex → multiple local minima possible → hard in general.

Important convex functions in engineering: any norm ‖·‖, squared norms ‖Ax − b‖², log-sum-exp (used in logistic regression), the negative log-likelihood of exponential family distributions. The sum of convex functions is convex; the composition f(Ax + b) is convex when f is convex; the maximum of convex functions is convex. These closure properties let you recognize convexity in complex objectives.

Quadratic Forms and Their Optimization

The most important class of optimization problems in linear algebra is the quadratic form. A general quadratic function of x ∈ ℝⁿ takes the form f(x) = xᵀAx + bᵀx + c, where A ∈ ℝⁿˣⁿ is a symmetric matrix, b ∈ ℝⁿ is the linear coefficient vector, and c ∈ ℝ is a constant. Quadratics arise everywhere: least-squares regression (A = XᵀX, b = −2Xᵀy), energy minimization in physics, regularized optimization, Gaussian likelihood functions.

Taking the gradient: ∇f(x) = 2Ax + b. Setting this to zero gives the optimality condition 2Ax + b = 0, or equivalently the linear system Ax = −b/2. When A is invertible, there is a unique critical point:

The connection to linear systems is deep: solving a positive definite linear system Ax = b and minimizing the quadratic f(x) = ½xᵀAx − bᵀx are equivalent problems. Iterative methods like the conjugate gradient method exploit this duality — they are simultaneously linear algebra algorithms (solving Ax = b) and optimization algorithms (minimizing f). The CG method is optimal in a certain sense: it minimizes f over Krylov subspaces, achieving convergence in at most n steps for n × n A.

The Ridge Regression Example

Ridge regression (ℓ₂ regularization) illustrates quadratic optimization perfectly. The objective is f(w) = ‖Xw − y‖² + λ‖w‖² = wᵀ(XᵀX + λI)w − 2yᵀXw + yᵀy. Here A = XᵀX + λI and b = −2Xᵀy. The regularization term λ‖w‖² adds λ to every eigenvalue of XᵀX, ensuring A is positive definite even when X is rank-deficient. The closed-form solution is w* = (XᵀX + λI)⁻¹Xᵀy — the regularized normal equations. As λ → 0 we recover ordinary least squares; as λ → ∞ we get w* → 0 (over-regularized).

Lagrange Multipliers — The Linear Algebra View

Many practical optimization problems have constraints. We want to minimize f(x) subject to one or more equality constraints g(x) = 0. For example: minimize power consumption subject to a data-rate constraint; minimize weight subject to a structural load constraint; minimize reconstruction error subject to a unit-norm constraint (as in PCA).

The method of Lagrange multipliers turns a constrained problem into an unconstrained one by introducing auxiliary variables λ (the multipliers) that penalize constraint violations. The key insight is geometric: at a constrained optimum, the gradient of f must be a linear combination of the gradients of the constraint functions — otherwise, there would be a feasible direction that decreases f while staying on the constraint surface. This geometric condition, ∇f = λᵀ∇g, is precisely the first-order optimality condition.

When both f and g are differentiable, the KKT conditions give a system of equations:

1. Stationarity: ∇_x f(x*) − λᵀ∇_x g(x*) = 0 (gradient of Lagrangian vanishes)
2. Primal feasibility: g(x*) = 0 (constraints are satisfied)

For a single equality constraint g(x) = 0 with m constraints, stacking these conditions gives a block linear system. With g linear (g(x) = Cx − d), the KKT system becomes:

[2A, Cᵀ; C, 0] [x; λ] = [−b; d]

This saddle-point linear system is the fundamental object in constrained quadratic programming. Its coefficient matrix is symmetric but indefinite (it has both positive and negative eigenvalues), which requires specialized solvers — not Cholesky. The block structure is exploited in interior-point methods, the dominant algorithms for large-scale convex optimization.

Eigenvalue Problems as Constrained Optimization

There is a beautiful connection to eigenvalue problems. Consider: maximize xᵀAx subject to ‖x‖ = 1 (a quadratic objective on the unit sphere). The Lagrangian is L(x, λ) = xᵀAx − λ(xᵀx − 1). Setting ∇_x L = 2Ax − 2λx = 0 gives Ax = λx — the eigenvalue equation! The maximum of xᵀAx on the unit sphere equals λ_max (the largest eigenvalue of A), achieved at the corresponding eigenvector. This proves the Rayleigh quotient theorem and shows that PCA (finding the direction of maximum variance) is a constrained optimization problem whose solution is the top eigenvector of the covariance matrix.

Applications

Machine Learning

Every ML model trained by gradient descent is a direct application of this lesson. In logistic regression, the loss L(w) = −Σ yᵢ log σ(wᵀxᵢ) − (1−yᵢ) log(1−σ(wᵀxᵢ)) is convex in w (the Hessian XᵀDX with D diagonal positive is PSD), so gradient descent finds the global optimum. In deep networks the loss is non-convex, but stochastic gradient descent (SGD) and its variants — momentum SGD, Adam, AdaGrad — have proven highly effective in practice. Adam maintains running estimates of the first and second moments of the gradient, effectively constructing a diagonal approximation to H⁻¹ that adapts per-parameter step sizes. The update rule is: θ ← θ − η m̂/(√v̂ + ε), where m̂ and v̂ are bias-corrected moment estimates — a sophisticated variant of Newton's method diagonal in disguise.

Signal Processing and the Wiener Filter

The Wiener filter is the optimal linear filter for signal estimation in the minimum mean-squared error (MMSE) sense. Given a desired signal d and a noisy observed signal x, the optimal filter w* minimizes E[‖d − wᵀx‖²] over all weight vectors w. Setting ∂/∂w = 0 gives the Wiener-Hopf equations: R_xx w = r_xd, where R_xx = E[xxᵀ] is the input autocorrelation matrix and r_xd = E[xd] is the cross-correlation vector. This is exactly the quadratic optimization problem Aw = b, with A = R_xx (PSD by construction) and b = r_xd. The solution w* = R_xx⁻¹ r_xd is the Wiener filter — a direct application of solving a PD linear system from quadratic optimization. The LMS algorithm is simply gradient descent on this quadratic: w_{k+1} = w_k + μ e_k x_k, where e_k = d_k − w_kᵀx_k is the current error — an online approximation to the Wiener filter using stochastic gradient descent.

Control Systems

Linear-quadratic regulators (LQR) minimize a quadratic cost J = Σ (xᵀQx + uᵀRu) over a control sequence u, subject to linear state dynamics x_{k+1} = Ax_k + Bu_k. This is a constrained quadratic program in the control variables. The solution — the LQR optimal controller — is obtained by solving the discrete algebraic Riccati equation, a nonlinear matrix equation whose solution P gives the optimal cost-to-go. The resulting controller is linear: u* = −Kx with K = (R + BᵀPB)⁻¹BᵀPA. This elegant result — a linear controller for a quadratic objective with linear constraints — is the crown jewel of linear systems theory and rests entirely on quadratic optimization via linear algebra.