CrunchAcademy · K-12

Math for Data Science · MDS-9

Mathematical Optimization

Turn 'find the best answer' into math you can solve — from convex sets and gradient descent to Lagrange multipliers, LPs, and the loss landscapes behind machine learning.

Optimization is how data science chooses the best model, the best decision, or the best allocation under constraints. This course builds the subject from the geometry of convex sets through first-order methods (gradient descent, SGD, momentum, Adam), Lagrangian duality and the KKT conditions, linear and quadratic programs, and the convex-optimization toolkit. Every idea is paired with hands-on solvers in Python (NumPy/SciPy/CVXPY) or Julia, ending with how optimization actually drives machine-learning training.

PrereqsMDS-2 Linear AlgebraMDS-4 Multivariable & Matrix Calculus
Code inPython (NumPy/SciPy/CVXPY)Julia
7Units
45Lessons
108Worked examples

Course Outline

MDS-9 — units

Every lesson is self-contained: a full explanation, worked step-by-step examples (in Python and R/Julia), common mistakes, a try-it-yourself, and an interactive quiz. Jump to a unit:

Weeks 1-2 Unit 1: Optimization Framing & Convexity
define-optimization-problemidentify-feasible-setdistinguish-local-globalrecognize-convex-setstest-function-convexityexplain-convexity-benefitmodel-standard-form
Lecture
What is an optimization problem? Objective, variables, constraints

Every modeling, fitting, and decision task in data science is secretly an optimization problem. Here we name its three pieces.

An optimization problem asks: choose decision variables x to make an objective function f(x) as small (or large) as possible, while respecting constraints. The variables are the knobs you control (model weights, allocations). The objective f(x) scores every candidate. Constraints g_i(x) <= 0 and h_j(x) = 0 carve out the feasible set of allowed choices. In ML, x is the parameter vector, f is the training loss, and constraints might cap weight norms or enforce probabilities summing to one. Recognizing these three pieces lets you hand a problem to a solver and reason about whether a solution exists, is unique, or is even attainable. Almost all of supervised learning, portfolio choice, and resource planning fits this template.

Worked Example 1

Problem. A bakery makes profit 3 per cake and 2 per pie. It has 60 units of flour; a cake uses 4, a pie uses 3. Identify objective, variables, constraints.

  1. Variables: c = cakes, p = pies, both >= 0.
  2. Objective: maximize 3c + 2p (profit).
  3. Resource constraint: 4c + 3p <= 60.
  4. Evaluate a candidate (c,p)=(9,8): `prof = 3*9 + 2*8` gives 43, flour `4*9+3*8 = 60` exactly feasible.

Answer. Variables (c,p)>=0; objective max 3c+2p; constraint 4c+3p<=60. Candidate (9,8) is feasible with profit 43.

Worked Example 2

Problem. Frame least-squares line fitting y ~ a x + b on data (xi, yi) as an optimization problem.

  1. Variables: slope a and intercept b.
  2. Objective: minimize sum of squared residuals f(a,b) = sum_i (a*xi + b - yi)^2.
  3. No constraints (unconstrained).
  4. In NumPy: `import numpy as np; A = np.c_[x, np.ones_like(x)]; coef,_,_,_ = np.linalg.lstsq(A, y, rcond=None)` returns (a, b) minimizing f.

Answer. min over (a,b) of sum_i (a xi + b - yi)^2; solved by np.linalg.lstsq giving the best-fit slope and intercept.

Common mistakes
  • Confusing the objective with a constraint. The thing you score and optimize is the objective; the rules you must obey are constraints.
  • Forgetting variable bounds like x >= 0. Implicit sign or range limits are real constraints that change the feasible set.
✎ Try it yourself

Problem. A factory makes profit 5 per unit of A and 4 per unit of B. Each A needs 2 labor-hours, each B needs 1; 100 labor-hours are available. Write objective, variables, constraints.

Solution. Variables a,b >= 0. Objective: maximize 5a + 4b. Constraint: 2a + b <= 100. The feasible set is the triangle bounded by a>=0, b>=0, 2a+b<=100. This is a linear program; its optimum will sit at a corner of that triangle (a foreshadow of Unit 5).

Minimize vs. maximize, feasible sets, and optimal value

Solvers expect minimization; here we convert maximize to minimize and pin down what 'feasible' and 'optimal value' mean.

By convention solvers minimize. Maximizing f(x) is identical to minimizing -f(x): argmax f = argmin(-f), and the optimal value flips sign. The feasible set is every x satisfying all constraints; if it is empty the problem is infeasible. The optimal value p* is the smallest objective over the feasible set; if the objective decreases without bound on the feasible set, p* = -infinity (unbounded). A minimizer x* is a point achieving p*. Distinguishing the optimal value (a number) from a minimizer (a point) matters because some problems have a finite value but no attaining point, or many minimizers with the same value. In ML this is the gap between the best achievable loss and the weights that realize it.

Worked Example 1

Problem. Convert: maximize f(x) = -x^2 + 4x over x in [0,3]. Give the equivalent minimization and the optimal value.

  1. Equivalent: minimize g(x) = x^2 - 4x over [0,3].
  2. g'(x) = 2x - 4 = 0 gives x = 2, inside [0,3].
  3. g(2) = 4 - 8 = -4, so min g = -4 and max f = +4.
  4. Check endpoints: g(0)=0, g(3)=-3; interior point x=2 wins.

Answer. Minimize x^2-4x; minimizer x*=2, optimal value of original maximize is f(2)=4.

Worked Example 2

Problem. Is minimize f(x)=x subject to x >= 1 bounded? What is the feasible set and optimal value?

  1. Feasible set: {x : x >= 1} = [1, infinity).
  2. f decreases as x decreases, so the smallest feasible x governs the optimum.
  3. Smallest feasible x is 1.
  4. `from scipy.optimize import linprog; res = linprog(c=[1], A_ub=[[-1]], b_ub=[-1])` returns res.fun = 1.0.

Answer. Feasible set [1, inf); optimal value p* = 1 at x* = 1; bounded.

Worked Example 3

Problem. Show minimize f(x) = -x over x >= 0 is unbounded.

  1. Feasible set is [0, infinity).
  2. f(x) = -x can be made arbitrarily negative by taking x large.
  3. As x -> infinity, f -> -infinity.
  4. No finite minimizer exists; solvers report 'unbounded'.

Answer. p* = -infinity; the problem is unbounded below, so no optimal point exists.

Common mistakes
  • Negating only the objective value but not flipping the sign back when reporting a maximum. After minimizing -f, the max of f is -(min value).
  • Assuming a finite optimal value always has an attaining point. Open feasible sets (e.g. x > 0) can have an infimum that is never reached.
✎ Try it yourself

Problem. Convert maximize 2x + 3y subject to x + y <= 5, x,y >= 0 into standard minimization form, and state the feasible set.

Solution. Minimize -(2x + 3y) = -2x - 3y subject to x + y <= 5, x >= 0, y >= 0. Feasible set is the triangle with corners (0,0), (5,0), (0,5). The optimum sits at a corner; (0,5) gives objective 15, the maximum. So min of -2x-3y is -15, attained at (0,5).

Global vs. local minima and why the difference matters

A method can get stuck at a point that only looks best nearby. Knowing the difference tells you whether to trust it.

A local minimum x* is best within some neighborhood: f(x*) <= f(x) for all x near x* in the feasible set. A global minimum is best over the entire feasible set. For nonconvex problems these differ — gradient methods may converge to a local minimum far worse than the global one. This matters enormously in ML: neural-network losses are nonconvex, so training finds a local (often good-enough) minimum, not provably the global one. The payoff of convexity, covered later, is that every local minimum is global, so any descent method that stops at a stationary point has found the best answer. Until then, you must reason about basins of attraction, restarts, and whether your landscape has one valley or many.

Worked Example 1

Problem. f(x) = x^4 - 3x^2 + 1. Find stationary points and classify global vs. local minima.

  1. f'(x) = 4x^3 - 6x = 2x(2x^2 - 3) = 0 gives x = 0 and x = +/- sqrt(1.5).
  2. f''(x) = 12x^2 - 6; at x=0, f''=-6 (local max); at x=+/-sqrt(1.5), f''=12*1.5-6=12 (minima).
  3. f(+/-sqrt(1.5)) = (1.5)^2 - 3(1.5) + 1 = 2.25 - 4.5 + 1 = -1.25, equal by symmetry.
  4. Both minima tie, so both are global; x=0 is a local max separating two basins.

Answer. Two global minima at x = +/- sqrt(1.5) with value -1.25; x=0 is a local maximum.

Worked Example 2

Problem. Run gradient descent on f(x)=x^4-3x^2+1 from x0=0.1 vs x0=-0.1 and observe which minimum each finds.

  1. Update: `x = x - 0.05*(4*x**3 - 6*x)`.
  2. From x0=0.1 (slightly positive), the negative gradient pushes x toward +sqrt(1.5)~1.225.
  3. From x0=-0.1, descent moves toward -sqrt(1.5)~-1.225.
  4. Each start lands in a different basin even though both are global here; with asymmetric functions one could be worse.

Answer. Starting point determines which basin/minimum gradient descent reaches; nonconvexity makes the result initialization-dependent.

Common mistakes
  • Assuming a converged gradient-descent point is the global optimum. On nonconvex problems it is only guaranteed local.
  • Ignoring initialization and restarts. For multi-basin landscapes, multiple random starts are needed to probe for a better minimum.
✎ Try it yourself

Problem. For f(x) = (x^2 - 1)^2, find all minima and say whether they are local, global, or both.

Solution. f'(x) = 2(x^2-1)(2x) = 4x(x^2-1) = 0 at x=0 and x=+/-1. f(0)=1, f(+/-1)=0. So x=+/-1 are minima with value 0 (both global, tied by symmetry) and x=0 is a local maximum with value 1. Any descent method lands at +1 or -1 depending on the starting side.

Convex sets and convex functions

Convexity is the single property that turns hard optimization into reliably solvable optimization. Define it precisely.

A set C is convex if, for any two points x, y in C, the whole segment tx + (1-t)y (for t in [0,1]) stays in C — no dents or holes. A function f is convex if its domain is convex and its graph lies on or below every chord: f(tx + (1-t)y) <= t f(x) + (1-t) f(y). Equivalently the region above the graph (the epigraph) is a convex set. Convex functions bowl upward, so they have a single connected set of minimizers. This is why convexity is the dividing line in optimization: convex problems have no spurious local minima, and efficient algorithms find the global optimum. Norms, squared error, log-sum-exp, and L2-regularized linear losses are all convex.

Worked Example 1

Problem. Is the set {(x,y) : x^2 + y^2 <= 1} convex? Verify numerically.

  1. Take two interior/boundary points, e.g. p=(1,0), q=(0,1).
  2. Midpoint m = (0.5, 0.5); check `0.5**2 + 0.5**2 = 0.5 <= 1`. Inside.
  3. General: a ball is convex because the norm is convex; any chord stays within radius 1.
  4. `import numpy as np; t=np.linspace(0,1,50); pts=np.outer(t,[1,0])+np.outer(1-t,[0,1]); (np.sum(pts**2,1)<=1).all()` -> True.

Answer. Yes, the unit disk is convex; every segment between its points stays inside.

Worked Example 2

Problem. Check whether f(x) = |x| is convex using the chord inequality at x=-2, y=4, t=0.5.

  1. LHS: f(0.5*(-2)+0.5*4) = f(1) = 1.
  2. RHS: 0.5*f(-2) + 0.5*f(4) = 0.5*2 + 0.5*4 = 3.
  3. 1 <= 3, inequality holds here.
  4. It holds for all x,y because absolute value is a norm; `f=lambda x: abs(x)` is convex everywhere.

Answer. f(x)=|x| is convex (chord lies above the graph for every pair of points).

Worked Example 3

Problem. Show the set {(x,y) : xy >= 1, x > 0} is NOT convex.

  1. Pick p=(0.5, 2) (product 1) and q=(2, 0.5) (product 1), both in the set.
  2. Midpoint m=(1.25, 1.25); product = 1.5625 >= 1 here, so try a different pair.
  3. Pick p=(0.1,10), q=(10,0.1); midpoint (5.05,5.05) product 25.5 ok — but the region is the union of a curved boundary; instead test p=(0.5,2),q=(3,0.34): segment dips below xy=1.
  4. The boundary xy=1 is curved (a hyperbola); its feasible region above it is convex, but the full constraint set with x>0 only is convex; the canonical nonconvex example is {xy>=1} WITHOUT sign restriction, which splits into two branches whose connecting segment leaves the set.

Answer. Without the x>0 restriction, {xy >= 1} has two disconnected branches; a segment from one branch to the other leaves the set, so it is nonconvex.

Common mistakes
  • Thinking any 'round-looking' set is convex. A set with a hole, a dent, or two disconnected pieces is not convex.
  • Confusing convex and concave. A concave function curves the other way (chord below graph); -f is convex iff f is concave.
✎ Try it yourself

Problem. Is f(x,y) = x^2 + y^2 convex? Verify the chord inequality at p=(0,0), q=(2,2), t=0.5.

Solution. LHS: f(0.5*(0,0)+0.5*(2,2)) = f(1,1) = 1+1 = 2. RHS: 0.5*f(0,0) + 0.5*f(2,2) = 0.5*0 + 0.5*8 = 4. Since 2 <= 4, the inequality holds. The function is a paraboloid bowl; its Hessian is 2I (positive definite), so it is convex everywhere.

First- and second-order conditions for convexity

Two practical tests tell you whether a differentiable function is convex without checking every chord.

For a differentiable f, the first-order condition says f is convex iff its graph lies above every tangent: f(y) >= f(x) + grad f(x).(y - x) for all x, y. The tangent line/plane underestimates the function. For a twice-differentiable f, the second-order condition is cleaner: f is convex iff its Hessian H(x) is positive semidefinite (all eigenvalues >= 0) everywhere on its domain; in one variable this is just f''(x) >= 0. Strict positive definiteness gives strict convexity and a unique minimizer. These tests are the workhorse for proving a loss is convex — checking the Hessian of squared error or logistic loss is how you certify that gradient descent will find the global optimum.

Worked Example 1

Problem. Use the second-order test to confirm f(x) = e^x is convex.

  1. f'(x) = e^x, f''(x) = e^x.
  2. e^x > 0 for all x.
  3. f'' > 0 everywhere means strictly convex.
  4. `import numpy as np; x=np.linspace(-3,3,100); (np.exp(x) > 0).all()` -> True.

Answer. e^x is strictly convex since f''(x) = e^x > 0 for all x.

Worked Example 2

Problem. Is f(x,y) = x^2 + 3y^2 - xy convex? Check the Hessian.

  1. Gradient: (2x - y, 6y - x). Hessian H = [[2, -1], [-1, 6]].
  2. Eigenvalues via `import numpy as np; np.linalg.eigvalsh([[2,-1],[-1,6]])` -> approximately [1.75, 6.25].
  3. Both eigenvalues > 0, so H is positive definite.
  4. Constant Hessian that is PD everywhere means strictly convex.

Answer. Yes; the Hessian [[2,-1],[-1,6]] has positive eigenvalues (~1.75, ~6.25), so f is strictly convex.

Worked Example 3

Problem. Show f(x,y) = x^2 - y^2 is NOT convex.

  1. Hessian H = [[2, 0], [0, -2]].
  2. Eigenvalues are 2 and -2.
  3. A negative eigenvalue means H is indefinite, not PSD.
  4. `np.linalg.eigvalsh([[2,0],[0,-2]])` -> [-2, 2]; the -2 disqualifies convexity.

Answer. Not convex; the Hessian is indefinite (eigenvalue -2), giving a saddle shape.

Common mistakes
  • Checking the Hessian at one point only. Convexity needs the Hessian PSD over the whole domain, not just at a candidate.
  • Confusing positive determinant with positive definiteness. A 2x2 Hessian with positive determinant but negative trace is negative definite, not PSD.
✎ Try it yourself

Problem. Determine whether f(x,y) = 2x^2 + 2y^2 + 2xy is convex using the Hessian eigenvalue test.

Solution. Hessian H = [[4, 2], [2, 4]]. Eigenvalues solve det(H - lambda I)=0: (4-lambda)^2 - 4 = 0, so 4-lambda = +/-2, giving lambda = 2 and 6. Both positive, so H is positive definite and f is strictly convex. In code: np.linalg.eigvalsh([[4,2],[2,4]]) -> [2., 6.].

Why convexity makes optimization tractable

This lesson states the central payoff: for convex problems, 'local' equals 'global' and solvers are reliable.

The defining benefit of convexity: in a convex optimization problem (convex objective over a convex feasible set), every local minimum is a global minimum, and the set of optimal points is itself convex. Therefore any algorithm that reaches a stationary point — where the gradient is zero, or KKT conditions hold — has found the global optimum; there are no spurious traps. This makes convex problems polynomial-time solvable to high accuracy with interior-point or first-order methods, and it is why disciplined-convex tools like CVXPY can guarantee a correct answer. In ML, framing a model so its training loss is convex (linear/logistic regression, SVM, lasso, ridge) means optimization is a solved problem; the difficulty in deep learning is precisely that its loss is nonconvex.

Worked Example 1

Problem. Explain why minimizing the convex f(x) = (x-3)^2 needs no global search.

  1. f''(x) = 2 > 0, so f is strictly convex.
  2. Set f'(x) = 2(x-3) = 0 -> x = 3, the unique stationary point.
  3. By convexity that stationary point is the global minimum.
  4. `from scipy.optimize import minimize; minimize(lambda x:(x-3)**2, x0=100).x` -> ~3.0 regardless of start.

Answer. x* = 3 is the unique global minimum; any starting point converges there because the problem is convex.

Worked Example 2

Problem. Contrast: minimize convex 0.5*||Ax - b||^2 vs. a nonconvex neural loss, in terms of guarantees.

  1. The least-squares objective has Hessian A^T A which is PSD, so it is convex.
  2. Any solver (normal equations, gradient descent) reaches the same global optimum: `x = np.linalg.lstsq(A,b,rcond=None)[0]`.
  3. A neural net's loss is nonconvex; different initializations yield different minima.
  4. Only the convex case carries a global-optimality certificate.

Answer. The convex least-squares problem has a guaranteed unique global solution; the nonconvex net loss does not.

Common mistakes
  • Believing convexity guarantees a unique minimizer. It guarantees local=global, but a flat-bottomed convex function can have many minimizers (a convex set of them).
  • Assuming any smooth function with one visible valley is convex. You must verify with the Hessian; appearances mislead in higher dimensions.
✎ Try it yourself

Problem. A loss L(w) is shown to have a positive-semidefinite Hessian everywhere. A gradient-descent run stops where the gradient is essentially zero. What can you conclude about that point?

Solution. Because the Hessian is PSD everywhere, L is convex. For a convex function a point with zero gradient is a global minimizer. So the gradient-descent stopping point is a global optimum (one of possibly several if the Hessian is only PSD, not PD). No restarts or global search are needed to certify optimality.

Modeling a real problem in standard form

Solvers expect a canonical layout. This lesson turns a word problem into minimize-subject-to standard form.

Standard form writes every problem as: minimize f(x) subject to g_i(x) <= 0 (inequalities) and h_j(x) = 0 (equalities). To convert a real problem: (1) name the decision variables, (2) write the objective and negate it if it was a maximization, (3) move every inequality to the form (something) <= 0, and (4) write equalities as (something) = 0. Variable bounds like x >= 0 become -x <= 0. Putting a problem in standard form is the interface to solvers like scipy.optimize and CVXPY, and it makes the structure (linear? quadratic? convex?) visible so you can pick the right algorithm. This discipline prevents sign errors and ensures the solver optimizes what you intend.

Worked Example 1

Problem. A diet must supply at least 50 units of protein. Food A gives 10/unit at cost 2; food B gives 5/unit at cost 1. Write in standard form.

  1. Variables: xA, xB >= 0 (servings).
  2. Objective (minimize cost): minimize 2*xA + 1*xB.
  3. Protein constraint 10 xA + 5 xB >= 50 -> standard form 50 - 10 xA - 5 xB <= 0.
  4. Bounds: -xA <= 0, -xB <= 0.

Answer. minimize 2xA + xB s.t. 50 - 10xA - 5xB <= 0, -xA <= 0, -xB <= 0.

Worked Example 2

Problem. Maximize return 0.08 a + 0.05 b with a + b = 1 (full investment), a,b >= 0. Put in standard form and solve.

  1. Negate objective: minimize -(0.08 a + 0.05 b).
  2. Equality: a + b - 1 = 0. Bounds: -a <= 0, -b <= 0.
  3. `from scipy.optimize import linprog; res = linprog(c=[-0.08,-0.05], A_eq=[[1,1]], b_eq=[1], bounds=[(0,None),(0,None)])`.
  4. res.x -> [1, 0]; put everything in the higher-return asset.

Answer. minimize -0.08a-0.05b s.t. a+b-1=0, a,b>=0; solution a=1,b=0, max return 0.08.

Common mistakes
  • Leaving a constraint as >= without flipping it. Standard form needs <= 0; rewrite g(x) >= c as c - g(x) <= 0.
  • Forgetting to negate the objective for a maximization. Solvers minimize, so a max becomes minimize -f.
✎ Try it yourself

Problem. A shop maximizes profit 6x + 4y subject to x + y <= 10 and x <= 6, with x,y >= 0. Write it in standard minimization form.

Solution. Minimize -(6x + 4y) = -6x - 4y subject to: x + y - 10 <= 0, x - 6 <= 0, -x <= 0, -y <= 0. This is a linear program; in scipy: linprog(c=[-6,-4], A_ub=[[1,1],[1,0]], b_ub=[10,6], bounds=[(0,None),(0,None)]), which returns x=6, y=4 with profit 52.

Key terms
  • Objective function — the quantity f(x) we minimize or maximize over the decision variables.
  • Feasible set — the set of points x satisfying all constraints; optimization searches only within it.
  • Optimal value — the smallest (or largest) objective value attainable over the feasible set.
  • Local minimum — a point that is best within a small neighborhood but not necessarily overall.
  • Global minimum — a point with the smallest objective value over the entire feasible set.
  • Convex set — a set where the line segment between any two members stays entirely inside the set.
  • Convex function — a function whose graph lies below any chord, equivalently f(tx+(1-t)y) <= t f(x)+(1-t) f(y).
  • Standard form — a canonical way to write a problem (minimize f subject to g_i(x) <= 0, h_j(x) = 0) that solvers expect.
Assignment · Spot the convexity

Pick five functions of one or two variables (mix convex, concave, and neither). In Python with NumPy/Matplotlib, plot each and numerically check the chord/midpoint inequality and the sign of the second derivative or Hessian eigenvalues. Then take one small real decision (for example, a budget allocation) and write it in standard form. Argue whether each function is convex and whether your modeled problem is a convex program.

Deliverable · A notebook with plots, your convexity verdict for each function with evidence, and the standard-form statement of your modeled problem.

Quiz · 5 questions
  1. 1. A set is convex if, for any two points in the set:

  2. 2. For a convex optimization problem, any local minimum is:

  3. 3. A twice-differentiable function of one variable is convex on an interval when:

  4. 4. Which is the standard form of a minimization problem?

  5. 5. Maximizing f(x) is equivalent to:

You'll be able to

I can express a word problem as a standard-form optimization with objective, variables, and constraints.

I can decide whether a set or function is convex and explain why convexity guarantees a global optimum.

I can distinguish local from global minima and identify when they coincide.

Weeks 3-4 Unit 2: Unconstrained Optimization
apply-first-order-conditionsclassify-with-hessianuse-line-searchapply-newton-methodrecognize-quasi-newtonuse-scipy-minimize
Lecture
Optimality conditions: stationary points and the gradient

Before running any solver, calculus tells you where minima can hide: where the gradient vanishes.

For a smooth, unconstrained function f, a necessary condition for an interior local minimum at x* is that the gradient is zero: grad f(x*) = 0. Such points are stationary points. The gradient points in the direction of steepest increase, so at a minimum there is no direction of first-order decrease left. This first-order condition is necessary but not sufficient — maxima and saddle points are also stationary — so it only narrows candidates. In ML, training a model means driving the gradient of the loss toward zero; convergence criteria check ||grad|| < tolerance. Finding stationary points by solving grad f = 0 analytically works for small problems and underlies every gradient-based optimizer used at scale.

Worked Example 1

Problem. Find the stationary point of f(x,y) = x^2 + y^2 - 2x - 4y + 5.

  1. grad f = (2x - 2, 2y - 4).
  2. Set each component to zero: 2x-2=0 -> x=1; 2y-4=0 -> y=2.
  3. Stationary point is (1, 2).
  4. `import numpy as np; from scipy.optimize import minimize; minimize(lambda v:(v[0]-1)**2+(v[1]-2)**2, [0,0]).x` -> ~[1,2].

Answer. Unique stationary point (1, 2); value f(1,2) = 0.

Worked Example 2

Problem. For the squared-error loss L(w) = sum_i (w*xi - yi)^2 with x=[1,2,3], y=[2,4,5], find the stationary w.

  1. L'(w) = sum_i 2 xi (w xi - yi) = 0.
  2. Solve w * sum xi^2 = sum xi yi.
  3. sum xi^2 = 1+4+9 = 14; sum xi yi = 2+8+15 = 25.
  4. `w = 25/14` -> ~1.786, the least-squares slope through the origin.

Answer. Stationary (optimal) w = 25/14 ~ 1.786; gradient is zero there.

Common mistakes
  • Treating grad f = 0 as sufficient for a minimum. It also flags maxima and saddles; you must check second-order info.
  • Forgetting that the condition is for interior points. At a boundary, the minimum can occur where the gradient is nonzero.
✎ Try it yourself

Problem. Find the stationary point of f(x,y) = 3x^2 + 2y^2 + 12x - 8y + 1.

Solution. grad f = (6x + 12, 4y - 8). Set to zero: 6x+12=0 -> x=-2; 4y-8=0 -> y=2. Stationary point (-2, 2). Since both pure second derivatives are positive and there is no cross term, the Hessian diag(6,4) is positive definite, so it is a minimum with value f(-2,2)=3*4+2*4-24-16+1 = 12+8-24-16+1 = -19.

Second-order conditions and the Hessian test

Once you have a stationary point, the Hessian decides whether it is a minimum, maximum, or saddle.

At a stationary point x*, the Hessian H (matrix of second partials) classifies it via its eigenvalues. If H is positive definite (all eigenvalues > 0), x* is a strict local minimum — the function curves up in every direction. If H is negative definite, it is a local maximum. If H has both positive and negative eigenvalues, it is a saddle point. If H is positive semidefinite with a zero eigenvalue, the test is inconclusive. This second-order sufficient condition is how you confirm a candidate is genuinely a minimizer. For ML losses, checking the Hessian's definiteness certifies convexity locally and predicts how fast curvature-aware methods like Newton's will converge.

Worked Example 1

Problem. Classify the stationary point of f(x,y) = x^2 + y^2 - xy.

  1. grad = (2x - y, 2y - x); zero at (0,0).
  2. Hessian H = [[2, -1], [-1, 2]].
  3. `import numpy as np; np.linalg.eigvalsh([[2,-1],[-1,2]])` -> [1, 3].
  4. Both eigenvalues positive -> positive definite.

Answer. (0,0) is a strict local (and global) minimum; Hessian eigenvalues 1 and 3 are positive.

Worked Example 2

Problem. Classify the stationary point of f(x,y) = x^2 - y^2.

  1. grad = (2x, -2y); zero at (0,0).
  2. Hessian H = [[2, 0], [0, -2]].
  3. Eigenvalues 2 and -2 (one positive, one negative).
  4. Indefinite Hessian indicates a saddle.

Answer. (0,0) is a saddle point: minimum along x, maximum along y.

Worked Example 3

Problem. For the Rosenbrock function f(x,y) = (1-x)^2 + 100(y - x^2)^2, classify the point (1,1).

  1. Gradient is zero at (1,1) (the known minimizer).
  2. Hessian at (1,1): H = [[2 - 400(y - 3x^2), -400x], [-400x, 200]] = [[802, -400], [-400, 200]].
  3. `import numpy as np; np.linalg.eigvalsh([[802,-400],[-400,200]])` -> approx [0.4, 1001.6], both > 0.
  4. Positive definite (though ill-conditioned).

Answer. (1,1) is the strict local minimum; Hessian is PD but highly ill-conditioned, foreshadowing slow plain gradient descent.

Common mistakes
  • Reading definiteness off the diagonal alone. Off-diagonal terms matter; use eigenvalues or leading-minor tests.
  • Calling an inconclusive (PSD with a zero eigenvalue) case a confirmed minimum. The second-order test is silent there.
✎ Try it yourself

Problem. Classify the stationary point of f(x,y) = 2x^2 + y^2 + 2xy at (0,0).

Solution. grad = (4x + 2y, 2y + 2x); zero at (0,0). Hessian H = [[4, 2], [2, 2]]. Eigenvalues solve (4-l)(2-l) - 4 = 0 -> l^2 - 6l + 4 = 0 -> l = 3 +/- sqrt(5) ~ 0.76 and 5.24, both positive. So H is positive definite and (0,0) is a strict local minimum.

Line search and step-size selection

Knowing which direction to move is half the battle; line search picks how far to step.

Given a descent direction d at point x (one with grad f . d < 0), a line search chooses the step length alpha in the update x <- x + alpha d. Exact line search minimizes f along the ray, but it is usually too costly. Backtracking (Armijo) line search instead starts with a guess alpha and shrinks it by a factor until a sufficient-decrease condition holds: f(x + alpha d) <= f(x) + c1 * alpha * grad f . d, with c1 small (e.g. 1e-4). The Wolfe conditions add a curvature requirement to avoid steps that are too short. Good step-size selection is what keeps gradient and Newton methods stable; in ML this is the principled cousin of tuning a learning rate.

Worked Example 1

Problem. Backtracking on f(x) = x^2 from x=2 with direction d = -f'(2) = -4, start alpha=1, factor 0.5, c1=1e-4.

  1. f(2)=4, grad=4, d=-4, grad.d = -16.
  2. Try alpha=1: x_new = 2 + 1*(-4) = -2, f=4. Need 4 <= 4 + 1e-4*1*(-16)=3.9984. Fails.
  3. Shrink alpha=0.5: x_new = 2 - 2 = 0, f=0. Need 0 <= 4 + 0.5*(-16)*1e-4 ~ 3.9992. Holds.
  4. Accept alpha=0.5; one step lands exactly at the minimum x=0.

Answer. Backtracking selects alpha = 0.5, moving from x=2 to the minimizer x=0.

Worked Example 2

Problem. Implement Armijo backtracking generically in Python and apply to f(x)=(x-3)^2 from x0=0.

  1. `f=lambda x:(x-3)**2; g=lambda x:2*(x-3)`
  2. `x=0.0; d=-g(x); a=1.0` # d = 6
  3. `while f(x+a*d) > f(x) + 1e-4*a*g(x)*d: a*=0.5`
  4. After loop, `x = x + a*d`; with a=1 the sufficient-decrease holds (f(6)=9 fails, so a halves to 0.5 -> x=3, f=0 passes).

Answer. Backtracking finds alpha=0.5, stepping to x=3, the exact minimum of (x-3)^2.

Common mistakes
  • Using a fixed step that ignores the function's scale. Too big diverges, too small crawls; line search adapts automatically.
  • Stepping along a direction that is not a descent direction. If grad f . d >= 0 the Armijo test can never be satisfied for small alpha.
✎ Try it yourself

Problem. From x=4 on f(x)=x^2 with direction d=-f'(4), starting alpha=1 and halving, find the accepted Armijo step (c1=1e-4).

Solution. f(4)=16, f'(4)=8, d=-8, grad.d=-64. alpha=1: x_new=4-8=-4, f=16; need 16 <= 16 + 1e-4*(-64)=15.9936, fails. alpha=0.5: x_new=4-4=0, f=0; need 0 <= 16 + 0.5*(-64)*1e-4 ~ 15.9968, holds. Accept alpha=0.5, landing at x=0, the minimum.

Newton's method and its quadratic convergence

Newton's method uses curvature to take smart steps, converging dramatically fast near the optimum.

Newton's method models f locally as a quadratic using the Hessian and jumps to that quadratic's minimum. The update is x <- x - H(x)^{-1} grad f(x). Whereas gradient descent only knows the slope, Newton rescales the step by the inverse Hessian, automatically correcting for ill-conditioning and pointing toward the true minimum. Near a minimizer with positive-definite Hessian it converges quadratically — the number of correct digits roughly doubles each iteration. The costs are computing and inverting the Hessian (O(n^3)) and needing it positive definite to guarantee descent. For small-to-medium problems and the inner loops of many ML methods (e.g. Newton steps in logistic regression / IRLS), this speed is worth it.

Worked Example 1

Problem. One Newton step on f(x) = x^2 - 4x + 4 from x0 = 0.

  1. f'(x) = 2x - 4, f''(x) = 2.
  2. Newton update: x1 = x0 - f'(x0)/f''(x0) = 0 - (-4)/2 = 2.
  3. f'(2) = 0, so x1 = 2 is the exact minimizer.
  4. Quadratic functions are solved in a single Newton step.

Answer. x1 = 2 in one step; Newton finds the minimum of a quadratic immediately.

Worked Example 2

Problem. Two Newton steps on f(x) = x^4 from x0 = 1 (note f'' vanishes at the min, slowing convergence).

  1. f'(x)=4x^3, f''(x)=12x^2.
  2. x1 = 1 - (4*1)/(12*1) = 1 - 1/3 = 0.6667.
  3. x2 = 0.6667 - (4*0.6667^3)/(12*0.6667^2) = 0.6667 - 0.6667/3 = 0.4444.
  4. Convergence here is only linear (ratio 2/3) because the Hessian degenerates at x=0; classic quadratic speed needs nonzero curvature at the optimum.

Answer. x2 ~ 0.444; Newton still converges but only linearly when f''=0 at the minimum.

Worked Example 3

Problem. Newton in 2D on f(x,y)=x^2+10y^2 from (1,1).

  1. grad = (2x, 20y) = (2, 20) at (1,1).
  2. Hessian H = [[2,0],[0,20]], H^{-1} = [[0.5,0],[0,0.05]].
  3. Step: (1,1) - H^{-1} grad = (1,1) - (0.5*2, 0.05*20) = (1,1) - (1,1) = (0,0).
  4. `import numpy as np; H=np.array([[2,0],[0,20]]); g=np.array([2,20]); np.array([1,1]) - np.linalg.solve(H,g)` -> [0,0].

Answer. Reaches the minimizer (0,0) in one step; Newton neutralizes the 10x conditioning that would slow gradient descent.

Common mistakes
  • Applying Newton when the Hessian is not positive definite. The 'step' may point uphill or toward a saddle; modified/damped Newton is needed.
  • Ignoring cost. Forming and inverting H is O(n^3); for high-dimensional ML, quasi-Newton or first-order methods are cheaper.
✎ Try it yourself

Problem. Take one Newton step on f(x) = x^2 + 2x + 1 from x0 = 3.

Solution. f'(x) = 2x + 2, f''(x) = 2. Newton step: x1 = 3 - f'(3)/f''(3) = 3 - (8)/2 = 3 - 4 = -1. Check f'(-1) = 0, so x1 = -1 is the exact minimizer (f = (x+1)^2). Newton solves any quadratic in a single step.

Quasi-Newton methods (BFGS) at a glance

BFGS gets Newton-like speed without ever forming the Hessian — the workhorse of scipy.optimize.

Quasi-Newton methods approximate the (inverse) Hessian using only gradient information gathered across iterations. BFGS, the most popular, updates an estimate B_k of the Hessian (or H_k of its inverse) from successive parameter changes s_k = x_{k+1} - x_k and gradient changes y_k = grad_{k+1} - grad_k, enforcing the secant equation B_{k+1} s_k = y_k. The step is then x <- x - H_k grad, like Newton but with H built up cheaply (O(n^2) per step, no second derivatives). BFGS achieves superlinear convergence and is the default for smooth unconstrained problems in scipy.optimize.minimize. For large-scale ML, the memory-limited L-BFGS variant stores only a few vectors, making it practical for logistic regression and conditional random fields.

Worked Example 1

Problem. Solve the Rosenbrock minimum with BFGS in scipy and report the iteration count.

  1. `from scipy.optimize import minimize; import numpy as np`
  2. `rosen = lambda v: (1-v[0])**2 + 100*(v[1]-v[0]**2)**2`
  3. `res = minimize(rosen, [-1.2, 1.0], method='BFGS')`
  4. `res.x` -> ~[1,1], `res.nit` -> roughly 20-40 iterations without any explicit Hessian.

Answer. BFGS converges to (1,1) in a few dozen iterations using only gradients (estimated numerically).

Worked Example 2

Problem. Why use L-BFGS over BFGS for a 1-million-parameter logistic regression?

  1. BFGS stores a dense n-by-n inverse-Hessian estimate: 1e6 x 1e6 is infeasible (terabytes).
  2. L-BFGS keeps only the last m (e.g. 10) pairs (s_k, y_k) and reconstructs the action of H implicitly.
  3. Memory is O(m n) instead of O(n^2).
  4. `minimize(loss, w0, jac=grad, method='L-BFGS-B')` scales to such sizes.

Answer. L-BFGS uses O(m n) memory by storing a short history, making quasi-Newton feasible at ML scale.

Common mistakes
  • Expecting BFGS to match Newton's exact quadratic step early on. The Hessian estimate needs several iterations to warm up.
  • Using full BFGS in very high dimensions. The dense O(n^2) inverse-Hessian estimate blows up memory; switch to L-BFGS.
✎ Try it yourself

Problem. Minimize f(x,y) = (x-2)^2 + (y+3)^2 with scipy BFGS and state the expected minimizer.

Solution. from scipy.optimize import minimize; res = minimize(lambda v:(v[0]-2)**2+(v[1]+3)**2, [0,0], method='BFGS'). The objective is a simple convex bowl with minimum at (2, -3). BFGS converges there in very few iterations (the gradient (2(x-2), 2(y+3)) vanishes at (2,-3)). res.x is approximately [2.0, -3.0] and res.fun approximately 0.

Solving small problems with scipy.optimize.minimize

Tie it together: a single unified interface drives gradient descent, Newton, and quasi-Newton solvers.

scipy.optimize.minimize is the standard entry point for unconstrained (and some constrained) smooth optimization. You pass the objective, a starting point x0, optionally the gradient via jac and Hessian via hess, and choose a method: 'BFGS' or 'L-BFGS-B' (quasi-Newton, gradient-based), 'Newton-CG' or 'trust-ncg' (use the Hessian), or 'Nelder-Mead' (derivative-free). It returns an OptimizeResult with .x (the minimizer), .fun (optimal value), .nit (iterations), and .success. Supplying an analytic gradient speeds convergence and accuracy versus finite differences. Knowing how to call this, read the result, and pick a method by problem size and smoothness is the practical skill that turns optimization theory into working code for fitting models and tuning small systems.

Worked Example 1

Problem. Minimize f(x) = (x-5)^2 + 2 with an analytic gradient using minimize.

  1. `from scipy.optimize import minimize`
  2. `f = lambda x: (x[0]-5)**2 + 2; g = lambda x: [2*(x[0]-5)]`
  3. `res = minimize(f, x0=[0.0], jac=g, method='BFGS')`
  4. `res.x` -> ~[5.0], `res.fun` -> ~2.0, `res.success` -> True.

Answer. Minimizer x ~ 5, optimal value ~ 2; supplying jac makes it converge in a couple of iterations.

Worked Example 2

Problem. Fit a least-squares slope by minimizing loss with Newton-CG, providing gradient and Hessian.

  1. Data x=[1,2,3], y=[2,4,5]; loss L(w)=sum (w*xi - yi)^2.
  2. `import numpy as np; X=np.array([1,2,3]); Y=np.array([2,4,5])`
  3. `L=lambda w: np.sum((w*X-Y)**2); g=lambda w: np.array([np.sum(2*X*(w*X-Y))]); H=lambda w: np.array([[np.sum(2*X**2)]])`
  4. `from scipy.optimize import minimize; minimize(L,[0.0],jac=g,hess=H,method='Newton-CG').x` -> ~[1.786] = 25/14.

Answer. Newton-CG returns w ~ 1.786 = 25/14, matching the closed-form least-squares slope.

Common mistakes
  • Ignoring res.success and reading res.x blindly. A failed or non-converged run can return a meaningless point; always check the status.
  • Forgetting that finite-difference gradients are noisy and slow. Pass an analytic jac when you can for speed and accuracy.
✎ Try it yourself

Problem. Use scipy.optimize.minimize to find the minimum of f(x,y) = (x-1)^2 + (y-2)^2 + 3, and predict .x and .fun.

Solution. from scipy.optimize import minimize; res = minimize(lambda v:(v[0]-1)**2+(v[1]-2)**2+3, [0,0], method='BFGS'). The bowl is centered at (1,2) with floor value 3. The gradient (2(x-1), 2(y-2)) vanishes at (1,2). So res.x is approximately [1.0, 2.0] and res.fun is approximately 3.0, with res.success True.

Key terms
  • Stationary point — a point where the gradient is zero; a necessary condition for an interior minimum.
  • Gradient — the vector of partial derivatives pointing in the direction of steepest increase of f.
  • Hessian — the matrix of second partial derivatives; its definiteness classifies a stationary point.
  • Second-order condition — a stationary point is a local min when the Hessian is positive definite there.
  • Line search — choosing the step size along a search direction, e.g. by satisfying the Armijo/Wolfe conditions.
  • Newton's method — uses the Hessian to take steps -H^{-1} g, giving fast quadratic local convergence.
  • Quasi-Newton (BFGS) — approximates the Hessian from gradient changes, avoiding costly second-derivative computation.
  • Descent direction — a direction d with gradient dot d < 0, along which the objective initially decreases.
Assignment · Newton vs. gradient steps on a bowl

Take a smooth convex test function such as a 2D quadratic or the Rosenbrock function. In Python, implement a simple backtracking line search and run both a basic steepest-descent loop and a Newton step loop from the same starting point. Track the objective and gradient norm each iteration. Then solve the same problem with scipy.optimize.minimize using BFGS and Newton-CG and compare iteration counts.

Deliverable · A notebook with convergence plots (objective and gradient norm vs. iteration) comparing your steepest-descent, Newton, and SciPy solver runs, plus a short note on why Newton converged faster.

Quiz · 5 questions
  1. 1. A necessary condition for x* to be an interior local minimum of a smooth f is:

  2. 2. At a stationary point, a positive definite Hessian indicates a:

  3. 3. Newton's method differs from steepest descent because it:

  4. 4. What problem do quasi-Newton methods like BFGS solve?

  5. 5. A backtracking line search is used to:

You'll be able to

I can use the gradient and Hessian to find and classify candidate minimizers.

I can explain how line search and Newton's method choose directions and step sizes.

I can solve an unconstrained problem with scipy.optimize and interpret the result.

Weeks 5-6 Unit 3: Gradient Descent & Variants
implement-batch-gdanalyze-convergenceimplement-sgdapply-momentumapply-adamtune-learning-ratebuild-optimizer
Lecture
Batch gradient descent and the learning rate

The single most important algorithm in ML training: step downhill against the full-data gradient.

Batch gradient descent minimizes a loss L(w) by repeatedly moving against the gradient: w <- w - eta * grad L(w), where eta is the learning rate (step size) and the gradient is computed over the entire training set. Each step reduces the loss provided eta is small enough; for a smooth loss with gradient Lipschitz constant beta, any eta < 2/beta converges. The learning rate is the central knob: too large and the iterates overshoot and diverge; too small and training crawls. In ML this exact rule trains linear and logistic regression and is the backbone of deep-learning optimizers. Convergence to the global optimum is guaranteed when L is convex; otherwise it finds a stationary point of whatever basin it starts in.

Worked Example 1

Problem. Two batch GD steps on L(w) = (w-3)^2 from w=0 with eta=0.1.

  1. grad L(w) = 2(w-3).
  2. Step 1: w = 0 - 0.1*2*(0-3) = 0 + 0.6 = 0.6.
  3. Step 2: w = 0.6 - 0.1*2*(0.6-3) = 0.6 + 0.48 = 1.08.
  4. Iterates march toward the minimizer w*=3; `w=0\nfor _ in range(50): w-=0.1*2*(w-3)` -> ~3.0.

Answer. After two steps w ~ 1.08, converging geometrically to w* = 3.

Worked Example 2

Problem. Batch GD for linear regression on x=[1,2,3], y=[2,4,5], one step from w=0, b=0, eta=0.01.

  1. Predictions p = w*x + b = [0,0,0]; residuals r = p - y = [-2,-4,-5].
  2. grad_w = (2/n) sum(r*x) = (2/3)(-2-8-15) = (2/3)(-25) = -16.67; grad_b = (2/3)sum(r) = (2/3)(-11) = -7.33.
  3. Update: w = 0 - 0.01*(-16.67) = 0.167; b = 0 - 0.01*(-7.33) = 0.073.
  4. `import numpy as np; x=np.array([1,2,3.]); y=np.array([2,4,5.]); w=b=0\nfor _ in range(2000):\n r=(w*x+b)-y; w-=0.01*(2/3)*np.sum(r*x); b-=0.01*(2/3)*np.sum(r)` -> w~1.5, b~0.667.

Answer. One step gives w~0.167, b~0.073; iterating converges to the least-squares fit w~1.5, b~0.667.

Common mistakes
  • Picking the learning rate by guesswork without scaling for the data. Standardizing features and tuning eta is essential to avoid divergence.
  • Confusing batch GD with SGD. Batch GD uses the full dataset for every step; the gradient is exact, not sampled.
✎ Try it yourself

Problem. Perform two batch GD steps on L(w) = (w-10)^2 from w=0 with eta=0.2. Where is it heading?

Solution. grad = 2(w-10). Step 1: w = 0 - 0.2*2*(0-10) = 0 + 4 = 4. Step 2: w = 4 - 0.2*2*(4-10) = 4 + 2.4 = 6.4. The contraction factor is (1 - 2*eta) = 0.6 per step on (w-10), so w -> 10. After two steps w = 6.4, converging to the minimizer w* = 10.

Convergence, step sizes, and conditioning

Why some losses train fast and others zigzag forever — it comes down to curvature and the step size.

For a convex quadratic with Hessian H, gradient descent converges at a rate governed by the condition number kappa = lambda_max / lambda_min of H. With the optimal fixed step, the error contracts by roughly (kappa - 1)/(kappa + 1) per iteration: a well-conditioned problem (kappa near 1) converges in a few steps, while an ill-conditioned one (large kappa) zigzags slowly across a narrow valley. The largest stable step is eta < 2/lambda_max; beyond it the iterates diverge. This is why feature scaling, normalization, and preconditioning matter in ML — they shrink kappa so plain gradient descent trains quickly. It also motivates momentum and adaptive methods, which partly cancel the effect of bad conditioning.

Worked Example 1

Problem. For L(w) = 0.5*(w1^2 + 100 w2^2), find the condition number and the max stable learning rate.

  1. Hessian H = diag(1, 100); lambda_min = 1, lambda_max = 100.
  2. Condition number kappa = 100/1 = 100 (ill-conditioned).
  3. Max stable step: eta < 2/lambda_max = 2/100 = 0.02.
  4. Convergence factor at best eta: (kappa-1)/(kappa+1) = 99/101 ~ 0.98 per step (slow).

Answer. kappa = 100; eta must be < 0.02; convergence is slow at ~0.98 contraction per iteration.

Worked Example 2

Problem. Show how rescaling fixes conditioning: substitute u2 = 10 w2 in the loss above.

  1. Original stretches by 100 along w2; let u2 = 10 w2 so the 100 w2^2 term becomes u2^2.
  2. New loss 0.5*(w1^2 + u2^2) has Hessian I, condition number 1.
  3. Gradient descent now converges in essentially one step at eta=1.
  4. In ML this is exactly what standardizing features (zero mean, unit variance) achieves.

Answer. Rescaling drops kappa from 100 to 1, turning slow zigzag into near-instant convergence.

Common mistakes
  • Choosing one global learning rate for badly-scaled features. Unscaled features create large kappa; standardize before training.
  • Setting eta just below the divergence threshold for speed. Near 2/lambda_max the iterates oscillate; a moderate eta is more robust.
✎ Try it yourself

Problem. A quadratic loss has Hessian eigenvalues 4 and 36. Give the condition number and the largest stable learning rate for gradient descent.

Solution. Condition number kappa = lambda_max/lambda_min = 36/4 = 9 (moderately ill-conditioned). Largest stable step eta < 2/lambda_max = 2/36 ~ 0.0556. The per-iteration contraction at the best fixed step is (kappa-1)/(kappa+1) = 8/10 = 0.8, so each step removes about 20% of the error — slower than a well-conditioned problem but far better than kappa=100.

Stochastic gradient descent (SGD) and mini-batches

Full-data gradients are too expensive at scale; SGD trades exactness for many cheap, noisy steps.

Stochastic gradient descent estimates the gradient from a random mini-batch of examples instead of the whole dataset: w <- w - eta * grad L_batch(w). Each step is far cheaper (cost scales with batch size, not dataset size), enabling many more updates per epoch, which is why all large-scale ML uses it. The gradient is an unbiased but noisy estimate of the true gradient; the noise actually helps escape shallow saddle points and flat regions, but it also means a fixed step never fully settles — you typically decay eta over time. Mini-batches (e.g. 32-256) balance gradient-estimate variance against vectorization efficiency on hardware. SGD's variance also has a mild regularizing effect that can improve generalization.

Worked Example 1

Problem. Why is one SGD epoch with batch size 1 cheaper per step than batch GD on n=1,000,000 examples?

  1. Batch GD: each step processes all 1e6 examples -> 1 update per epoch.
  2. SGD (batch=1): each step processes 1 example -> 1e6 updates per epoch.
  3. Per-step cost ratio is 1e6 : 1, so SGD makes vastly more progress per pass.
  4. The trade is noisier individual gradients, mitigated by averaging over many steps.

Answer. SGD does ~1,000,000 cheap updates per epoch vs. one expensive batch-GD update, dramatically speeding early progress.

Worked Example 2

Problem. Implement a minimal SGD loop for linear regression with mini-batches in NumPy.

  1. `import numpy as np; X=...; y=...; n=len(y); w=np.zeros(X.shape[1]); eta=0.01; bs=32`
  2. `for epoch in range(20):\n idx=np.random.permutation(n)\n for s in range(0,n,bs):\n b=idx[s:s+bs]\n grad = (2/len(b)) * X[b].T @ (X[b]@w - y[b])\n w -= eta*grad`
  3. Shuffling each epoch decorrelates batches; the per-batch gradient approximates the full gradient.
  4. Decaying eta (e.g. eta/(1+0.01*epoch)) helps the noisy iterate settle.

Answer. The loop performs n/bs noisy updates per epoch; w converges near the least-squares solution with far less per-step cost than batch GD.

Common mistakes
  • Forgetting to shuffle data between epochs. Sorted or correlated batches bias the gradient estimates and slow or skew convergence.
  • Keeping a large fixed learning rate forever. SGD noise prevents settling; decay eta (or use a schedule) to converge tightly.
✎ Try it yourself

Problem. You have 10,000 examples and use mini-batches of size 50. How many SGD updates happen per epoch, and how does the gradient variance compare to batch size 500?

Solution. Updates per epoch = 10000 / 50 = 200. Batch size 500 would give 20 updates per epoch. The gradient-estimate variance scales like 1/batch_size, so the size-50 batches have about 10x the variance of size-500 batches (noisier per step) but yield 10x more updates per epoch. Smaller batches mean noisier but more frequent steps; larger batches mean smoother but fewer steps.

Momentum and Nesterov acceleration

Momentum gives the optimizer inertia, smoothing oscillation and accelerating through long valleys.

Momentum augments gradient descent with a velocity that accumulates past gradients: v <- mu*v - eta*grad; w <- w + v, with momentum coefficient mu (e.g. 0.9). On a stretched, ill-conditioned loss, plain gradient descent zigzags across the valley; momentum cancels the oscillating components and reinforces the consistent downhill direction, effectively reducing the dependence on the condition number. Nesterov accelerated gradient evaluates the gradient at the look-ahead point w + mu*v rather than at w, giving a more accurate correction and better theoretical convergence. Both are standard in deep learning because they speed training and help carry the iterate through small bumps and saddle regions where the raw gradient is tiny.

Worked Example 1

Problem. Two momentum steps on L(w)=(w-2)^2 from w=0, eta=0.1, mu=0.9, v0=0.

  1. grad = 2(w-2). Step 1: v = 0.9*0 - 0.1*2*(0-2) = 0.4; w = 0 + 0.4 = 0.4.
  2. Step 2: grad at 0.4 = 2*(0.4-2) = -3.2; v = 0.9*0.4 - 0.1*(-3.2) = 0.36 + 0.32 = 0.68; w = 0.4 + 0.68 = 1.08.
  3. Compare plain GD step 2 from 0.4: w = 0.4 - 0.1*(-3.2) = 0.72.
  4. Momentum (1.08) advances faster than plain GD (0.72) toward w*=2 thanks to accumulated velocity.

Answer. After two steps momentum reaches w ~ 1.08 vs. plain GD's 0.72; the velocity term accelerates progress.

Worked Example 2

Problem. State the Nesterov update and contrast it with classical momentum.

  1. Classical: v = mu*v - eta*grad(w); w = w + v.
  2. Nesterov: w_look = w + mu*v; v = mu*v - eta*grad(w_look); w = w + v.
  3. Nesterov measures the gradient where momentum is about to carry the iterate, so it 'corrects sooner'.
  4. `v = mu*v - eta*g(w + mu*v); w += v` is the one-line NumPy form.

Answer. Nesterov evaluates the gradient at the look-ahead point w+mu*v, yielding a more responsive, better-converging update than classical momentum.

Common mistakes
  • Setting momentum mu too high (e.g. 0.99) with a large eta. The iterate can overshoot and oscillate; lower eta when mu is large.
  • Forgetting to scale the effective step. With momentum the effective learning rate is roughly eta/(1-mu), so 0.9 momentum amplifies steps ~10x.
✎ Try it yourself

Problem. Run one momentum step on L(w)=(w-5)^2 from w=1 with eta=0.1, mu=0.9, v0=0. Then one more step.

Solution. grad = 2(w-5). Step 1: grad at 1 = 2*(1-5) = -8; v = 0.9*0 - 0.1*(-8) = 0.8; w = 1 + 0.8 = 1.8. Step 2: grad at 1.8 = 2*(1.8-5) = -6.4; v = 0.9*0.8 - 0.1*(-6.4) = 0.72 + 0.64 = 1.36; w = 1.8 + 1.36 = 3.16. The accumulated velocity makes the iterate climb toward w*=5 faster than plain GD would (which would give 1.8 then 2.44).

Adaptive methods: AdaGrad, RMSProp, Adam

These optimizers give each parameter its own learning rate, scaled by its gradient history — Adam is the deep-learning default.

Adaptive methods divide the step by an estimate of each parameter's gradient magnitude, so rarely-updated or steep coordinates get appropriately sized steps. AdaGrad accumulates squared gradients G += g^2 and steps eta*g/sqrt(G+eps); its denominator only grows, so the rate decays, sometimes too aggressively. RMSProp fixes this with an exponential moving average of squared gradients, keeping the scale responsive. Adam combines RMSProp's per-parameter scaling with momentum: it keeps an EMA of gradients m (first moment) and of squared gradients v (second moment), applies bias correction, and steps eta * m_hat / (sqrt(v_hat) + eps). Adam's robustness to learning-rate choice and bad conditioning makes it the default optimizer for training neural networks.

Worked Example 1

Problem. One Adam step for a single parameter with g=0.5, eta=0.1, beta1=0.9, beta2=0.999, eps=1e-8, t=1, m0=v0=0.

  1. m = 0.9*0 + 0.1*0.5 = 0.05; v = 0.999*0 + 0.001*0.25 = 0.00025.
  2. Bias correct: m_hat = 0.05/(1-0.9^1) = 0.05/0.1 = 0.5; v_hat = 0.00025/(1-0.999^1) = 0.00025/0.001 = 0.25.
  3. Update: w -= 0.1 * 0.5 / (sqrt(0.25) + 1e-8) = 0.1*0.5/0.5 = 0.1.
  4. On the first step Adam's effective step equals eta because bias correction recovers the raw gradient direction.

Answer. w decreases by 0.1 on the first step; bias correction makes the initial effective step ~eta independent of gradient scale.

Worked Example 2

Problem. Write the from-scratch Adam update in NumPy for a parameter vector.

  1. `m = beta1*m + (1-beta1)*g; v = beta2*v + (1-beta2)*g*g`
  2. `mhat = m/(1-beta1**t); vhat = v/(1-beta2**t)`
  3. `w -= eta * mhat / (np.sqrt(vhat) + eps)`
  4. Increment t each step; the per-coordinate sqrt(vhat) gives each weight its own adaptive rate.

Answer. Four lines implement Adam: EMA of g and g^2, bias-correct both, then take a per-parameter scaled step.

Common mistakes
  • Letting AdaGrad run too long. Its monotonically growing denominator can shrink the learning rate to near zero; RMSProp/Adam avoid this.
  • Omitting Adam's bias correction. Without dividing by (1 - beta^t), the early m and v are biased toward zero and steps are too small.
✎ Try it yourself

Problem. Compute the first RMSProp step for a parameter with g=2, eta=0.01, decay rho=0.9, eps=1e-8, accumulator s0=0.

Solution. Update accumulator: s = rho*s0 + (1-rho)*g^2 = 0.9*0 + 0.1*4 = 0.4. Step: delta = eta * g / (sqrt(s) + eps) = 0.01 * 2 / (sqrt(0.4) + 1e-8) = 0.02 / 0.6325 ~ 0.0316. So w decreases by about 0.0316. RMSProp normalized the raw gradient by its running RMS magnitude, giving a scale-robust step.

Learning-rate schedules and diagnosing training

A constant learning rate is rarely best; schedules and loss curves are how practitioners control and debug training.

A learning-rate schedule changes eta over training: step decay (drop by a factor every k epochs), exponential decay, cosine annealing, and warmup (ramp eta up early to stabilize Adam) are common. Decaying eta lets training take big early steps and fine, settling steps near the optimum, which is essential for SGD's noise to average out. Diagnosing training means reading the loss curve: a loss that explodes or oscillates signals eta too large; a loss that barely moves signals eta too small or vanishing gradients; a training loss that keeps falling while validation loss rises signals overfitting. Plotting loss and gradient norm versus iteration is the primary tool for tuning the rate and catching divergence early.

Worked Example 1

Problem. Step decay: eta0=0.1, halved every 10 epochs. What is eta at epochs 0, 10, 25?

  1. Schedule: eta = eta0 * 0.5^(floor(epoch/10)).
  2. Epoch 0: 0.1 * 0.5^0 = 0.1.
  3. Epoch 10: 0.1 * 0.5^1 = 0.05.
  4. Epoch 25: 0.1 * 0.5^2 = 0.025 (floor(25/10)=2).

Answer. eta = 0.1, 0.05, 0.025 at epochs 0, 10, 25 respectively.

Worked Example 2

Problem. A loss curve rises and oscillates wildly after a few steps. Diagnose and fix.

  1. Rising/oscillating loss is the signature of an overlarge learning rate overshooting the minimum.
  2. First fix: reduce eta (e.g. divide by 10) and re-run.
  3. Also check feature scaling (high condition number amplifies overshoot) and gradient correctness.
  4. `for lr in [1.0, 0.1, 0.01]: run(lr)` then keep the largest lr whose loss decreases monotonically.

Answer. Diagnosis: learning rate too large; remedy: lower eta (and standardize features), choosing the largest stable rate.

Worked Example 3

Problem. Cosine annealing: eta(t) = 0.5*eta_max*(1 + cos(pi*t/T)) with eta_max=0.2, T=100. Find eta at t=0, 50, 100.

  1. t=0: 0.5*0.2*(1+cos 0) = 0.1*(1+1) = 0.2.
  2. t=50: 0.5*0.2*(1+cos(pi/2)) = 0.1*(1+0) = 0.1.
  3. t=100: 0.5*0.2*(1+cos(pi)) = 0.1*(1-1) = 0.0.
  4. `import numpy as np; 0.5*0.2*(1+np.cos(np.pi*np.array([0,50,100])/100))` -> [0.2, 0.1, 0.0].

Answer. eta = 0.2, 0.1, 0.0 at t = 0, 50, 100; cosine annealing smoothly cools the rate to zero.

Common mistakes
  • Reading only the training loss. Overfitting hides there; always track validation loss to catch a rising generalization gap.
  • Decaying the learning rate too early or too fast. Premature decay stalls progress before the model reaches a good region.
✎ Try it yourself

Problem. Exponential decay eta(t) = eta0 * gamma^t with eta0=0.05, gamma=0.95. What is the learning rate at t=0, t=10, and t=20?

Solution. t=0: 0.05 * 0.95^0 = 0.05. t=10: 0.05 * 0.95^10 = 0.05 * 0.5987 ~ 0.0299. t=20: 0.05 * 0.95^20 = 0.05 * 0.3585 ~ 0.0179. The rate decays smoothly and multiplicatively, so big early steps shrink to small refining steps as training proceeds.

Implementing optimizers from scratch

Coding GD, SGD, momentum, and Adam yourself cements exactly what each update rule does — the capstone of this unit.

Implementing optimizers from scratch means writing the parameter-update rule directly given a function that returns the loss gradient. The shared skeleton is a loop that computes grad = grad_loss(w) and applies an update; the optimizer differs only in how grad becomes a step. Plain GD: w -= eta*grad. Momentum: maintain v. Adam: maintain m, v, and a step counter t. Building these makes the mechanics concrete and lets you verify against a known optimum, test learning-rate sensitivity, and confirm momentum/Adam beat plain GD on ill-conditioned losses. This is also how you debug training in practice: a hand-written optimizer that matches a library's behavior proves you understand the math driving model fitting.

Worked Example 1

Problem. Write a generic gradient-descent optimizer and run it on a quadratic to verify the minimum.

  1. `import numpy as np`
  2. `def gd(grad, w, eta=0.1, steps=100):\n for _ in range(steps): w = w - eta*grad(w)\n return w`
  3. `grad = lambda w: 2*(w - np.array([3.0, -1.0]))`
  4. `gd(grad, np.zeros(2))` -> ~[3, -1], the known minimizer of ||w - [3,-1]||^2.

Answer. The from-scratch GD converges to [3, -1], matching the analytic minimizer.

Worked Example 2

Problem. Add momentum and Adam variants that share the same gradient interface.

  1. Momentum: `v = mu*v - eta*grad(w); w = w + v` inside the loop (init v = zeros).
  2. Adam: track `m, v, t`; `m=b1*m+(1-b1)*g; v=b2*v+(1-b2)*g*g; t+=1; mh=m/(1-b1**t); vh=v/(1-b2**t); w-=eta*mh/(np.sqrt(vh)+1e-8)`.
  3. Run all three from the same w0 on an ill-conditioned quadratic (Hessian diag(1,100)).
  4. Track ||grad|| each step; Adam and momentum reach a small gradient norm in far fewer steps than plain GD.

Answer. Sharing the gradient interface, momentum and Adam converge on the diag(1,100) bowl in many fewer iterations than plain GD, demonstrating their conditioning advantage.

Common mistakes
  • Reusing stale optimizer state across separate runs. Reset velocity v, moments m/v, and the step counter t before each new training run.
  • Mismatching the gradient sign. The update subtracts the gradient (w -= eta*grad); a sign flip turns descent into ascent and diverges.
✎ Try it yourself

Problem. Write pseudocode for a from-scratch momentum optimizer and apply it to minimize L(w)=(w-4)^2 from w=0, eta=0.1, mu=0.9. Give the first two w values.

Solution. Pseudocode: v=0; repeat { g = 2*(w-4); v = 0.9*v - 0.1*g; w = w + v }. Step 1: g=2*(0-4)=-8; v=0.9*0 - 0.1*(-8)=0.8; w=0+0.8=0.8. Step 2: g=2*(0.8-4)=-6.4; v=0.9*0.8 - 0.1*(-6.4)=0.72+0.64=1.36; w=0.8+1.36=2.16. The iterate accelerates toward w*=4 via accumulated velocity.

Key terms
  • Gradient descent — the iterative update x <- x - eta * gradient, moving downhill against the gradient.
  • Learning rate (eta) — the step-size multiplier; too large diverges, too small crawls.
  • Stochastic gradient descent (SGD) — uses a random sample (mini-batch) to estimate the gradient each step.
  • Mini-batch — a small subset of data used to compute a noisy but cheap gradient estimate.
  • Momentum — accumulates a velocity from past gradients to smooth and accelerate descent through valleys.
  • Nesterov acceleration — a look-ahead momentum variant that evaluates the gradient at the predicted next point.
  • Adam — an adaptive optimizer combining momentum with per-parameter scaling from squared-gradient estimates.
  • Condition number — the ratio of largest to smallest Hessian eigenvalue; high values slow plain gradient descent.
Assignment · Build and race your own optimizers

On a logistic-regression or least-squares loss with a known minimum, implement plain gradient descent, SGD with mini-batches, momentum, and Adam from scratch in Python (NumPy only). Run each from the same initialization and plot loss vs. iteration and vs. wall-clock. Sweep the learning rate to show one value that diverges and one that converges. Explain why momentum and Adam handle an ill-conditioned (stretched) loss better than plain GD.

Deliverable · A NumPy notebook with from-scratch optimizers, convergence plots comparing all four, and a short write-up on learning-rate sensitivity and conditioning.

Quiz · 5 questions
  1. 1. In gradient descent x <- x - eta * grad, a learning rate that is too large typically causes:

  2. 2. What is the main difference between batch gradient descent and SGD?

  3. 3. Momentum improves gradient descent primarily by:

  4. 4. Adam is best described as combining:

  5. 5. A high condition number of the loss Hessian causes plain gradient descent to:

You'll be able to

I can implement batch gradient descent and explain how the learning rate affects convergence.

I can describe how SGD, momentum, and Adam modify the basic update rule and when each helps.

I can diagnose divergence or slow training and adjust the step size or schedule.

Weeks 7-8 Unit 4: Constrained Optimization (Lagrange, KKT)
apply-lagrange-multipliersinterpret-multiplier-geometryhandle-inequality-constraintsstate-kkt-conditionsapply-complementary-slacknesssolve-constrained-in-code
Lecture
Equality constraints and the method of Lagrange multipliers

When you must optimize subject to an equation, Lagrange multipliers turn it back into solving for stationary points.

To minimize f(x) subject to an equality constraint h(x) = 0, form the Lagrangian L(x, lambda) = f(x) + lambda*h(x). At a constrained optimum the gradient of L with respect to x vanishes, giving grad f(x) = -lambda*grad h(x): the objective gradient is parallel to the constraint gradient. Solving grad_x L = 0 together with h(x) = 0 yields the candidate point and the multiplier lambda. The multiplier is not just bookkeeping — it measures the sensitivity of the optimal value to the constraint level (a shadow price). This method generalizes to several equalities (one lambda each) and is the foundation for KKT with inequalities. In ML it underlies constrained estimators and equality-constrained least squares.

Worked Example 1

Problem. Minimize f(x,y) = x^2 + y^2 subject to x + y = 1.

  1. L = x^2 + y^2 + lambda(x + y - 1).
  2. grad_x L: 2x + lambda = 0, 2y + lambda = 0 -> x = y = -lambda/2.
  3. Constraint: x + y = 1 -> 2x = 1 -> x = y = 0.5, lambda = -1.
  4. `import numpy as np; from scipy.optimize import minimize; minimize(lambda v:v@v, [0,0], constraints={'type':'eq','fun':lambda v:v[0]+v[1]-1}).x` -> [0.5, 0.5].

Answer. Minimizer (0.5, 0.5) with value 0.5; multiplier lambda = -1.

Worked Example 2

Problem. Find the rectangle of maximum area with perimeter 20 (maximize xy s.t. 2x + 2y = 20).

  1. Maximize xy <=> minimize -xy; constraint x + y = 10.
  2. L = -xy + lambda(x + y - 10); grad: -y + lambda = 0, -x + lambda = 0 -> x = y = lambda.
  3. Constraint: x + y = 10 -> x = y = 5.
  4. Area = 25; the square is optimal, a classic Lagrange result.

Answer. x = y = 5 (a 5x5 square), maximum area 25; lambda = 5.

Common mistakes
  • Forgetting the constraint equation when solving. You need both grad_x L = 0 AND h(x)=0 to pin down x and lambda.
  • Treating the multiplier sign as meaningless. Its sign and magnitude carry the sensitivity (shadow-price) information.
✎ Try it yourself

Problem. Minimize f(x,y) = x^2 + 4y^2 subject to x + 2y = 8 using Lagrange multipliers.

Solution. L = x^2 + 4y^2 + lambda(x + 2y - 8). grad_x: 2x + lambda = 0 -> x = -lambda/2. grad_y: 8y + 2lambda = 0 -> y = -lambda/4. Constraint: x + 2y = 8 -> -lambda/2 - lambda/2 = -lambda = 8 -> lambda = -8. Then x = 4, y = 2. Check: 4 + 2*2 = 8. Minimum value = 16 + 16 = 32, attained at (4, 2).

Geometric intuition: gradients align at the optimum

The picture behind Lagrange: at the best feasible point, you cannot improve without violating the constraint.

Geometrically, the constraint h(x) = 0 defines a surface (a curve in 2D), and grad h is perpendicular to it. The objective's level sets f(x) = c are nested curves; moving along the constraint surface changes f unless the surface is tangent to a level set. At the constrained optimum the level set is tangent to the constraint, which means grad f is perpendicular to the surface too — hence parallel to grad h. That is exactly grad f = -lambda*grad h. If the gradients were not aligned, there would be a direction along the surface that decreases f, so you could improve. This tangency picture is the intuition that makes the Lagrange equations memorable and explains why the multiplier scales the two parallel gradients.

Worked Example 1

Problem. On x^2 + y^2 = 1, where is f(x,y) = x maximized? Use gradient alignment.

  1. grad f = (1, 0); constraint grad h = (2x, 2y).
  2. Alignment: (1,0) = lambda(2x, 2y) -> y = 0, 1 = 2 lambda x.
  3. On the circle with y=0: x = +/-1; x=1 maximizes f.
  4. At (1,0) the level set x=1 is tangent to the circle — gradients are parallel.

Answer. Maximum of f=x on the unit circle is at (1,0), where grad f and grad h align.

Worked Example 2

Problem. Why does a feasible point with non-parallel grad f and grad h fail to be optimal?

  1. Decompose grad f into a part along grad h and a part tangent to the constraint surface.
  2. If the tangential part is nonzero, moving along the surface (staying feasible) changes f.
  3. You can move in the direction that decreases f, so the point is not a minimum.
  4. Only when the tangential component is zero (gradients parallel) is no feasible improvement possible.

Answer. A nonzero tangential component of grad f means a feasible descent direction exists, so the point cannot be optimal; optimality forces gradient alignment.

Common mistakes
  • Thinking grad f must be zero at a constrained optimum. It need not be; it only must be parallel to the constraint gradient.
  • Ignoring the constraint-gradient regularity. If grad h = 0 at the point, the alignment argument can break down (a constraint-qualification issue).
✎ Try it yourself

Problem. Maximize f(x,y) = 3x + 4y on the circle x^2 + y^2 = 25 using gradient alignment, and interpret geometrically.

Solution. grad f = (3,4), grad h = (2x, 2y). Alignment: (3,4) = lambda(2x,2y) -> x = 3/(2lambda), y = 4/(2lambda) = 2/lambda. On the circle: x^2+y^2 = (9 + 16)/(4lambda^2) = 25/(4lambda^2) = 25 -> lambda^2 = 1/4 -> lambda = 1/2 (take positive for max). Then x = 3, y = 4. The level line 3x+4y is tangent to the circle at (3,4); the maximum value is 3*3+4*4 = 25.

Inequality constraints and feasibility

Real problems cap and floor things with inequalities; here we handle <= constraints and the idea of active vs. slack.

An inequality constraint g(x) <= 0 either binds at the optimum (g(x*) = 0, active) or has slack (g(x*) < 0, inactive). If a constraint is inactive, the optimum is the same as if it were absent — its multiplier is zero. If it is active, it behaves like an equality and gets a non-negative multiplier. The first task is feasibility: the feasible set is the intersection of all g_i <= 0 and h_j = 0 regions, and you must check candidates lie inside it. Reasoning about which inequalities are active is the key step in solving constrained problems by hand; you guess an active set, solve the resulting equality-constrained system, and verify the multipliers are non-negative and the inactive constraints hold. This active-set logic powers QP solvers used throughout ML.

Worked Example 1

Problem. Minimize f(x) = (x-2)^2 subject to x <= 5. Is the constraint active?

  1. Unconstrained minimizer is x = 2.
  2. Check feasibility: 2 <= 5 holds with slack.
  3. Since the unconstrained optimum is feasible, the constraint is inactive.
  4. Its multiplier is 0; the answer is x* = 2.

Answer. x* = 2; the constraint x<=5 is inactive (slack), multiplier 0.

Worked Example 2

Problem. Minimize f(x) = (x-2)^2 subject to x <= 1. Now is the constraint active?

  1. Unconstrained minimizer x=2 violates x<=1, so it is infeasible.
  2. The objective decreases toward x=2, so on the feasible set [.., 1] the best is the boundary x=1.
  3. Constraint is active: g(x*)=x-1=0.
  4. Multiplier: grad f + mu*grad g = 2(x-1=1? ) -> 2(1-2) + mu*1 = -2 + mu = 0 -> mu = 2 >= 0, valid.

Answer. x* = 1; the constraint is active with multiplier mu = 2 (non-negative, as required).

Common mistakes
  • Assuming every constraint is active. Inactive constraints have zero multiplier and can be dropped at the optimum.
  • Accepting a candidate without checking the multiplier sign. For g<=0, an active constraint must have mu >= 0; a negative mu signals the constraint should be inactive.
✎ Try it yourself

Problem. Minimize f(x) = (x-4)^2 subject to x <= 6 and x >= 0. Which constraints are active, and what is x*?

Solution. Unconstrained minimizer is x=4. Check feasibility: 0 <= 4 <= 6 holds with slack on both sides. So both constraints are inactive (multipliers 0) and x* = 4 with f=0. Neither bound binds because the free optimum already lies inside the feasible interval [0,6].

The Karush-Kuhn-Tucker (KKT) conditions

KKT packages everything — equalities, inequalities, and optimality — into four conditions that characterize a solution.

For minimize f(x) s.t. g_i(x) <= 0 and h_j(x) = 0, the KKT conditions are the first-order necessary conditions (and, for convex problems with a constraint qualification, sufficient). They are: (1) Stationarity — grad f + sum mu_i grad g_i + sum lambda_j grad h_j = 0; (2) Primal feasibility — g_i(x) <= 0 and h_j(x) = 0; (3) Dual feasibility — mu_i >= 0; (4) Complementary slackness — mu_i*g_i(x) = 0. A point satisfying all four (with a qualification like Slater's) is optimal for a convex problem. KKT generalizes both the unconstrained gradient=0 condition and Lagrange's equality condition, and it is the certificate that solvers like CVXPY verify, returning the multipliers as dual values.

Worked Example 1

Problem. Write and check the KKT conditions for minimize x^2 s.t. x >= 1 (i.e. 1 - x <= 0).

  1. g(x) = 1 - x <= 0; Lagrangian stationarity: 2x + mu*(-1) = 0 -> mu = 2x.
  2. Guess active: 1 - x = 0 -> x = 1, then mu = 2 >= 0 (dual feasible).
  3. Complementary slackness: mu*(1-x) = 2*0 = 0, satisfied.
  4. Primal feasibility x=1>=1 holds. All four KKT conditions met.

Answer. x* = 1, mu = 2; all four KKT conditions hold, so x=1 is optimal.

Worked Example 2

Problem. Verify KKT for minimize (x-3)^2 s.t. x <= 5 where the constraint is inactive.

  1. g(x) = x - 5 <= 0; stationarity 2(x-3) + mu = 0.
  2. Guess inactive: mu = 0 -> 2(x-3) = 0 -> x = 3.
  3. Primal feasibility: 3 - 5 = -2 <= 0, holds with slack.
  4. Complementary slackness mu*g = 0*(-2) = 0, and dual feasibility mu=0>=0 hold.

Answer. x* = 3, mu = 0; KKT confirmed with the constraint inactive.

Worked Example 3

Problem. Set up KKT for minimize x^2 + y^2 s.t. x + y = 1 and x >= 0.

  1. Constraints: h = x+y-1 = 0, g = -x <= 0. Stationarity: (2x, 2y) + lambda(1,1) + mu(-1, 0) = 0.
  2. So 2x + lambda - mu = 0 and 2y + lambda = 0.
  3. Guess x>=0 inactive (mu=0): then 2x+lambda=0, 2y+lambda=0 -> x=y=0.5 from h, lambda=-1.
  4. Check x=0.5 >= 0 (slack), mu=0>=0, complementary slackness 0*(-0.5)=0 — all KKT hold.

Answer. x*=y*=0.5, lambda=-1, mu=0; the inequality is inactive and all KKT conditions are satisfied.

Common mistakes
  • Dropping dual feasibility. A solution needs mu_i >= 0 for g_i <= 0 constraints; a negative multiplier invalidates the KKT point.
  • Treating KKT as sufficient for nonconvex problems. Without convexity (and a constraint qualification), KKT is only necessary.
✎ Try it yourself

Problem. State and verify the KKT conditions for minimize (x-5)^2 subject to x <= 2.

Solution. g(x) = x - 2 <= 0. Stationarity: 2(x-5) + mu = 0. The free optimum x=5 violates x<=2, so guess active: x=2, then mu = -2(2-5) = 6 >= 0 (dual feasible). Complementary slackness: mu*(x-2) = 6*0 = 0, holds. Primal feasibility x=2<=2 holds. All four KKT conditions are satisfied, so x* = 2 with multiplier mu = 6.

Complementary slackness and active constraints

This single condition links which constraints bind to which multipliers are nonzero — the heart of solving KKT by hand.

Complementary slackness states mu_i * g_i(x*) = 0 for every inequality. It forces an either/or: either the constraint is active (g_i = 0) so its multiplier may be positive, or the constraint is inactive (g_i < 0) so its multiplier must be zero. They cannot both be nonzero. This is what makes the active-set strategy work: you hypothesize a set of active constraints, set their slacks to zero and the inactive multipliers to zero, solve the reduced equality-constrained system, then verify dual feasibility (mu >= 0) and primal feasibility of the dropped constraints. The nonzero multipliers identify the binding constraints and quantify their shadow prices — how much the optimum would improve if you relaxed each binding limit.

Worked Example 1

Problem. A problem has constraints g1: x <= 4 and g2: y <= 10. At the optimum x=4, y=3. Which multipliers can be nonzero?

  1. g1 active: 4 - 4 = 0, so mu1 may be > 0.
  2. g2 slack: 3 - 10 = -7 < 0, so by complementary slackness mu2 = 0.
  3. Only the binding constraint carries a nonzero multiplier.
  4. Interpretation: relaxing x<=4 (raising the cap) would change the optimum; relaxing y<=10 would not.

Answer. mu1 can be positive (g1 active); mu2 must be 0 (g2 has slack).

Worked Example 2

Problem. Use the active-set idea: minimize x^2 + y^2 s.t. x + y >= 4 (i.e. 4 - x - y <= 0). Solve.

  1. Free optimum (0,0) violates the constraint, so guess active: x + y = 4.
  2. By symmetry/Lagrange on the equality: x = y = 2.
  3. Stationarity: (2x,2y) + mu*(-1,-1)=0 -> mu = 2x = 4 >= 0, dual feasible.
  4. Complementary slackness mu*(4-x-y) = 4*0 = 0 holds; solution (2,2).

Answer. x* = y* = 2, mu = 4 (constraint active); complementary slackness confirmed.

Common mistakes
  • Allowing a positive multiplier on a slack constraint. If g_i < 0 strictly, mu_i must be exactly 0; otherwise complementary slackness is violated.
  • Guessing the wrong active set and not re-checking. If a multiplier comes out negative, that constraint should be inactive — revise the active set.
✎ Try it yourself

Problem. Minimize x^2 + y^2 subject to x >= 3 (i.e. 3 - x <= 0) and y >= 0. Use complementary slackness to find the optimum and active set.

Solution. Free optimum is (0,0), which violates x>=3, so guess g1: 3-x active. Then x=3. For y: the y>=0 constraint — the free optimum y=0 is feasible and minimizes y^2, so y=0 with that constraint exactly at its boundary but the unconstrained y optimum is already 0, multiplier 0. Stationarity in x: 2x + mu1*(-1)=0 -> mu1 = 2*3 = 6 >= 0. Complementary slackness: mu1*(3-x)=6*0=0. So x*=3, y*=0, mu1=6. Only x>=3 is the binding active constraint.

Worked examples and solving with SciPy/CVXPY

Bring KKT to code: solvers compute the same primal solution and return the multipliers as dual values you can interpret.

Constrained problems are solved in practice with scipy.optimize.minimize (using a constraints list of dicts for 'eq'/'ineq') or, for convex problems, CVXPY, which is more reliable and exposes dual values directly. In SciPy, inequality constraints are written as fun >= 0; in CVXPY you write natural relations like x >= 0 and read constraint.dual_value. The dual values are exactly the KKT multipliers, so you can confirm hand-derived lambdas/mus and interpret them as shadow prices. The practical workflow is: model the objective and constraints, solve, check the solver converged and the solution is feasible, then read the duals to understand which constraints bind and how sensitive the optimum is. This closes the loop between KKT theory and reliable computation.

Worked Example 1

Problem. Solve minimize x^2 + y^2 s.t. x + y = 1 in CVXPY and read the dual.

  1. `import cvxpy as cp; x = cp.Variable(2)`
  2. `con = [cp.sum(x) == 1]; prob = cp.Problem(cp.Minimize(cp.sum_squares(x)), con)`
  3. `prob.solve()` -> x.value ~ [0.5, 0.5], prob.value ~ 0.5.
  4. `con[0].dual_value` ~ -1, matching the hand-derived multiplier lambda = -1.

Answer. CVXPY returns x=[0.5,0.5], value 0.5, dual ~ -1, agreeing with the Lagrange solution.

Worked Example 2

Problem. Solve minimize (x-2)^2 s.t. x <= 1 with scipy.optimize.minimize and verify the multiplier.

  1. `from scipy.optimize import minimize`
  2. `con = {'type':'ineq', 'fun': lambda v: 1 - v[0]}` # encodes x <= 1 as 1 - x >= 0
  3. `res = minimize(lambda v:(v[0]-2)**2, [0.0], constraints=[con])`
  4. `res.x` -> ~[1.0]; the active constraint gives KKT multiplier mu = 2 (since 2(x-2)+mu*(-1)=0 at x=1).

Answer. x* = 1 (constraint active); multiplier mu = 2, matching the hand KKT derivation.

Common mistakes
  • Mixing up sign conventions between SciPy and CVXPY. SciPy 'ineq' means fun >= 0; CVXPY uses natural <=/>= relations — write constraints accordingly.
  • Reading dual values without confirming convergence and feasibility. Check res.success / prob.status before trusting the multipliers.
✎ Try it yourself

Problem. Set up CVXPY code to minimize x^2 + y^2 subject to x + y >= 4, and predict the solution and dual sign.

Solution. import cvxpy as cp; x = cp.Variable(2); con = [cp.sum(x) >= 4]; prob = cp.Problem(cp.Minimize(cp.sum_squares(x)), con); prob.solve(). The optimum is x = [2, 2] (value 8) since the constraint binds at x+y=4 and symmetry gives equal coordinates. The dual value of the binding x+y>=4 constraint is positive (about 4), reflecting that tightening the requirement raises the optimal cost — a positive shadow price.

Key terms
  • Lagrange multiplier — a scalar lambda introduced per equality constraint so that gradients balance at the optimum.
  • Lagrangian — L(x, lambda, mu) = f(x) + sum lambda_j h_j(x) + sum mu_i g_i(x), combining objective and constraints.
  • Stationarity (KKT) — the gradient of the Lagrangian with respect to x is zero at an optimum.
  • Primal feasibility — the candidate point satisfies all the original equality and inequality constraints.
  • Dual feasibility — the inequality multipliers mu_i are non-negative.
  • Complementary slackness — for each inequality, either the constraint is active or its multiplier is zero (mu_i g_i = 0).
  • Active constraint — an inequality constraint that holds with equality at the optimum, binding the solution.
  • Shadow price — the multiplier's interpretation as the marginal change in optimal value per unit relaxation of a constraint.
Assignment · From Lagrange to KKT by hand and in code

Choose a small problem such as minimizing a quadratic subject to one linear equality and one inequality constraint. Solve it by hand with the Lagrangian and verify all four KKT conditions, identifying which constraint is active. Then reproduce the solution in Python with both scipy.optimize.minimize (using constraints) and CVXPY, and confirm the solver's dual values match your multipliers. Interpret one multiplier as a shadow price.

Deliverable · A notebook showing the hand KKT derivation, the SciPy and CVXPY solutions, and a paragraph interpreting the multipliers and the active constraint.

Quiz · 5 questions
  1. 1. At the optimum of an equality-constrained problem, the gradient of the objective is:

  2. 2. Which is NOT one of the KKT conditions?

  3. 3. Complementary slackness states that for each inequality constraint:

  4. 4. For inequality constraints written as g_i(x) <= 0, dual feasibility requires the multipliers mu_i to be:

  5. 5. An inequality constraint that holds with equality at the optimum is called:

You'll be able to

I can solve equality-constrained problems with Lagrange multipliers and interpret the multipliers.

I can write and check the KKT conditions for a problem with inequality constraints.

I can explain complementary slackness and identify which constraints are active at the optimum.

Weeks 9-10 Unit 5: Linear & Quadratic Programming
formulate-lpexplain-lp-algorithmsapply-lp-dualityformulate-qpmodel-as-lp-qpsolve-lp-qp-in-code
Lecture
Linear programs: standard form and geometry of polytopes

Linear programs optimize a linear objective over a region carved by linear constraints; their geometry explains why the optimum lives at a corner.

A linear program (LP) minimizes a linear objective c^T x subject to linear constraints A x <= b, A_eq x = b_eq, and bounds. The feasible region is a polytope: an intersection of half-spaces with flat faces and sharp corners (vertices). Because both the objective contours and the constraints are linear (flat), the objective can only be pushed lower until it hits the boundary; the lowest point is always attained at a vertex (or along an edge connecting two equally optimal vertices). This corner principle is the foundation of LP algorithms and matters in ML for resource allocation, L1-norm fitting recast as LPs, and optimal transport. Recognizing an LP and writing it in standard form (minimize c^T x s.t. linear constraints) lets you hand it to a guaranteed-optimal solver.

Worked Example 1

Problem. Put 'maximize 3x + 2y subject to x + y <= 4, x <= 3, x,y >= 0' into LP standard form and identify the polytope vertices.

  1. Negate the objective for minimization: minimize -3x - 2y.
  2. Constraints stay linear: x + y <= 4, x <= 3, x >= 0, y >= 0.
  3. Vertices are corners where constraints intersect: (0,0), (3,0), (3,1), (0,4).
  4. Evaluate -3x-2y at each: (0,0)->0, (3,0)->-9, (3,1)->-11, (0,4)->-8; the minimum -11 is at vertex (3,1).

Answer. Standard form: minimize -3x-2y s.t. x+y<=4, x<=3, x,y>=0; optimum at vertex (3,1) with objective value -11 (original max 11).

Worked Example 2

Problem. Solve the LP from Example 1 with SciPy and confirm the corner solution.

  1. `from scipy.optimize import linprog`
  2. `c = [-3, -2]; A_ub = [[1,1],[1,0]]; b_ub = [4, 3]; bnds = [(0,None),(0,None)]`
  3. `res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bnds)`
  4. `res.x` -> [3., 1.]; `res.fun` -> -11.0, matching the vertex enumeration.

Answer. linprog returns x=[3,1], fun=-11, confirming the optimum sits at the polytope vertex (3,1).

Worked Example 3

Problem. Explain why an LP can have infinitely many optimal solutions, and give a condition.

  1. If the objective contour c^T x is parallel to an edge (face) of the polytope, every point on that face is equally optimal.
  2. Example: minimize -x - y s.t. x + y <= 1, x,y >= 0; the contour x+y is parallel to the binding edge x+y=1.
  3. Both vertices (1,0) and (0,1) and the entire segment between them give objective -1.
  4. `linprog([-1,-1], A_ub=[[1,1]], b_ub=[1], bounds=[(0,None)]*2).fun` -> -1.0 with x at a vertex, but the whole edge is optimal.

Answer. When the objective gradient is parallel to a constraint face, the optimal set is that entire face (an edge or facet), giving infinitely many optima with the same value.

Common mistakes
  • Forgetting to negate the objective for a maximization. linprog minimizes, so max c^T x becomes min -c^T x and the reported value flips sign.
  • Assuming the optimum is interior. For an LP the optimum is always on the boundary, at a vertex (or face); never search the interior for an LP optimum.
✎ Try it yourself

Problem. Write 'maximize 5x + 4y subject to 6x + 4y <= 24, x + 2y <= 6, x,y >= 0' in LP standard form, then solve it with scipy.optimize.linprog.

Solution. Standard form: minimize -5x - 4y s.t. 6x+4y<=24, x+2y<=6, x,y>=0. In code: from scipy.optimize import linprog; res = linprog([-5,-4], A_ub=[[6,4],[1,2]], b_ub=[24,6], bounds=[(0,None),(0,None)]). The optimum is the vertex where 6x+4y=24 and x+2y=6 meet: solving gives x=3, y=1.5. res.x ~ [3, 1.5], res.fun ~ -21, so the maximum profit is 21 at the corner (3, 1.5).

The simplex method and interior-point methods (intuition)

Two families of algorithms solve LPs: simplex walks the corners, interior-point cuts through the middle.

The simplex method exploits the corner principle: it starts at a feasible vertex and repeatedly moves along an edge to an adjacent vertex that improves the objective, stopping when no neighbor is better — that vertex is globally optimal. Each pivot swaps which constraints are 'active'. Simplex is fast in practice though worst-case exponential. Interior-point methods take a different route: they stay strictly inside the feasible region and follow a smooth 'central path' toward the optimum, solving a sequence of barrier-penalized problems with Newton steps. Interior-point methods have polynomial-time guarantees and scale well to large, sparse LPs and to QPs/conic programs. SciPy's linprog defaults to an interior-point-style 'highs' solver. Knowing both pictures explains solver options and why one may outperform the other on a given problem.

Worked Example 1

Problem. Trace simplex on: maximize 3x + 2y s.t. x + y <= 4, x <= 3, x,y >= 0, starting at vertex (0,0).

  1. At (0,0) objective = 0; both x and y have positive objective coefficients, so increasing either improves.
  2. Pivot along the x-edge: increase x until x<=3 binds, reaching (3,0), objective 9.
  3. From (3,0), increase y until x+y<=4 binds, reaching (3,1), objective 11.
  4. At (3,1) no adjacent vertex improves the objective -> optimal.

Answer. Simplex path (0,0)->(3,0)->(3,1); optimum 11 at (3,1), reached in two pivots.

Worked Example 2

Problem. Solve the same LP with SciPy's interior-point/HiGHS solver and contrast the route.

  1. `from scipy.optimize import linprog`
  2. `res = linprog([-3,-2], A_ub=[[1,1],[1,0]], b_ub=[4,3], bounds=[(0,None)]*2, method='highs-ipm')`
  3. `res.x` -> [3., 1.], `res.fun` -> -11.0 (same optimum).
  4. Interior-point reaches (3,1) by following the central path from inside, not by visiting (0,0) and (3,0).

Answer. Interior-point returns the same optimum x=[3,1], value -11, but approaches it through the feasible interior rather than corner-to-corner.

Common mistakes
  • Believing simplex examines all vertices. It follows an improving path and typically visits only a few, ignoring most corners.
  • Assuming interior-point lands exactly on a vertex. It converges toward the optimal face; a cleanup/crossover step is used to recover an exact vertex when one is needed.
✎ Try it yourself

Problem. For maximize x + y s.t. x <= 2, y <= 3, x,y >= 0, list the simplex vertex path from (0,0) and give the optimum.

Solution. Vertices: (0,0), (2,0), (2,3), (0,3). From (0,0), increase x to its bound -> (2,0), objective 2; then increase y to its bound -> (2,3), objective 5; no neighbor improves, so (2,3) is optimal. Path: (0,0)->(2,0)->(2,3). Optimum value 5 at (2,3). (linprog([-1,-1], bounds=[(0,2),(0,3)]).fun -> -5.)

LP duality and sensitivity analysis

Every LP has a twin dual LP; its solution reveals shadow prices that tell you the value of relaxing each constraint.

Every LP (the primal) has an associated dual LP built from the same data: a maximization with one dual variable per primal constraint. Strong duality for LPs guarantees their optimal values are equal, so the dual gives a certificate of optimality and a lower bound on the primal at every step. The dual variables are shadow prices: each equals the marginal change in the optimal objective per unit increase of the corresponding constraint's right-hand side. Sensitivity analysis reads these prices to answer 'which resource is worth buying more of?' without re-solving. In ML and operations, duality underlies SVMs (the dual form), optimal transport, and resource pricing. Recognizing the primal-dual pair and interpreting multipliers is what turns an LP solution into actionable insight.

Worked Example 1

Problem. Write the dual of: minimize 4x1 + 3x2 s.t. x1 + x2 >= 2, x1 + 3x2 >= 3, x >= 0.

  1. Primal is a minimization with >= constraints; assign dual variables y1, y2 >= 0.
  2. Dual is a maximization: maximize 2 y1 + 3 y2 (using the RHS as objective coefficients).
  3. Dual constraints come from primal columns: y1 + y2 <= 4 (for x1), y1 + 3 y2 <= 3 (for x2).
  4. So the dual is: maximize 2y1 + 3y2 s.t. y1 + y2 <= 4, y1 + 3y2 <= 3, y >= 0.

Answer. Dual: maximize 2y1 + 3y2 s.t. y1 + y2 <= 4, y1 + 3y2 <= 3, y1,y2 >= 0; its optimal value equals the primal's by strong duality.

Worked Example 2

Problem. Read shadow prices from a solved LP using SciPy's marginals.

  1. `from scipy.optimize import linprog`
  2. `res = linprog([4,3], A_ub=[[-1,-1],[-1,-3]], b_ub=[-2,-3], bounds=[(0,None)]*2)`
  3. `res.fun` gives the optimal cost; `res.ineqlin.marginals` gives the per-constraint shadow prices.
  4. A nonzero marginal on a constraint means tightening/relaxing its RHS changes the optimal cost at that marginal rate.

Answer. The constraint marginals returned by HiGHS are the dual values (shadow prices): each tells how much the optimal cost moves per unit change in that constraint's bound.

Worked Example 3

Problem. A factory's profit LP gives a shadow price of 5 on the labor constraint (100 hours). What does buying 10 more hours do?

  1. Shadow price = marginal change in optimal objective per unit of the resource = 5 per labor-hour.
  2. Adding 10 hours (within the validity range) raises optimal profit by about 5 * 10 = 50.
  3. If extra hours cost less than 5 each, buying them is profitable; if more than 5, it is not.
  4. Validity holds only while the same constraints stay active (the sensitivity range).

Answer. 10 extra labor-hours increase optimal profit by ~50 (price 5 each), so it pays to buy them only if the cost per hour is below 5, and only within the range where the active set is unchanged.

Common mistakes
  • Confusing the sign/direction of shadow prices. For a >= resource constraint, relaxing means lowering the bound; track which direction increases feasibility.
  • Extrapolating a shadow price arbitrarily far. The marginal value is valid only over the sensitivity range where the optimal basis (active set) does not change.
✎ Try it yourself

Problem. State the dual of: maximize 3x1 + 5x2 s.t. x1 <= 4, 2x2 <= 12, 3x1 + 2x2 <= 18, x >= 0, and explain what the dual variable on the third constraint means.

Solution. Assign duals y1,y2,y3 >= 0 to the three <= constraints. The dual is: minimize 4 y1 + 12 y2 + 18 y3 s.t. y1 + 3 y3 >= 3 (column x1), 2 y2 + 2 y3 >= 5 (column x2), y >= 0. By strong duality the dual optimum equals the primal optimum (36 for this classic problem). The dual variable y3 on 3x1+2x2<=18 is the shadow price of that resource: it equals the increase in maximum profit per extra unit of the 18-unit budget, valid while the same constraints remain active.

Quadratic programs and their structure

Adding a quadratic objective to linear constraints gives the QP, the workhorse behind SVMs, portfolio optimization, and least squares.

A quadratic program (QP) minimizes (1/2) x^T P x + q^T x subject to linear constraints A x <= b, A_eq x = b_eq. The objective is quadratic, the feasible set still a polytope. The QP is convex exactly when the matrix P is positive semidefinite; then it has a unique optimal value and efficient solvers reach the global optimum. Unlike an LP, the optimum need not be at a vertex — the quadratic bowl can be minimized in the interior of a face or even strictly inside the feasible region if unconstrained there. QPs are central in ML: ridge regression, the soft-margin SVM dual, Markowitz portfolio selection, and model-predictive control are all QPs. Checking that P is PSD is the key step that certifies the problem is convex and solvable.

Worked Example 1

Problem. Identify P, q for the QP: minimize x1^2 + 2 x2^2 + x1 x2 - 3 x1, and confirm convexity.

  1. Match (1/2) x^T P x: the x1^2 term needs P11=2 (since 1/2*2=1), x2^2 needs P22=4, cross term x1 x2 needs P12=P21=1.
  2. So P = [[2,1],[1,4]], q = [-3, 0].
  3. Eigenvalues: `import numpy as np; np.linalg.eigvalsh([[2,1],[1,4]])` -> ~[1.59, 4.41], both > 0.
  4. P is positive definite, so the QP is (strictly) convex.

Answer. P=[[2,1],[1,4]], q=[-3,0]; P is positive definite, so this is a convex QP with a unique minimizer.

Worked Example 2

Problem. Solve the unconstrained convex QP minimize (1/2) x^T P x + q^T x with P=[[2,0],[0,2]], q=[-4,-6].

  1. Stationarity: gradient P x + q = 0 -> x = -P^{-1} q.
  2. `import numpy as np; P=np.array([[2,0],[0,2.]]); q=np.array([-4,-6.])`
  3. `x = np.linalg.solve(P, -q)` -> [2., 3.].
  4. Since P is PD, this stationary point is the global minimum.

Answer. Minimizer x = -P^{-1}q = [2, 3]; for a convex QP without active constraints the solution solves P x = -q.

Worked Example 3

Problem. Why is ridge regression a QP? Write min_w ||Xw - y||^2 + lambda ||w||^2 in QP form.

  1. Expand: ||Xw-y||^2 + lambda||w||^2 = w^T(X^T X + lambda I) w - 2 y^T X w + const.
  2. Match (1/2) w^T P w + q^T w: P = 2(X^T X + lambda I), q = -2 X^T y.
  3. X^T X is PSD and lambda I (lambda>0) is PD, so P is positive definite -> convex QP.
  4. Solution: w = (X^T X + lambda I)^{-1} X^T y, the closed-form ridge estimator.

Answer. Ridge regression is a convex QP with P = 2(X^T X + lambda I) (PD for lambda>0) and q = -2 X^T y; its optimum is the ridge normal-equation solution.

Common mistakes
  • Assuming every QP is convex. If P has a negative eigenvalue the QP is nonconvex and may have multiple local minima; check P is PSD first.
  • Expecting the QP optimum at a vertex like an LP. The quadratic objective can be minimized in the interior of a face or the whole region.
✎ Try it yourself

Problem. For the QP minimize 2x1^2 + 2x2^2 + 2 x1 x2 - 6 x1 - 4 x2, find P and q, verify convexity, and solve the unconstrained minimum.

Solution. Write as (1/2)x^T P x + q^T x. The x1^2 coefficient 2 needs P11=4, x2^2 needs P22=4, the cross term 2 x1 x2 needs P12=P21=2. So P=[[4,2],[2,4]], q=[-6,-4]. Eigenvalues of P are 2 and 6 (both >0), so P is PD and the QP is convex. Solve P x = -q: [[4,2],[2,4]] x = [6,4]. From 4x1+2x2=6 and 2x1+4x2=4: multiply the second by 2 (4x1+8x2=8), subtract the first to get 6x2=2 -> x2=1/3; then 4x1=6-2/3 -> x1=4/3. Minimizer x=(4/3, 1/3).

Modeling problems as LPs and QPs

The practical skill: recognize when a real decision is linear or quadratic and cast it in solver-ready form.

Modeling means translating a real objective and rules into the LP/QP template. Use an LP when both objective and constraints are linear: production planning, blending/diet, transportation, scheduling, and L1 (absolute-error) fitting after a standard reformulation. Use a QP when the objective has a quadratic term but constraints stay linear: least-squares with bounds, ridge regression, portfolio variance minimization, and the SVM margin. The core moves are: name decision variables, write a linear or quadratic objective, express limits as linear inequalities, and handle absolute values or max terms by introducing auxiliary variables (e.g. minimize t with -t <= expr <= t turns |expr| into an LP). Good modeling exposes the structure so you pick linprog, a QP solver, or CVXPY and get a globally optimal answer.

Worked Example 1

Problem. Model minimizing portfolio risk: minimize variance w^T Sigma w s.t. sum w = 1, w >= 0, given covariance Sigma.

  1. Decision variables: portfolio weights w (one per asset).
  2. Objective is quadratic: minimize w^T Sigma w (Sigma is a covariance matrix, hence PSD -> convex QP).
  3. Constraints are linear: sum(w) = 1 (fully invested), w >= 0 (no short-selling).
  4. QP form: P = 2 Sigma, q = 0, A_eq = ones, b_eq = 1, bounds w >= 0.

Answer. Minimum-variance portfolio is a convex QP: minimize w^T Sigma w s.t. 1^T w = 1, w >= 0; P=2 Sigma is PSD so it solves globally.

Worked Example 2

Problem. Cast L1 (least-absolute-deviations) regression min_w sum_i |x_i^T w - y_i| as an LP.

  1. Introduce auxiliary variables t_i >= 0 to bound each absolute residual.
  2. Add linear constraints: x_i^T w - y_i <= t_i and -(x_i^T w - y_i) <= t_i, so t_i >= |residual_i|.
  3. Objective becomes linear: minimize sum_i t_i over (w, t).
  4. This is a standard LP in the stacked variable [w; t], solvable by linprog.

Answer. L1 regression becomes an LP by replacing each |r_i| with t_i and the constraints -t_i <= r_i <= t_i, minimizing sum t_i — a classic absolute-value linearization.

Common mistakes
  • Trying to put |x| or max(...) directly into an LP. These are nonlinear; linearize with auxiliary variables and pairs of inequalities first.
  • Using a QP solver on a nonconvex objective. If the quadratic matrix is not PSD the problem is nonconvex; reformulate or use a global solver instead.
✎ Try it yourself

Problem. A diet must provide >= 60 protein and >= 40 fiber. Food A: 3 protein, 1 fiber, cost 2. Food B: 1 protein, 2 fiber, cost 3. Model the minimum-cost diet as an LP.

Solution. Variables: a, b >= 0 (servings). Objective (linear): minimize 2a + 3b. Constraints (linear): 3a + b >= 60 (protein), a + 2b >= 40 (fiber), a,b >= 0. This is an LP. In scipy: linprog([2,3], A_ub=[[-3,-1],[-1,-2]], b_ub=[-60,-40], bounds=[(0,None),(0,None)]). The optimum sits at the vertex where both constraints bind: solving 3a+b=60 and a+2b=40 gives a=16, b=12, cost 2*16+3*12=68.

Solving LPs and QPs with SciPy and CVXPY

Two tools: scipy.optimize.linprog for LPs and CVXPY for both LPs and QPs with readable, near-mathematical syntax.

SciPy's linprog solves LPs given c, A_ub, b_ub, A_eq, b_eq, and bounds, returning an OptimizeResult with .x, .fun, and constraint marginals (duals). For QPs and richer convex problems, CVXPY lets you declare variables and write the objective and constraints almost as you would on paper; it then compiles the problem, dispatches to an appropriate solver, and returns the optimum. CVXPY enforces disciplined convex programming, so if it accepts your model it is provably convex and the answer is globally optimal. The practical workflow is: write the model, call prob.solve(), then read variable.value and prob.value, and check prob.status == 'optimal'. Mastering these tools is what makes the modeling skills of this unit actually executable for real ML and operations problems.

Worked Example 1

Problem. Solve a small LP with CVXPY: maximize 3x + 2y s.t. x + y <= 4, x <= 3, x,y >= 0.

  1. `import cvxpy as cp; x = cp.Variable(); y = cp.Variable()`
  2. `prob = cp.Problem(cp.Maximize(3*x + 2*y), [x + y <= 4, x <= 3, x >= 0, y >= 0])`
  3. `prob.solve()`
  4. `x.value, y.value` -> ~(3.0, 1.0); `prob.value` -> ~11.0; `prob.status` -> 'optimal'.

Answer. CVXPY returns x=3, y=1, objective 11 with status 'optimal' — note CVXPY handles the maximize directly, no manual negation.

Worked Example 2

Problem. Solve a convex QP with CVXPY: minimize (1/2)(x1^2 + x2^2) + x1 s.t. x1 + x2 = 1, x >= 0.

  1. `import cvxpy as cp; x = cp.Variable(2, nonneg=True)`
  2. `obj = cp.Minimize(0.5*cp.sum_squares(x) + x[0])`
  3. `prob = cp.Problem(obj, [cp.sum(x) == 1]); prob.solve()`
  4. `x.value` -> ~[0., 1.]; the linear +x1 term pushes mass onto x2; `prob.value` -> ~0.5.

Answer. CVXPY solves the QP to x ~ [0, 1] with optimal value ~0.5; sum_squares keeps the model disciplined-convex.

Worked Example 3

Problem. Solve the minimum-variance portfolio QP in CVXPY for Sigma=[[0.1,0.02],[0.02,0.08]].

  1. `import cvxpy as cp, numpy as np; S = np.array([[0.1,0.02],[0.02,0.08]])`
  2. `w = cp.Variable(2, nonneg=True)`
  3. `prob = cp.Problem(cp.Minimize(cp.quad_form(w, S)), [cp.sum(w) == 1]); prob.solve()`
  4. `w.value` -> weights summing to 1 that minimize variance; the lower-variance asset (asset 2) gets the larger weight.

Answer. CVXPY's quad_form encodes w^T Sigma w; it returns the fully-invested, long-only weights minimizing portfolio variance, with status 'optimal'.

Common mistakes
  • Reading .x or .value without checking status. A problem can be 'infeasible' or 'unbounded'; always confirm prob.status == 'optimal' (or res.success) first.
  • Writing a non-DCP objective in CVXPY (e.g. a nonconvex quadratic). CVXPY rejects it; use cp.quad_form with a PSD matrix or cp.sum_squares to stay disciplined-convex.
✎ Try it yourself

Problem. Use CVXPY to solve: minimize x^2 + y^2 subject to x + y >= 2, x,y >= 0, and report the optimal point and value.

Solution. import cvxpy as cp; x = cp.Variable(nonneg=True); y = cp.Variable(nonneg=True); prob = cp.Problem(cp.Minimize(x**2 + y**2), [x + y >= 2]); prob.solve(). The constraint binds at x+y=2 and symmetry/convexity gives x=y=1. So x.value ~ 1.0, y.value ~ 1.0, prob.value ~ 2.0, prob.status 'optimal'. (Check: the unconstrained min is (0,0) but it violates x+y>=2, so the optimum lies on the active line x+y=2 at its closest point to the origin, (1,1).)

Key terms
  • Linear program (LP) — minimize a linear objective subject to linear equality and inequality constraints.
  • Polytope — the feasible region of an LP, a region bounded by flat faces (hyperplanes).
  • Vertex (extreme point) — a corner of the polytope; an LP optimum is always attained at a vertex if one exists.
  • Simplex method — an algorithm that walks vertex-to-vertex along edges to an optimal corner.
  • Interior-point method — an algorithm that approaches the optimum through the interior of the feasible region.
  • LP duality — every LP has a dual LP whose optimal value equals the primal's (strong duality).
  • Quadratic program (QP) — minimize a convex quadratic objective subject to linear constraints.
  • Sensitivity analysis — how the optimal value changes as constraint bounds or costs are perturbed, read from the duals.
Assignment · Model and solve a planning problem

Formulate a small resource-allocation problem (for example, a diet or production-mix problem) as a linear program in standard form. Solve it with scipy.optimize.linprog and again with CVXPY, confirming the objective values match. Then add a quadratic penalty (turning it into a QP, e.g. least-squares with bounds) and solve it in CVXPY. Use the dual values to do a one-line sensitivity analysis on a binding constraint.

Deliverable · A notebook with the LP and QP formulations, both solver solutions, and a short sensitivity interpretation from the dual values.

Quiz · 5 questions
  1. 1. The optimum of a bounded, feasible linear program is always attained:

  2. 2. A quadratic program differs from a linear program because it has:

  3. 3. Strong duality for a feasible LP means:

  4. 4. The simplex method finds an LP optimum by:

  5. 5. Which SciPy function is intended for solving linear programs?

You'll be able to

I can put a linear program in standard form and describe its feasible region as a polytope.

I can model a resource-allocation or fitting problem as an LP or QP.

I can solve LPs and QPs with SciPy or CVXPY and interpret the dual/sensitivity output.

Weeks 11-12 Unit 6: Convex Optimization Essentials
define-convex-programapply-convexity-rulesrecognize-convex-objectivesform-dual-problemexplain-duality-gapmodel-with-cvxpy
Lecture
The convex optimization problem and disciplined convex programming

A convex problem is the gold standard of tractability; disciplined convex programming (DCP) is the rule system that keeps your models in that safe zone.

A convex optimization problem minimizes a convex objective over a convex feasible set: minimize f0(x) s.t. fi(x) <= 0 with each fi convex and equalities Ax = b affine. Its great property is that any local optimum is global, so solvers return a certified best answer in polynomial time. Disciplined convex programming (DCP), the ruleset behind CVXPY, builds expressions only from a library of atoms with known curvature (convex, concave, affine) and sign, combined by rules that provably preserve convexity. If a model passes DCP, it is guaranteed convex and solvable; if it fails, CVXPY refuses rather than returning a wrong answer. This discipline matters in ML because it lets you compose losses, norms, and constraints confidently and get globally optimal estimators without proving convexity by hand each time.

Worked Example 1

Problem. Is 'minimize cp.square(x) + cp.abs(y) s.t. x + y == 1' DCP-valid? Explain.

  1. cp.square(x) is a convex atom; cp.abs(y) is a convex atom.
  2. A sum of convex functions is convex, so the objective is convex.
  3. The constraint x + y == 1 is affine (allowed as an equality).
  4. `import cvxpy as cp; x=cp.Variable(); y=cp.Variable(); cp.Problem(cp.Minimize(cp.square(x)+cp.abs(y)),[x+y==1]).is_dcp()` -> True.

Answer. Yes, it is DCP-valid: a sum of convex atoms minimized over an affine equality, so CVXPY accepts and solves it globally.

Worked Example 2

Problem. Why does 'minimize cp.square(x) * cp.square(y)' fail DCP even though each factor is convex?

  1. DCP needs a known-curvature composition rule; a product of two convex functions has no general curvature.
  2. square(x)*square(y) = x^2 y^2 is nonconvex (e.g. its Hessian is indefinite away from the origin).
  3. DCP only certifies products when one factor is a known constant/affine sign-definite term.
  4. `cp.Problem(cp.Minimize(cp.square(x)*cp.square(y))).is_dcp()` -> False; CVXPY rejects it.

Answer. It fails DCP because a product of two nonconstant convex terms is not provably convex; CVXPY returns is_dcp() False rather than risk a wrong solve.

Common mistakes
  • Assuming any expression that 'looks bowl-shaped' is DCP. DCP follows strict composition rules; verify with .is_dcp() rather than intuition.
  • Writing an equality constraint that is nonlinear. Convex problems require equality constraints to be affine; a nonlinear equality breaks convexity.
✎ Try it yourself

Problem. Decide whether 'minimize cp.norm(A @ x - b, 2) subject to cp.norm(x, 1) <= 1' is a DCP convex problem, and say why.

Solution. cp.norm(A@x-b, 2) is the composition of the convex 2-norm atom with an affine argument A@x-b, so the objective is convex. The constraint uses the convex 1-norm of x bounded above by a constant (cp.norm(x,1) <= 1), which defines a convex sublevel set. Objective convex, constraint set convex -> it is a DCP-valid convex problem (a LASSO-type / norm-ball problem). In code, cp.Problem(cp.Minimize(cp.norm(A@x-b,2)), [cp.norm(x,1) <= 1]).is_dcp() returns True, so CVXPY solves it to a global optimum.

Operations that preserve convexity

You rarely prove convexity from scratch; instead you build convex functions from convex pieces using a small set of safe operations.

Convexity is preserved under specific operations, and chaining them is how DCP works. Key rules: a nonnegative weighted sum of convex functions is convex; a pointwise maximum (or supremum) of convex functions is convex; composition with an affine map, f(Ax + b), preserves convexity; and certain scalar compositions are convex (e.g. h(g(x)) is convex when h is convex-and-nondecreasing and g is convex). Norms are convex; the max of linear functions is convex (piecewise-linear). Conversely, a pointwise minimum of convex functions is generally not convex. Knowing these rules lets you certify complex objectives — like a regularized hinge loss sum or a max over linear scores — as convex by inspection, which is exactly how losses in SVMs, robust regression, and many ML models are shown solvable.

Worked Example 1

Problem. Show f(x) = max(0, 1 - x) (the hinge) is convex.

  1. 0 is constant (convex/affine); 1 - x is affine (convex).
  2. The pointwise maximum of convex functions is convex.
  3. Therefore max(0, 1 - x) is convex (it is piecewise-linear with a single kink at x=1).
  4. `import cvxpy as cp; x=cp.Variable(); cp.maximum(0, 1-x).is_convex()` -> True.

Answer. The hinge max(0, 1-x) is convex as a pointwise max of two affine functions; this is why SVM training is a convex problem.

Worked Example 2

Problem. Is g(x) = (sum_i x_i^2)^2 (square of a convex function) convex? Use the composition rule.

  1. Inner g0(x) = sum x_i^2 is convex and nonnegative.
  2. Outer h(u) = u^2 is convex AND nondecreasing for u >= 0.
  3. Composition h(g0(x)) with h convex-nondecreasing and g0 convex (and nonnegative) is convex.
  4. Numerically, its Hessian is PSD on the domain; CVXPY accepts cp.square(cp.sum_squares(x)).

Answer. Yes, convex: h(u)=u^2 is convex and nondecreasing on u>=0 and the inner sum-of-squares is convex and nonnegative, so the composition is convex.

Worked Example 3

Problem. Show that the pointwise MINIMUM of two convex functions can fail to be convex.

  1. Take f1(x) = (x+2)^2 and f2(x) = (x-2)^2, both convex bowls.
  2. Their minimum m(x) = min(f1, f2) is a 'W'-shaped curve with two valleys at x=-2 and x=2.
  3. The chord between (-2, 0) and (2, 0) lies below the bump at x=0 where m(0)=4 — so a chord can dip below the function: m is nonconvex.
  4. Rule: max of convex is convex, but min of convex is not in general.

Answer. min((x+2)^2,(x-2)^2) has two separate minima and a hump between them, so it is nonconvex — the minimum operation does not preserve convexity.

Common mistakes
  • Applying the scalar composition rule without the monotonicity condition. h(g(x)) is convex only if h is convex AND nondecreasing (for convex g); skipping that can give a wrong verdict.
  • Assuming min of convex functions is convex like max. Pointwise maximum preserves convexity; pointwise minimum generally does not.
✎ Try it yourself

Problem. Use convexity-preserving rules to argue that f(x) = ||Ax - b||^2 + lambda ||x||_1 (lambda >= 0) is convex.

Solution. ||Ax-b|| is the composition of the convex 2-norm with the affine map Ax-b, hence convex; squaring it via h(u)=u^2 (convex, nondecreasing on u>=0) keeps it convex, so ||Ax-b||^2 is convex. ||x||_1 is a norm, hence convex. A nonnegative weighted sum of convex functions is convex, and lambda>=0, so ||Ax-b||^2 + lambda||x||_1 is convex. This is the LASSO objective, certified convex purely by composition rules, so any convex solver (or CVXPY) finds its global optimum.

Norms, regularizers, and common convex objectives

Norms are the building blocks of regularizers; their geometry explains why L1 gives sparsity and L2 gives smooth shrinkage.

A norm measures size and is always convex; the common ones are the L2 (Euclidean) norm sqrt(sum x_i^2), the L1 norm sum |x_i|, and the L-infinity norm max |x_i|. Used as regularizers added to a loss, they pull solutions toward smaller, simpler parameter vectors. The geometry matters: the L1 ball is a diamond with sharp corners on the axes, so the optimum of loss-plus-L1 tends to land at a corner where some coordinates are exactly zero — that is why the LASSO produces sparse models (feature selection). The L2 ball is a smooth sphere, so ridge regularization shrinks all coordinates smoothly without forcing zeros. Common convex objectives — least squares, ridge, LASSO, logistic loss, hinge loss — are all convex sums of these pieces, which is why they train reliably to a global optimum.

Worked Example 1

Problem. Compute the L1, L2, and L-infinity norms of x = (3, -4, 0).

  1. L1: |3| + |-4| + |0| = 7.
  2. L2: sqrt(9 + 16 + 0) = sqrt(25) = 5.
  3. L-infinity: max(|3|, |-4|, |0|) = 4.
  4. `import numpy as np; x=np.array([3,-4,0.]); np.linalg.norm(x,1), np.linalg.norm(x,2), np.linalg.norm(x,np.inf)` -> (7.0, 5.0, 4.0).

Answer. L1 = 7, L2 = 5, L-infinity = 4; each norm is convex and emphasizes a different notion of vector size.

Worked Example 2

Problem. Explain geometrically why minimizing loss + lambda||w||_1 produces zeros while loss + lambda||w||_2^2 does not.

  1. Picture the loss contours (ellipses) expanding until they touch the regularizer's level set (norm ball).
  2. The L1 ball is a diamond whose corners lie on the axes (where coordinates are zero).
  3. Expanding ellipses most often first touch a corner -> the solution has exact zeros (sparse).
  4. The L2 ball is round with no corners, so contact happens at a generic point with all coordinates small but nonzero (shrinkage, not sparsity).

Answer. L1's diamond corners sit on the axes, so the optimum is pushed onto a zero coordinate (sparsity); L2's smooth sphere yields uniform shrinkage with no exact zeros.

Worked Example 3

Problem. Set up the LASSO objective in CVXPY and solve a tiny instance.

  1. `import cvxpy as cp, numpy as np; A=np.array([[1,0],[0,1],[1,1.]]); b=np.array([1,0,0.8]); lam=0.5`
  2. `w = cp.Variable(2)`
  3. `prob = cp.Problem(cp.Minimize(cp.sum_squares(A@w - b) + lam*cp.norm(w,1))); prob.solve()`
  4. `w.value` -> small/possibly-zero coefficients; the L1 term drives weak coordinates to exactly 0 as lam grows.

Answer. CVXPY minimizes sum_squares(Aw-b) + lam*||w||_1; the L1 penalty yields a sparse w, and increasing lam zeros out more coefficients — the LASSO feature-selection effect.

Common mistakes
  • Penalizing the bias/intercept term. Regularizers should usually apply to weights only; shrinking the intercept biases predictions.
  • Forgetting to scale features before L1/L2 regularization. Unscaled features make the penalty hit large-scale coordinates unfairly; standardize first.
✎ Try it yourself

Problem. For w = (0, -2, 1, -1), compute ||w||_1, ||w||_2, and ||w||_inf, and say which regularizer would most encourage exact zeros.

Solution. ||w||_1 = 0 + 2 + 1 + 1 = 4. ||w||_2 = sqrt(0 + 4 + 1 + 1) = sqrt(6) ~ 2.449. ||w||_inf = max(0,2,1,1) = 2. The L1 norm (sum of absolute values) most encourages exact zeros: its diamond-shaped ball has corners on the coordinate axes, so loss-plus-L1 tends to set small coefficients to exactly 0, producing a sparse solution, whereas L2 only shrinks them toward zero without eliminating them.

Lagrangian duality and the dual problem

Duality builds a companion problem whose optimum bounds the original; for convex problems the two values often coincide exactly.

Given minimize f0(x) s.t. fi(x) <= 0, hj(x) = 0, the Lagrangian is L(x, lambda, nu) = f0(x) + sum lambda_i fi(x) + sum nu_j hj(x) with lambda_i >= 0. The Lagrange dual function g(lambda, nu) = inf_x L is the minimum of the Lagrangian over x; it is always concave (an infimum of affine functions of the multipliers) regardless of the primal. The dual problem maximizes g over lambda >= 0, nu. Its optimal value d* is always a lower bound on the primal p* (weak duality). Duality is powerful: it produces certificates of optimality, turns constrained problems into often-simpler dual problems, and reveals structure — the SVM is usually solved in its dual form, and the dual multipliers are the support-vector weights. Building and reading the dual is a core convex-optimization skill.

Worked Example 1

Problem. Form the dual function for minimize x^2 s.t. x >= 1 (i.e. 1 - x <= 0).

  1. Lagrangian: L(x, lambda) = x^2 + lambda(1 - x), lambda >= 0.
  2. Minimize over x: dL/dx = 2x - lambda = 0 -> x = lambda/2.
  3. Dual function g(lambda) = (lambda/2)^2 + lambda(1 - lambda/2) = lambda - lambda^2/4.
  4. Maximize g over lambda >= 0: g'(lambda) = 1 - lambda/2 = 0 -> lambda = 2, g(2) = 2 - 1 = 1.

Answer. Dual function g(lambda) = lambda - lambda^2/4, maximized at lambda=2 giving d* = 1, which matches the primal optimum x=1, f=1.

Worked Example 2

Problem. Why is the dual function always concave, and what does that buy us?

  1. For each fixed x, L(x, lambda, nu) is affine (linear) in (lambda, nu).
  2. g(lambda, nu) = inf_x L is a pointwise infimum of affine functions.
  3. An infimum of affine functions is always concave.
  4. So maximizing the dual is a convex problem even when the primal is not — giving a tractable lower bound for hard problems.

Answer. Because g is an infimum of affine functions of the multipliers, it is concave; thus the dual is always a convex (maximize-concave) problem, useful even for nonconvex primals.

Common mistakes
  • Dropping the sign requirement lambda >= 0 on inequality multipliers. The dual is maximized only over lambda >= 0 (equality multipliers nu are free).
  • Forgetting the dual gives only a lower bound in general. Equality of primal and dual values needs strong duality (next lesson), not just weak duality.
✎ Try it yourself

Problem. Form and solve the Lagrange dual of: minimize x^2 + y^2 subject to x + y = 1.

Solution. Lagrangian L(x,y,nu) = x^2 + y^2 + nu(x + y - 1) (equality multiplier nu is free). Minimize over x,y: dL/dx = 2x + nu = 0 -> x = -nu/2; similarly y = -nu/2. Substitute into L: g(nu) = 2*(nu/2)^2 + nu(-nu - 1) = nu^2/2 - nu^2 - nu = -nu^2/2 - nu. Maximize: g'(nu) = -nu - 1 = 0 -> nu = -1, g(-1) = -0.5 + 1 = 0.5. So d* = 0.5, matching the primal optimum x=y=0.5 with value 0.5 (strong duality holds for this convex problem).

Weak vs. strong duality and Slater's condition

Weak duality always holds; strong duality (no gap) is what makes the dual exactly solve the primal, and Slater's condition guarantees it for convex problems.

Weak duality always holds: the dual optimum d* is a lower bound on the primal optimum p*, so p* - d* = the duality gap >= 0. Strong duality means the gap is zero, p* = d*, so solving the dual solves the primal. Strong duality is not automatic, but for a convex problem it holds under a constraint qualification — most commonly Slater's condition: there exists a strictly feasible point where all inequality constraints hold strictly (fi(x) < 0) and equalities are met. When Slater holds, p* = d* and the optimal multipliers exist, which is the theoretical backbone of KKT-based and interior-point solvers. This is why convex problems are so well-behaved: you can solve whichever of primal or dual is easier and trust the values match — exploited heavily in SVM training and convex ML estimators.

Worked Example 1

Problem. Check Slater's condition for minimize x^2 s.t. x >= 1, and conclude about strong duality.

  1. Rewrite the constraint as 1 - x <= 0; it is convex (affine), objective x^2 is convex.
  2. Slater needs a strictly feasible point: pick x = 2, then 1 - 2 = -1 < 0 (strict).
  3. A strictly feasible point exists, so Slater's condition holds.
  4. Therefore strong duality holds: indeed p* = 1 (at x=1) equals d* = 1 found earlier.

Answer. Slater holds (x=2 is strictly feasible), so strong duality holds and the dual value d*=1 equals the primal p*=1 — zero duality gap.

Worked Example 2

Problem. Give the duality gap for a problem where p* = 5 and the best dual value found is d* = 3, and interpret it.

  1. Duality gap = p* - d* = 5 - 3 = 2.
  2. By weak duality the gap is >= 0, consistent here.
  3. A positive gap means strong duality fails (or the dual is not yet maximized).
  4. If the problem were convex with Slater satisfied, the gap would be 0; a gap of 2 signals nonconvexity or a missing constraint qualification.

Answer. The duality gap is 2; since it is nonzero, strong duality does not hold here (the problem is likely nonconvex or fails a constraint qualification).

Common mistakes
  • Assuming strong duality for any problem. It is guaranteed for convex problems only under a constraint qualification like Slater's; nonconvex problems can have a positive gap.
  • Treating an interior point as Slater-satisfying when an inequality holds only with equality. Slater requires strict inequality fi(x) < 0, not fi(x) <= 0, at the witness point.
✎ Try it yourself

Problem. Does Slater's condition hold for minimize x^2 + y^2 subject to x + y <= -1? If so, what does it imply?

Solution. The objective x^2+y^2 is convex and the constraint x+y+1 <= 0 is affine (convex). Slater needs a strictly feasible point: pick (x,y)=(-1,-1), then x+y = -2 < -1, so x+y+1 = -1 < 0 strictly. A strictly feasible point exists, so Slater's condition holds. This implies strong duality: the primal and dual optimal values are equal (zero duality gap) and optimal Lagrange multipliers exist. (The primal optimum is the point on x+y=-1 nearest the origin, (-0.5,-0.5), with value 0.5.)

Modeling and solving with CVXPY

The capstone tool: express a convex problem in near-mathematical CVXPY syntax, solve it, and read a certified-optimal answer.

CVXPY is a Python modeling language for convex optimization. You declare cp.Variable objects, write a cp.Minimize or cp.Maximize objective from DCP atoms (sum_squares, norm, quad_form, log, entr, etc.), list constraints, build a cp.Problem, and call .solve(). CVXPY verifies the model is DCP-convex, transforms it to a standard conic form, dispatches to a solver (ECOS, OSQP, SCS, Clarabel), and returns the optimum; you read variable.value, prob.value, and check prob.status == 'optimal'. The workflow scales from least squares and LASSO to SVMs, portfolio optimization, and control. Mastering CVXPY ties together the whole unit: norms and regularizers become atoms, convexity-preserving operations become valid compositions, and duality is available via constraint.dual_value. It is the practical bridge from convex theory to working, globally-optimal ML and decision models.

Worked Example 1

Problem. Solve regularized least squares (ridge): minimize ||Ax - b||^2 + ||x||^2 in CVXPY for a small A, b.

  1. `import cvxpy as cp, numpy as np; A=np.array([[1,1],[1,-1],[2,1.]]); b=np.array([2,0,3.])`
  2. `x = cp.Variable(2)`
  3. `prob = cp.Problem(cp.Minimize(cp.sum_squares(A@x - b) + cp.sum_squares(x))); prob.solve()`
  4. `x.value` -> the ridge solution (X^T X + I)^{-1} X^T b; `prob.status` -> 'optimal'.

Answer. CVXPY returns the ridge estimate x = (A^T A + I)^{-1} A^T b with status 'optimal'; sum_squares keeps the model disciplined-convex.

Worked Example 2

Problem. Model a soft-margin SVM-style objective in CVXPY: minimize (1/2)||w||^2 + C*sum(hinge) for labels y, features X.

  1. `import cvxpy as cp, numpy as np; w=cp.Variable(d); b=cp.Variable(); C=1.0`
  2. `scores = cp.multiply(y, X@w + b)` # y in {-1,+1}
  3. `loss = cp.sum(cp.pos(1 - scores))` # hinge: max(0, 1 - margin)
  4. `prob = cp.Problem(cp.Minimize(0.5*cp.sum_squares(w) + C*loss)); prob.solve()`

Answer. The SVM objective is convex (squared norm plus a sum of hinges); CVXPY's cp.pos encodes the hinge and returns the optimal (w, b) — note SVMs are convex, so the solution is global.

Worked Example 3

Problem. Retrieve a dual value (shadow price) from a solved CVXPY problem.

  1. `import cvxpy as cp; x = cp.Variable(); con = (x >= 1)`
  2. `prob = cp.Problem(cp.Minimize(cp.square(x)), [con]); prob.solve()`
  3. `x.value` -> ~1.0; `prob.value` -> ~1.0.
  4. `con.dual_value` -> ~2.0, the Lagrange multiplier (sensitivity of the optimum to relaxing x >= 1).

Answer. CVXPY exposes each constraint's optimal multiplier via constraint.dual_value (here ~2), giving the shadow price directly from the solve.

Common mistakes
  • Reading variable.value before checking prob.status. If the status is 'infeasible' or 'unbounded', the values are None or meaningless.
  • Forcing a nonconvex model through CVXPY. If .is_dcp() is False, reformulate (use sum_squares, quad_form with a PSD matrix, pos for hinge) instead of hacking around the DCP error.
✎ Try it yourself

Problem. Use CVXPY to solve: minimize ||x||_1 subject to A x = b, with A = [[1, 1, 1]], b = [1] (an underdetermined sparse-recovery problem), and describe the expected solution.

Solution. import cvxpy as cp, numpy as np; A = np.array([[1.,1.,1.]]); b = np.array([1.]); x = cp.Variable(3); prob = cp.Problem(cp.Minimize(cp.norm(x,1)), [A@x == b]); prob.solve(). The objective ||x||_1 is convex and the constraint A x = b is affine, so this is a convex (DCP) problem. Minimizing the L1 norm subject to a single equality gives a sparse solution: it places all the mass on one coordinate, e.g. x ~ [1, 0, 0] (any single-coordinate solution achieves ||x||_1 = 1, the minimum). prob.value ~ 1.0, status 'optimal'. This is the basis-pursuit / compressed-sensing principle that L1 minimization recovers sparse vectors.

Key terms
  • Convex program — minimizing a convex objective over a convex feasible set; local optima are global.
  • Disciplined convex programming (DCP) — a rule set (used by CVXPY) that builds problems guaranteed to be convex.
  • Convexity-preserving operations — non-negative weighted sums, pointwise maxima, and composition rules that keep functions convex.
  • Norm — a convex measure of size (such as L1 or L2) widely used as objective or regularizer.
  • Lagrangian dual — the function obtained by minimizing the Lagrangian over x, giving a lower bound on the primal.
  • Weak duality — the dual optimal value is always a lower bound on the primal optimal value.
  • Strong duality — the primal and dual optimal values are equal, so the duality gap is zero.
  • Slater's condition — a strictly feasible point exists, which guarantees strong duality for a convex problem.
Assignment · Convex modeling with CVXPY

Pick a convex problem such as L1-regularized (lasso) or L2-regularized least squares, or a robust regression. Express it in CVXPY using DCP-compliant atoms and solve it for several regularization strengths. Verify your model is recognized as convex by CVXPY, plot how the solution changes with the regularizer, and inspect the dual variables of any constraints. Briefly state why each term keeps the problem convex.

Deliverable · A CVXPY notebook with the convex model, a regularization-path plot, the reported dual values, and a short convexity justification.

Quiz · 5 questions
  1. 1. In a convex optimization problem, the objective is convex and the feasible set is:

  2. 2. Which operation preserves convexity?

  3. 3. Weak duality guarantees that the dual optimal value is:

  4. 4. Slater's condition is important because it:

  5. 5. CVXPY enforces disciplined convex programming (DCP) in order to:

You'll be able to

I can recognize a convex program and use convexity-preserving rules to build one.

I can form the Lagrangian dual and explain the duality gap and when it is zero.

I can express and solve a convex problem in CVXPY using disciplined convex programming.

Weeks 13-14 Unit 7: Optimization in Machine Learning
frame-ermselect-loss-functiondescribe-loss-landscapeapply-regularizationreason-saddle-localdiagnose-trainingtrain-and-tune-model
Lecture
Learning as optimization: empirical risk minimization

Training a model is literally an optimization problem — minimize average loss over the data — which is what ties this whole course to machine learning.

Empirical risk minimization (ERM) frames supervised learning as optimization: choose parameters w to minimize the average loss over the training set, R_emp(w) = (1/n) sum_i L(f_w(x_i), y_i). The true goal is to minimize the expected (population) risk, but since the data distribution is unknown we minimize the empirical risk as a surrogate, often adding a regularizer to control overfitting. This makes every learning algorithm an instance of the optimization machinery built earlier: pick a model f_w, a loss L, run gradient-based or convex solvers to drive grad R_emp toward zero. Whether the resulting problem is easy (convex, like linear/logistic regression) or hard (nonconvex, like deep nets) is exactly the convexity question from Unit 1. ERM is the unifying lens: data + loss + regularizer = an optimization problem.

Worked Example 1

Problem. Write the ERM objective for linear regression with MSE loss on data (x_i, y_i), i=1..n.

  1. Model: f_w(x) = w^T x (absorb bias into x).
  2. Loss: squared error L = (f_w(x_i) - y_i)^2.
  3. Empirical risk: R(w) = (1/n) sum_i (w^T x_i - y_i)^2 = (1/n)||Xw - y||^2.
  4. `import numpy as np; w = np.linalg.lstsq(X, y, rcond=None)[0]` minimizes R(w) in closed form.

Answer. ERM objective R(w) = (1/n)||Xw - y||^2; it is convex (Hessian (2/n)X^T X is PSD), so least squares solves it globally.

Worked Example 2

Problem. Add L2 regularization (ridge) to the ERM objective and write the gradient for gradient descent.

  1. Regularized risk: R(w) = (1/n)||Xw - y||^2 + lambda ||w||^2.
  2. Gradient: grad R = (2/n) X^T (Xw - y) + 2 lambda w.
  3. Gradient-descent step: `w -= eta * ((2/len(y))*X.T@(X@w - y) + 2*lam*w)`.
  4. Setting grad=0 gives the closed form w = (X^T X + n*lam*I)^{-1} X^T y.

Answer. Regularized ERM: R(w) = (1/n)||Xw-y||^2 + lambda||w||^2 with gradient (2/n)X^T(Xw-y) + 2 lambda w; convex, so gradient descent reaches the global ridge solution.

Common mistakes
  • Confusing empirical risk with true risk. Minimizing training loss perfectly can overfit; the goal is low expected risk, approximated via regularization and validation.
  • Forgetting the 1/n averaging. Using a sum instead of a mean changes the effective learning rate and the scale of the regularization trade-off.
✎ Try it yourself

Problem. Write the empirical-risk objective for logistic regression with the log-loss on labels y_i in {0,1}, and state whether it is convex.

Solution. Model probability p_i = sigma(w^T x_i) with sigma the logistic function. Log-loss: L_i = -[y_i log p_i + (1-y_i) log(1-p_i)]. Empirical risk R(w) = (1/n) sum_i L_i. This objective is convex in w: the negative log-likelihood of logistic regression has Hessian X^T D X with D a diagonal of p_i(1-p_i) >= 0, which is positive semidefinite, so R(w) is convex. Therefore gradient descent (or Newton/IRLS) converges to the global minimum, and adding lambda||w||^2 keeps it strictly convex.

Common loss functions and their shapes

The loss you choose decides the optimization landscape and what the model is rewarded for; knowing each shape is essential.

Different tasks use different losses, each with a characteristic shape that drives optimization. Squared error (regression) is a smooth convex parabola that penalizes large residuals quadratically — sensitive to outliers. Absolute error (L1) is convex but kinked at zero, robust to outliers, and yields a non-smooth problem. The Huber loss blends the two: quadratic near zero, linear far out, giving robustness with smoothness. For classification, log-loss (cross-entropy) is smooth and convex, heavily penalizing confident wrong predictions, while hinge loss (SVM) is convex and piecewise-linear, zero once the margin is satisfied. All these losses are convex in the model's score, which is why linear models train reliably; the nonconvexity in deep learning comes from composing them with nonlinear networks. Matching loss shape to the task and to outlier behavior is a core modeling decision.

Worked Example 1

Problem. Compute squared, absolute, and Huber (delta=1) loss for a residual r = 3.

  1. Squared: r^2 = 9.
  2. Absolute: |r| = 3.
  3. Huber with delta=1: since |r|=3 > delta, use delta*(|r| - 0.5*delta) = 1*(3 - 0.5) = 2.5.
  4. Squared (9) >> absolute (3) > Huber (2.5): squared loss most amplifies the large residual.

Answer. Squared = 9, absolute = 3, Huber = 2.5; squared loss penalizes the outlier far more, motivating robust losses like Huber/L1.

Worked Example 2

Problem. Compute the hinge loss for a correctly-classified point with margin y*score = 1.5 and for a violated one with y*score = -0.2.

  1. Hinge loss = max(0, 1 - y*score).
  2. Margin 1.5: max(0, 1 - 1.5) = max(0, -0.5) = 0 (no penalty, margin satisfied).
  3. Margin -0.2: max(0, 1 - (-0.2)) = max(0, 1.2) = 1.2 (penalized).
  4. `import numpy as np; np.maximum(0, 1 - np.array([1.5, -0.2]))` -> [0., 1.2].

Answer. Hinge loss is 0 for the well-classified point (margin >= 1) and 1.2 for the violated one; only margin violations are penalized.

Worked Example 3

Problem. Compute log-loss for a true label y=1 when the model predicts probability p=0.9 vs. p=0.1.

  1. Log-loss for y=1 is -log(p).
  2. p=0.9: -log(0.9) ~ 0.105 (small penalty, confident and correct).
  3. p=0.1: -log(0.1) ~ 2.303 (large penalty, confident and wrong).
  4. `import numpy as np; -np.log([0.9, 0.1])` -> [0.105, 2.303].

Answer. Log-loss is ~0.105 for the confident correct prediction and ~2.303 for the confident wrong one; cross-entropy harshly punishes confident mistakes.

Common mistakes
  • Using squared error when data has outliers. Its quadratic growth lets a few outliers dominate the fit; prefer Huber or L1 for robustness.
  • Treating hinge and log-loss as interchangeable. Hinge is zero past the margin (sparse support vectors); log-loss always has a nonzero gradient, giving calibrated probabilities.
✎ Try it yourself

Problem. For residuals r = [0.5, 4], compute the total squared loss and total Huber loss (delta=1), and explain the difference.

Solution. Squared: 0.5^2 + 4^2 = 0.25 + 16 = 16.25. Huber (delta=1): for r=0.5, |r| <= delta so loss = 0.5*r^2 = 0.5*0.25 = 0.125; for r=4, |r| > delta so loss = delta*(|r| - 0.5*delta) = 1*(4 - 0.5) = 3.5. Total Huber = 0.125 + 3.5 = 3.625. The squared loss (16.25) is dominated by the single large residual (16 of it), whereas Huber grows only linearly past delta, so the outlier contributes 3.5 instead of 16. Huber is far more robust to the outlier while staying smooth near zero.

Loss landscapes: convex models vs. neural networks

The geometry of the loss surface determines whether training is a solved problem or a search — this is the convex/nonconvex divide in action.

A loss landscape is the surface of loss values over the parameter space. For convex models (linear/logistic regression, SVM, LASSO) the landscape is a single convex bowl: one connected set of global minima, no traps, so any descent method finds the best parameters and the answer is essentially unique up to flat directions. Neural networks produce highly nonconvex landscapes: many local minima, saddle points, plateaus, and symmetries (permuting hidden units leaves the loss unchanged). Remarkably, in high dimensions most critical points are saddles rather than bad local minima, and many local minima have similar low loss, so SGD with good initialization and momentum reliably finds a good (not provably global) solution. Understanding which regime you are in tells you whether to trust convergence (convex) or to care about initialization, restarts, and architecture (nonconvex).

Worked Example 1

Problem. Contrast the landscape of L(w) = ||Xw - y||^2 with a two-layer net loss, in terms of minima.

  1. Least squares: Hessian 2 X^T X is PSD everywhere -> convex bowl, single global minimum (or a flat affine set of them).
  2. Any solver from any start reaches the same optimal loss.
  3. Two-layer net: composing weights through a nonlinearity makes the loss nonconvex with many critical points.
  4. Different initializations land in different basins, though many achieve comparably low loss.

Answer. The least-squares landscape is a single convex bowl (global optimum guaranteed); the neural-net landscape is nonconvex with many minima/saddles, so the outcome depends on initialization.

Worked Example 2

Problem. Why do weight-permutation symmetries create many equivalent minima in a neural net?

  1. A hidden layer's units can be reordered without changing the function the net computes.
  2. Permuting unit i with unit j (and their weights) gives different parameter vectors w but identical loss.
  3. For h hidden units there are h! such permutations, each a distinct global-equivalent minimum.
  4. So the 'global minimum' is really a large set of symmetric copies — one reason the landscape is non-unique.

Answer. Reordering hidden units leaves the network function (and loss) unchanged, so each minimum has h! symmetric copies, making neural-net minima inherently non-unique.

Common mistakes
  • Expecting a unique optimum for a neural network. Symmetries and nonconvexity give many equivalent or distinct minima; 'the' solution is not unique.
  • Assuming all local minima in deep nets are bad. In high dimensions most low-loss critical points are saddles or near-global minima; getting stuck in a terrible local min is rarer than naive intuition suggests.
✎ Try it yourself

Problem. You train the same network twice from two random seeds and get two different weight vectors with nearly equal training loss. Explain why this is expected and consistent with the loss-landscape picture.

Solution. Neural-network loss landscapes are nonconvex and riddled with symmetries (e.g. permuting hidden units, or sign flips with odd activations) plus many distinct low-loss basins. Different random initializations start the optimizer in different basins of attraction, so SGD converges to different parameter vectors. Because many of these minima have comparably low loss (a known empirical property of overparameterized nets), the two runs reach nearly equal training loss despite different weights. This is fully consistent with the landscape view: there is no single global optimum but a large set of good solutions, and initialization selects which one you find.

Regularization (L1, L2) and the bias-variance tradeoff

Regularization adds a penalty that trades a little bias for a lot less variance — the lever that controls overfitting.

Regularization adds a penalty to the loss to discourage overly complex models: ridge (L2) adds lambda||w||^2, LASSO (L1) adds lambda||w||_1. This implements the bias-variance tradeoff: an unregularized flexible model has low bias but high variance (it overfits noise), while a heavily regularized model has higher bias but lower variance. The regularization strength lambda tunes this balance — larger lambda shrinks weights more, raising bias and lowering variance. L2 shrinks all coefficients smoothly and stabilizes ill-conditioned problems (it adds lambda to the Hessian's eigenvalues, improving conditioning). L1 drives some coefficients to exactly zero, performing feature selection and yielding sparse, interpretable models. From the optimization view, both keep the convex problem well-posed; choosing lambda by cross-validation is how you minimize the expected (test) error rather than just training error.

Worked Example 1

Problem. Show how ridge (L2) improves conditioning: compare the Hessian of least squares vs. ridge.

  1. Least-squares Hessian: 2 X^T X, eigenvalues can be near zero if features are collinear (ill-conditioned).
  2. Ridge objective ||Xw-y||^2 + lambda||w||^2 has Hessian 2(X^T X + lambda I).
  3. Adding lambda I shifts every eigenvalue up by lambda, so the smallest eigenvalue >= 2 lambda > 0.
  4. Condition number drops, so gradient descent converges faster and the solution is stable.

Answer. Ridge adds lambda I to X^T X, lifting all eigenvalues by lambda; this guarantees positive definiteness and a smaller condition number, stabilizing the fit.

Worked Example 2

Problem. As lambda increases from 0 to large, describe what happens to bias, variance, and sparsity (for L1).

  1. lambda = 0: no penalty, lowest bias but highest variance (overfits).
  2. Moderate lambda: weights shrink, variance falls faster than bias rises -> better test error.
  3. Large lambda: weights heavily shrunk toward 0, high bias, low variance (underfits).
  4. For L1 specifically, increasing lambda zeros out more coefficients, increasing sparsity.

Answer. Increasing lambda raises bias and lowers variance; an intermediate lambda minimizes test error, and for L1 larger lambda yields more exact zeros (more sparsity).

Worked Example 3

Problem. Select lambda by cross-validation in scikit-learn for ridge regression (sketch the code).

  1. `from sklearn.linear_model import RidgeCV; import numpy as np`
  2. `alphas = np.logspace(-3, 3, 13)` # candidate lambda values
  3. `model = RidgeCV(alphas=alphas, cv=5).fit(X, y)`
  4. `model.alpha_` -> the lambda minimizing cross-validated error; this targets test error, not training error.

Answer. RidgeCV searches alphas via 5-fold CV and picks model.alpha_ that minimizes validation error — the principled way to set the bias-variance lever.

Common mistakes
  • Choosing lambda to minimize training error. Training error always favors lambda=0; tune lambda on a validation set or by cross-validation to control overfitting.
  • Using L1 expecting smooth shrinkage or L2 expecting sparsity. L1 gives exact zeros (sparsity); L2 gives smooth shrinkage with no zeros — pick by whether you want feature selection.
✎ Try it yourself

Problem. A model has near-perfect training accuracy but poor test accuracy. Explain in bias-variance terms and say how L2 regularization helps and how you would choose its strength.

Solution. Near-perfect training accuracy with poor test accuracy is the signature of high variance / overfitting: the model has fit noise in the training data (low bias, high variance). L2 regularization adds lambda||w||^2, shrinking the weights toward zero, which reduces the model's effective complexity and variance at the cost of a small increase in bias. From the optimization side it also adds lambda I to the Hessian, improving conditioning and stability. To choose lambda, use k-fold cross-validation (e.g. RidgeCV over a logspace grid of alphas) and pick the value that minimizes validation error — the lambda that best balances bias and variance for generalization, not the one that minimizes training error (which would be lambda=0).

Saddle points, local minima, and overparameterization

In high-dimensional deep nets the real obstacle is saddle points, not bad local minima — and overparameterization surprisingly helps.

In nonconvex training the gradient vanishes not only at minima but at saddle points, where the Hessian has both positive and negative eigenvalues. In high dimensions saddles vastly outnumber local minima, and they slow training because the gradient is tiny near them — but they are escapable: a direction of negative curvature exists, and SGD's noise plus momentum push the iterate out. Bad local minima (much worse than global) turn out to be rare in large nets. Overparameterization — using more parameters than data points — reshapes the landscape so that almost all minima achieve near-zero training loss and are connected, making optimization easier rather than harder, and (with implicit regularization from SGD) often generalizing well. Practically this means: don't fear saddles unduly, use stochastic methods with momentum, and large models can be easier to train.

Worked Example 1

Problem. Classify the critical point (0,0) of f(x,y) = x^2 - y^2 and explain why gradient descent stalls there.

  1. Gradient (2x, -2y) = 0 at (0,0); it is a critical point.
  2. Hessian [[2,0],[0,-2]] has eigenvalues +2 and -2 -> indefinite -> saddle point.
  3. Along x the surface curves up (min), along y it curves down (max).
  4. Plain gradient descent approaching along the x-axis sees a vanishing gradient and stalls, until noise/curvature nudges it along the negative-curvature y direction to escape.

Answer. (0,0) is a saddle (Hessian eigenvalues +2, -2); the small gradient near it stalls plain GD, but a negative-curvature escape direction exists, which SGD noise/momentum exploit.

Worked Example 2

Problem. Why does overparameterization tend to make training loss reach (near) zero?

  1. With more parameters than constraints (data points), the system w that fits all examples is underdetermined.
  2. There exist many parameter settings achieving zero (or near-zero) training loss — a high-dimensional solution manifold.
  3. Gradient descent only needs to reach any point on that manifold, which is easy when it is large and connected.
  4. SGD's path implicitly selects a low-norm / smooth solution among them, often generalizing well.

Answer. Overparameterization creates a large, connected set of zero-loss solutions, so descent easily reaches one; implicit regularization from SGD then favors a generalizing member of that set.

Common mistakes
  • Confusing a saddle point with a local minimum. Both have zero gradient; only the Hessian (mixed-sign eigenvalues = saddle) distinguishes them.
  • Assuming more parameters always overfit. In the overparameterized regime with SGD's implicit regularization, larger models can train more easily and still generalize (the 'double descent' phenomenon).
✎ Try it yourself

Problem. At a critical point of a loss, the Hessian eigenvalues are {3, 0.5, -2, 4}. Is it a minimum, maximum, or saddle, and what does it imply for training?

Solution. The Hessian has both positive eigenvalues (3, 0.5, 4) and a negative one (-2), so it is indefinite -> the point is a saddle point, not a minimum or maximum. The negative eigenvalue -2 indicates a direction of negative curvature along which the loss decreases. For training, the gradient is zero here so plain gradient descent can stall and waste iterations near the saddle, but because an escape direction exists, stochastic gradients (SGD noise) and momentum will perturb the iterate along the negative-curvature direction and let it continue descending. This is why saddles, though they slow training, are not fatal in practice.

Practical training: schedules, early stopping, and diagnostics

Beyond the update rule, practical training is about controlling the learning rate over time and knowing when to stop — read the curves.

Practical training wraps the optimizer in control logic. A learning-rate schedule (step decay, cosine annealing, warmup) takes big early steps then refines, which is essential for SGD's noise to settle. Early stopping monitors validation loss and halts when it stops improving (with a patience window), acting as implicit regularization that prevents overfitting without changing the loss. Diagnostics are how you steer: plot training and validation loss versus epoch — a rising/oscillating curve means the learning rate is too high; a flat curve means it is too low or gradients vanish; a widening train-vs-validation gap means overfitting; both still falling means keep training. Gradient-norm tracking catches vanishing/exploding gradients. Together, schedules + early stopping + curve reading turn the raw optimizer into a reliable training recipe.

Worked Example 1

Problem. Validation loss: epochs [1.0, 0.7, 0.55, 0.54, 0.56, 0.58] with early-stopping patience 2. When does training stop and which weights are kept?

  1. Best validation loss so far improves at epochs 1->4 (0.54 at epoch 4 is the minimum).
  2. Epoch 5: 0.56 > 0.54 (no improvement, patience counter = 1).
  3. Epoch 6: 0.58 > 0.54 (no improvement, patience counter = 2 -> patience exhausted).
  4. Stop after epoch 6; restore the weights from epoch 4 (the best validation loss).

Answer. Training stops after epoch 6 (2 epochs without improvement); the kept weights are from epoch 4, where validation loss was lowest (0.54).

Worked Example 2

Problem. Diagnose: training loss keeps falling but validation loss has been rising for several epochs. What is happening and what do you do?

  1. Falling train loss + rising validation loss = a widening generalization gap = overfitting.
  2. The model is memorizing training noise rather than learning generalizable structure.
  3. Remedies: stop early (use the best-validation checkpoint), increase regularization (L2/dropout), or get more data.
  4. Confirm by plotting both curves; the divergence point is roughly where to stop.

Answer. It is overfitting (generalization gap widening); apply early stopping at the validation minimum and/or add regularization rather than training further.

Worked Example 3

Problem. Cosine schedule with warmup: linear warmup for 5 steps to eta_max=0.1, then cosine decay over 95 steps. Give eta at steps 0, 5, 100.

  1. Warmup: eta(step) = 0.1 * step/5 for step <= 5. Step 0: 0.1*0/5 = 0.
  2. End of warmup, step 5: eta = 0.1 (the peak).
  3. Cosine decay from step 5 to 100 (span 95): eta = 0.5*0.1*(1 + cos(pi*(step-5)/95)).
  4. Step 100: 0.5*0.1*(1 + cos(pi)) = 0.05*(1 - 1) = 0.

Answer. eta = 0 at step 0 (warmup start), 0.1 at step 5 (peak), and 0 at step 100 (end of cosine decay); warmup stabilizes early training, cosine cools it to zero.

Common mistakes
  • Monitoring only training loss for early stopping. Training loss almost always keeps falling; stop based on validation loss to catch overfitting.
  • Forgetting to restore best weights after early stopping. Keep a checkpoint of the best-validation model; the final-epoch weights may be worse.
✎ Try it yourself

Problem. Validation loss over epochs is [0.9, 0.6, 0.5, 0.52, 0.49, 0.51, 0.53] with early-stopping patience 2. State when training stops and which epoch's weights are restored.

Solution. Track the best (lowest) validation loss and a no-improvement counter. Epoch 1:0.9 (best). Epoch 2:0.6 (best). Epoch 3:0.5 (best). Epoch 4:0.52 (>0.5, counter=1). Epoch 5:0.49 (new best, counter resets to 0). Epoch 6:0.51 (>0.49, counter=1). Epoch 7:0.53 (>0.49, counter=2 -> patience exhausted). Training stops after epoch 7, and the restored weights are from epoch 5, where validation loss was the minimum (0.49). The dip at epoch 5 mattered: a naive stop after the epoch-4 uptick would have missed the better epoch-5 model, which is why patience > 0 is used.

Putting it together: train and tune a model

The capstone: combine ERM, a loss, an optimizer, regularization, and diagnostics into one end-to-end training-and-tuning workflow.

End-to-end model training ties the whole course together. The workflow: (1) frame the task as ERM — pick a model and a loss; (2) split data into train/validation/test; (3) choose an optimizer (convex solver or SGD/Adam) and a learning-rate schedule; (4) add regularization (L1/L2) and tune its strength plus the learning rate via a validation set or cross-validation; (5) use early stopping and read loss curves to control training; (6) report final performance on the held-out test set only once. Each step draws on earlier units: convexity tells you whether the optimum is guaranteed, gradient methods do the fitting, regularization controls the bias-variance tradeoff, and diagnostics keep training honest. This disciplined loop is how optimization theory becomes a working, well-generalizing model.

Worked Example 1

Problem. Outline an end-to-end logistic-regression training pipeline with hyperparameter tuning in scikit-learn.

  1. Split: `from sklearn.model_selection import train_test_split; Xtr,Xte,ytr,yte = train_test_split(X,y,test_size=0.2)`.
  2. Model + regularization: logistic regression with L2 penalty, strength C = 1/lambda.
  3. Tune C by CV: `from sklearn.linear_model import LogisticRegressionCV; m = LogisticRegressionCV(Cs=10, cv=5).fit(Xtr, ytr)`.
  4. Evaluate once on the held-out test set: `m.score(Xte, yte)`.

Answer. Pipeline: split data, fit L2-regularized logistic regression with C chosen by 5-fold CV on the training set, then report accuracy on the untouched test set — a complete tune-and-evaluate loop.

Worked Example 2

Problem. Why must the test set be touched only once, at the very end?

  1. Hyperparameters (lambda, learning rate, architecture) are chosen using validation/CV performance.
  2. If you tune on the test set, its performance becomes optimistic — you have leaked test information into model selection.
  3. The test set's job is an unbiased estimate of generalization on truly unseen data.
  4. So evaluate on it exactly once after all tuning is final; otherwise the reported score overstates real-world performance.

Answer. Tuning against the test set leaks information and inflates the score; the test set must stay untouched until the final, single evaluation to give an honest estimate of generalization.

Worked Example 3

Problem. Sketch a from-scratch training loop combining mini-batch SGD, L2 regularization, a decaying learning rate, and early stopping.

  1. `for epoch in range(max_epochs):` shuffle data, iterate mini-batches.
  2. `grad = batch_grad(w) + 2*lam*w; eta = eta0/(1 + decay*epoch); w -= eta*grad` # L2 + decay.
  3. After each epoch compute validation loss; track the best and a no-improvement counter.
  4. `if counter >= patience: w = best_w; break` # early stop and restore best weights.

Answer. The loop fuses mini-batch SGD, an L2 gradient term, a per-epoch learning-rate decay, and validation-based early stopping with best-weight restore — the full practical training recipe in one place.

Common mistakes
  • Tuning hyperparameters on the test set. This leaks information and inflates reported accuracy; tune on validation/CV and reserve the test set for one final evaluation.
  • Reporting only training accuracy. Always report validation/test metrics, since training accuracy alone hides overfitting and says nothing about generalization.
✎ Try it yourself

Problem. You must train and tune a regularized linear classifier. Describe the full workflow from raw data to a reported generalization estimate, naming the role of each optimization concept from this course.

Solution. 1) Frame as ERM: choose a convex loss (logistic/hinge) and a linear model, so by convexity the optimizer reaches a global minimum. 2) Split data into train/validation/test. 3) Standardize features so the loss is well-conditioned (low condition number -> fast, stable gradient descent). 4) Add an L2 (or L1) regularizer to manage the bias-variance tradeoff; this also makes the Hessian PD (well-posed convex problem). 5) Fit with gradient descent / SGD / Adam using a learning-rate schedule (warmup + decay) and monitor training and validation loss curves. 6) Tune the regularization strength lambda and learning rate via k-fold cross-validation, picking values that minimize validation error. 7) Use early stopping at the validation minimum to prevent overfitting. 8) Finally, evaluate the chosen model once on the untouched test set to report an unbiased generalization estimate. Every step maps to a course concept: convexity (guaranteed optimum), conditioning (feature scaling), gradient methods (fitting), regularization (bias-variance), and diagnostics (curves, early stopping).

Key terms
  • Empirical risk minimization (ERM) — choosing parameters to minimize average loss over the training data.
  • Loss function — the per-example penalty (e.g. squared error, cross-entropy) whose average we minimize.
  • Loss landscape — the surface of the loss as a function of parameters; convex for linear models, nonconvex for deep nets.
  • Regularization — adding a penalty (L1, L2) to the loss to control complexity and improve generalization.
  • L2 (ridge) penalty — adds lambda * ||w||^2, shrinking weights and keeping the objective smooth and strongly convex.
  • L1 (lasso) penalty — adds lambda * ||w||_1, encouraging sparse solutions with many zero weights.
  • Saddle point — a stationary point that is a minimum in some directions and a maximum in others, common in deep nets.
  • Bias-variance tradeoff — the balance where more regularization raises bias but lowers variance, tuned for best generalization.
Assignment · Optimization behind a real model

Train a regularized logistic-regression or linear model on a small dataset by minimizing its loss with gradient-based optimization (your own loop or scikit-learn/CVXPY). Sweep the L2 (and optionally L1) regularization strength, plotting training vs. validation loss to expose the bias-variance tradeoff and any sparsity from L1. Then take a tiny neural net and visualize a 1D or 2D slice of its loss landscape to contrast its nonconvexity with the convex linear model.

Deliverable · A notebook with the regularization-path and train/validation curves, an L1 sparsity comparison, a loss-landscape slice for the neural net, and a paragraph linking regularization to optimization and generalization.

Quiz · 5 questions
  1. 1. Empirical risk minimization trains a model by minimizing:

  2. 2. Compared with linear/logistic regression, a deep neural network's loss landscape is:

  3. 3. Adding an L1 penalty to the loss tends to:

  4. 4. Increasing regularization strength generally:

  5. 5. Why are saddle points relevant when training deep networks?

You'll be able to

I can frame supervised learning as empirical risk minimization with a chosen loss and regularizer.

I can explain why convex-model loss surfaces differ from neural-network loss landscapes.

I can apply L1/L2 regularization and connect it to the bias-variance tradeoff and optimization behavior.

Where this leads

Course milestones

Convex modeling checkpoint: take a word problem, write it in standard form, and prove whether it is a convex program (Unit 1).
Solver foundations: solve an unconstrained problem with gradient/Newton methods and verify against scipy.optimize (Unit 2).
Optimizer build-off: implement gradient descent, SGD, momentum, and Adam from scratch and compare convergence on a real loss (Unit 3).
Constrained-and-duality lab: derive and verify the KKT conditions by hand, then reproduce the solution and duals in CVXPY (Units 4-6).
Capstone — optimization behind a model: train a regularized model via ERM, expose the bias-variance tradeoff, and contrast convex vs. neural-net loss landscapes (Unit 7).

Free, forever. Math you can actually use.

Part of the open-source Crunch Academy tech-education path. Pair it with the data-science and ML tracks when you're ready.

All Math for Data Science courses Crunch Academy home