Convex Optimization School

How to use this course

1. One day per sitting. Read the lesson from top to bottom. Each day builds on the one before.
2. Play with every widget. Each one is a small optimization running in your browser, not a picture.
3. Take the quiz, then mark the day done. Progress is saved on this device only.
4. Use the other tabs. Algorithms has one card per method. Papers explains the robotics papers from Week 3. Glossary lets you look up any term.
5. Open the sources. The small chips after sentences (like this one CO4A) link to the original.

What you need to know first

Some high-school algebra, what a derivative is, and the basics of vectors and matrices: adding them, the dot product, and multiplying a matrix by a vector. You do not need to read proofs. The few you will meet are optional and sit in grey "Going deeper" boxes. Keep this notation table open while you read:

The map

What the coloured boxes mean

About the sources

Weeks 1 and 2 follow the Korean open textbook 모두를 위한 컨벡스 최적화 (Convex Optimization for All, "CO4A"), released under the MIT License (© 2021 CONVEX-OPTIMIZATION-FOR-ALL Team). CO4A is based on CMU's 10-725 course (Fall 2016, taught by Javier Peña and Ryan Tibshirani) and on Boyd & Vandenberghe's Convex Optimization (2004). This course rephrases and simplifies that material in English. It reorders some topics, moves proofs into optional boxes, and adds widgets and worked numbers. Week 3 draws on the original papers and on official solver documentation. Every lesson links the CO4A chapters it covers. All sources were checked on 2026-10-26, and fast-moving facts are marked "as of Oct 2026".

1. Knobs, a score, and rules

Every optimization problem has the same three parts. Think of planning a drive home. You choose the turns, you want the trip to be short, and you must obey the one-way streets.

The things you choose are the optimization variables. We stack them into one vector x = (x₁, …, x_n) in ℝⁿ. The score you want small is the objective function f(x), also called the cost. The rules are the constraints.

A point that obeys every rule is called feasible. The set of all feasible points is the feasible set. Written out, every problem in this course has the same shape. "s.t." is short for "subject to". CO4A 1-1 CMU 10-725 2016

minimize f(x)
s.t. g_i(x) ≤ 0, i = 1, …, m
h_j(x) = 0, j = 1, …, r

Each g_i is an inequality constraint ("at most"). Each h_j is an equality constraint ("exactly"). Any rule can be moved into this form. For example, x ≤ 1 becomes g(x) = x − 1 ≤ 0.

The best score you can reach is the optimal value, written p*. A feasible point that reaches it is an optimal solution, or minimizer, written x*. Keep the two apart: p* is a number, and x* is a point.

Here is a problem with one knob. Minimize (x − 2)² subject to x ≤ 1. Without the rule, the best choice would be x = 2, with a score of 0. The rule forbids it, so the best allowed point sits on the edge: x* = 1, with p* = (1 − 2)² = 1.

Some rules are not written down at all. They hide inside the formula. Take the problem: minimize f(x) = x − log x. The logarithm only exists for x > 0, so x > 0 is a hidden rule, called an implicit constraint. Setting the slope f′(x) = 1 − 1/x to zero gives x* = 1 and p* = 1 − log 1 = 1. CO4A 1-1

A problem with no written-down constraints is an unconstrained problem. The example above counts as unconstrained, even though it hides the rule x > 0.

In practice you describe the three parts to a modelling tool and let a solver do the rest. Here is the one-knob problem in CVXPY, a free Python library. Day 4 explains how it checks that a problem is convex before it solves it.

import cvxpy as cp

x = cp.Variable()                          # the knob
objective = cp.Minimize(cp.square(x - 2))  # the score
constraints = [x <= 1]                     # the rule
problem = cp.Problem(objective, constraints)
problem.solve()
print(x.value, problem.value)              # about 1.0 and 1.0

2. Optimization is everywhere

Once you look for knobs, a score and rules, you find optimization problems everywhere. The table lists some from daily life, driving, robotics and machine learning. The last column previews today's main question: is the problem convex? Sections 3 and 4 explain how to tell.

Fitting a line. Suppose you measure three points: (0, 1), (1, 3) and (2, 4). You want the line y = ax + b that passes closest to them. The score adds up the squared vertical misses:

f(a, b) = (b − 1)² + (a + b − 3)² + (2a + b − 4)²

This score is a bowl over the two knobs a and b. Its lowest point is a = 1.5 and b ≈ 1.17, where the total squared miss is 0.167. Every "least squares" fit in statistics and machine learning is a bigger version of this bowl. Day 6 shows gradient descent finding its bottom.

Planning a car's speed. A self-driving car re-plans its speed many times a second. The knobs are its speeds over the next few seconds. The score adds up how hard it brakes and accelerates, so the ride stays smooth. The rules keep it under the speed limit and a safe distance behind the car ahead.

Distance travelled is a sum of speed × time, so "stay 10 m behind the car ahead" is a linear rule in the speeds. With linear rules and a squared score, the problem is convex. It is called a quadratic program, and Day 17 builds one step by step. Boyd and Vandenberghe's course slides say that convex solvers generated as code now run in milliseconds, and cite the landing of SpaceX's Falcon 9 rocket. Boyd EE364a slides

A robot arm, which is not convex. Picture an arm with two links, each 1 m long, trying to touch the point (1, 1). Two poses reach it exactly. "Elbow down" puts the shoulder at 0° and bends the elbow 90°. "Elbow up" puts the shoulder at 90° and bends the elbow −90°.

Now blend the two answers. The halfway pose has the shoulder at 45° and a straight elbow (0°). It puts the hand at (1.41, 1.41), which misses the target by 0.59 m. Two perfect answers, and the point halfway between them is bad. Section 4 shows that this cannot happen in a convex problem.

Machine learning. Training a model means choosing its weights to make the average error on the data small. For linear regression and logistic regression this problem is convex. Boyd & Vandenberghe 2004 For a deep neural network it is not, which is one reason training one needs careful starting points and tuning. Day 16 and Day 20 come back to this.

3. Convex sets and convex functions, one picture each

Two shapes carry the whole subject. Today you only need a rough picture of each. Day 2 and Day 3 treat them properly.

A set with no dents. Pick any two points inside a shape and draw the straight path between them. If the path never leaves the shape, whichever two points you pick, the shape is a convex set. A disk passes. A crescent moon fails, because the path between its two tips cuts across the empty bite.

To write "every point on the path", blend the two points x and y with a weight θ between 0 and 1. θ = 1 gives x, θ = 0 gives y, and θ = 0.5 gives the midpoint. Such a weighted average is called a convex combination. A set C is convex when CO4A 1-2

x, y ∈ C and θ ∈ [0, 1] ⇒ θx + (1 − θ)y ∈ C

A number check. The halfspace of points (p, q) with p + 2q ≤ 4 contains x = (0, 0) and y = (4, 0). With θ = 0.25 the blend is 0.25·(0, 0) + 0.75·(4, 0) = (3, 0). Since 3 + 2·0 = 3 ≤ 4, the blend is inside too.

A function shaped like a bowl. Stretch a string between any two points on the graph of a function. If the string never dips below the graph, the function is a convex function. On a bumpy curve, some strings pass under a hump.

f(θx + (1 − θ)y) ≤ θf(x) + (1 − θ)f(y)

The left side is the height of the graph at the blended point. The right side is the height of the string there. Take f(x) = x², x = −1, y = 3 and θ = 0.5. The blended point is 1, where the graph has height f(1) = 1. The string has height (1 + 9)/2 = 5 there. Since 1 ≤ 5, the string is above the graph, by 4.

The two ideas are linked. Paint everything on or above the graph of f. That painted region is the epigraph of f. A function is convex exactly when its epigraph is a convex set. CO4A 1-2 Boyd & Vandenberghe 2004

Check it on f(x) = |x|. The points (−1, 1) and (1, 1) lie on its graph. Their midpoint (0, 1) has height 1, which is above |0| = 0, so it lies inside the painted region, as convexity requires.

4. What makes a problem convex

Now put the pieces together. A convex optimization problem needs three things. The objective f is a convex function. Each inequality function g_i is convex. Each equality function h_j is affine, meaning flat: h_j(x) = a_j^Tx + b_j, a linear part plus a constant. CO4A 1-2 CMU 10-725 2016

Why these rules? A convex g_i is a bowl, and the region where it is at most 0 is the part of the bowl under the water line: a convex blob. A flat h_j makes the region where it equals 0 a flat line or plane. Overlapping convex pieces gives a convex feasible set (Day 2 shows why). A bowl-shaped score over a convex feasible set is the whole recipe.

Example A, convex. Minimize x² subject to 1 − x ≤ 0, which means x ≥ 1. The objective is a bowl and the rule is flat. The feasible set is the ray [1, ∞), with no gaps. The answer is x* = 1 with p* = 1.

Example B, not convex. Minimize x² subject to x² − 1 = 0. The equality is curved, not flat. Its feasible set is just two points, −1 and 1. The path between them runs through 0, which is not allowed.

Example C, not convex. Minimize x² subject to x² ≥ 1, written as 1 − x² ≤ 0. Here g(x) = 1 − x² is an upside-down bowl. Its feasible set, everything outside (−1, 1), has a gap in the middle.

Two more facts. Maximizing a concave function, an upside-down bowl, is also a convex problem, because it is the same as minimizing its negative. And modelling tools check the rules before solving. CVXPY raises an error if you ask it to solve example B, because "x² == 1" is not an affine equality. CVXPY docs

The robot arm from section 2 fails for the same kind of reason. In a convex problem, the set of best points is itself convex, so the blend of two best answers is also a best answer. The arm's two perfect poses blend into a pose that misses by 0.59 m, so its problem cannot be convex.

5. The payoff: every local minimum is global

Imagine a hiker in thick fog who wants the lowest point of a valley. She can only feel the ground under her boots. She steps downhill until every direction goes up, and then she stops.

Where she stops is a local minimum: no nearby point is lower. What she wants is a global minimum: no point anywhere is lower. In mountains these can differ, and she may stop in a small hollow high on a slope.

On a convex bowl they cannot differ. There is only one valley, so "nothing nearby is lower" already means "nothing anywhere is lower". This is the most useful fact in the course: CMU 10-725 2016 Boyd & Vandenberghe 2004

if f(x) ≤ f(y) for every feasible y with ‖y − x‖₂ ≤ ρ,
then f(x) ≤ f(y) for every feasible y

Here x is a feasible point, ρ > 0 is the radius of "nearby", and ‖y − x‖₂ is the ordinary straight-line distance. The reason fits in one line. If a better point y existed far away, a short step from x toward y would already be better, because the string to y lies above the bowl. The box at the end of this section works it through with numbers.

For a smooth function with no constraints, the test gets even simpler. If f is convex and differentiable, then ∇f(x) = 0 holds exactly at the global minimizers. CMU 10-725 2016 For f(x) = (x − 1)², setting f′(x) = 2(x − 1) to zero gives x = 1, and that is the bottom of the bowl.

The simplest computer version of the hiker is gradient descent. It repeats one move: step against the slope, scaled by a step size η (the Greek letter eta). See the gradient descent card for more.

x_k+1 = x_k − η f′(x_k)

Run it on f(x) = (x − 1)² with η = 0.1, starting at x₀ = 2. The slope there is 2, so x₁ = 2 − 0.1·2 = 1.8. Each step shrinks the distance to the bottom by a factor of 0.8. After 21 steps the ball is within 0.01 of x = 1.

Try it: tap either curve to start both balls there and watch them roll, or set x₀ with the slider and press Roll. Both balls use the same rule from the same start. Then move x₀ left of 0.13 and roll again. Finally, push the bowl's step size above 1.

Notice: from the default start x₀ = 2.00, the bowl's ball settles at x = 1.000 after 66 steps. The bumpy ball settles after 48 steps at x = 1.347, where g = −2.619. That is only a local minimum. The global minimum, g = −5.444 at x = −1.473, lies beyond the hilltop at x ≈ 0.126, and the ball never sees it.

6. Why convex problems are the "easy" ones

Most algorithms are like the hiker. They only look at the ground near the current point. On a nonconvex problem that local view can mislead them, as the bumpy ball showed.

Boyd and Vandenberghe put it bluntly. For general nonlinear problems there is no effective general method, and even some problems with about ten variables are very hard. Boyd & Vandenberghe 2004

The usual fallback is local optimization: start somewhere, go downhill, and hope. It needs a good starting guess and some tuning. It also cannot tell you how far you are from the true best. The book calls it "more art than technology". Global search does find the true best, but its worst-case cost grows exponentially with the number of variables.

Convex problems escape this trap. Usually there is no formula for the answer. But interior-point methods solve them reliably, almost always in 10 to 100 iterations. In the book's words, once you have formulated a problem as convex, solving it is "(almost) technology". Boyd & Vandenberghe 2004 Day 12 opens up these methods.

So the hard part moves from solving to formulating. Your job is to recognise, or rewrite, a problem as convex. That skill is what Week 1 trains.

This also corrects an old rule of thumb, "linear is easy, nonlinear is hard". The better dividing line is convex versus nonconvex. A curved bowl is just as easy as a flat plane, and only the bumps make problems hard. Boyd EE364a slides CMU 10-725 2016

Convex does not mean any method will do, though. Tibshirani's first lecture cleans up a noisy image by solving one convex problem with three methods. A specialised ADMM needed 20 iterations, proximal gradient needed 1,000, and coordinate descent needed 10,000 cycles. CMU 10-725 2016 CO4A 1-3

The lesson is not "ADMM is best". ADMM won here because its steps matched the problem's structure. On another problem a different method wins. Choosing the method is the theme of Week 2.

The ideas are older than computers, and the table gives the main steps. CO4A 1-4 Boyd EE364a slides

7. The road ahead

The course has three weeks. Week 1 teaches you to recognise convex problems: convex sets (Day 2), convex functions (Day 3), convex problems (Day 4), and the standard families from linear programs up to semidefinite programs (Day 5). Week 2 opens up the solvers. You will meet gradient descent, subgradients, proximal methods, duality, the KKT conditions, Newton's method, interior-point methods, ADMM, Frank–Wolfe and coordinate descent (Days 6–14). Week 3 leaves the safe bowl. It covers yes-or-no decisions, general nonlinear problems, and robots and cars: model predictive control, trajectory optimization, collision avoidance, safety filters, and solvers in practice (Days 15–21).

1. Blending two points: lines, segments and rays

Picture a slider between two points x and y. At one end you sit on x, at the other on y. Let the slider run past both ends and you trace the whole line through them.

The slider position is a number θ. The point it gives is a blend of x and y. CO4A 2-1-1

z = θx + (1 − θ)y = y + θ(x − y)

Read the second form as a recipe: start at y and move θ of the way toward x. θ = 0 gives y, θ = 1 gives x, and θ = 0.5 gives the midpoint.

The range of θ decides the shape. If θ can be any real number, you get the line through x and y. If 0 ≤ θ ≤ 1, you get the line segment between them. If θ ≥ 0, you get a ray: it starts at y, passes through x, and goes on forever.

A number check with x = (1, 0) and y = (0, 1). θ = 0.5 gives 0.5·(1, 0) + 0.5·(0, 1) = (0.5, 0.5), the midpoint, on the segment. θ = 3 gives 3·(1, 0) − 2·(0, 1) = (3, −2). That point is on the line, far past x, but not on the segment.

2. Affine sets and convex sets

Some sets are flat and endless, like a point, a line or a plane. An affine set contains the whole line through any two of its points: the blend θx + (1 − θ)y stays inside for every real θ. The solutions of linear equations, {x : Ax = b}, are the typical example. CO4A 2-1-2

Blends of more than two points work the same way. A sum θ₁x₁ + … + θ_kx_k whose weights add up to 1 is an affine combination. The weights may be negative. An affine set contains every affine combination of its points.

Take C = {x ∈ ℝ² : x₁ + x₂ = 1}, a line. It contains (1, 0) and (0, 1), and also their blend (3, −2) from above, since 3 − 2 = 1. Slide C by (1, 0) toward the origin and you get the line x₁ + x₂ = 0 through the origin. Every affine set is a subspace slid away from the origin like this.

A convex set asks for less. Only the segment has to stay inside, not the whole line. This is the test from Day 1: no dents and no holes. CO4A 2-1-3 Boyd & Vandenberghe 2004

x, y ∈ C and θ ∈ [0, 1] ⇒ θx + (1 − θ)y ∈ C

A weighted average θ₁x₁ + … + θ_kx_k with every θ_i ≥ 0 and weights adding up to 1 is a convex combination. A convex set contains every convex combination of its points, not only blends of two at a time.

Stretch a rubber band around a set of points and let it snap tight. The region inside the band is the convex hull, the smallest convex set that contains them. It is also the set of all their convex combinations.

For example, the convex hull of (0, 0), (2, 0), (0, 2) and (1, 1) is the triangle with corners (0, 0), (2, 0) and (0, 2). The fourth point adds nothing. It is 0.5·(2, 0) + 0.5·(0, 2), so it already lies on an edge.

Every affine set is convex, because the segment is part of the line. The reverse fails: a disk is convex, but the line through two of its points runs out of the disk. The edge matters too. A square that keeps its top corners but not the middle of its top edge is not convex.

Try it: pick a shape, then drag P and Q (or tap to move the nearer one). The segment is solid green where it stays inside and dashed red where it leaves. "Find a bad pair" lets the computer search 2,000 random pairs for a leak.

Notice: in the default crescent both handles are inside, P = (−0.20, 0.90) and Q = (−0.20, −0.90). Yet the points θP + (1 − θ)Q leave the set for θ from 0.285 to 0.715, where the segment crosses the bite. One leaking segment proves the crescent is not convex. On the disk, square or cut disk, no pair you drag ever leaks.

3. Cones: rays from the origin

Think of a flashlight beam with the bulb at the origin. If a point is lit, so is every point farther out along the same ray. A set C is a cone if x ∈ C and θ ≥ 0 always give θx ∈ C. CO4A 2-1-4

A convex cone also contains all nonnegative sums: x, y ∈ C and θ₁, θ₂ ≥ 0 give θ₁x + θ₂y ∈ C. In the plane it looks like a pie slice with its tip at the origin, running on forever.

The quarter-plane ℝ²₊, where both coordinates are ≥ 0, is a convex cone. From (1, 2) and (3, 0), the sum 2·(1, 2) + 1·(3, 0) = (5, 4) is still in it.

A cone does not have to be convex. The two coordinate axes together form a cone, since stretching a point on an axis keeps it there. But (1, 0) + (0, 1) = (1, 1) is on neither axis, so this cone is not convex. CMU 10-725 2016

4. A gallery of basic convex sets

A few shapes appear in almost every problem. Learn to spot them the way you spot a circle. The figure shows them, and the paragraphs below give formulas and number checks. CO4A 2-2-1 CMU 10-725 2016

Hyperplanes and halfspaces. A hyperplane {x : a^Tx = b}, with a ≠ 0, is a flat wall. The vector a, called the normal, sticks straight out of the wall. A halfspace {x : a^Tx ≤ b} is one side of the wall, wall included. It is convex but not affine. Take a = (1, 2) and b = 4. The point (2, 1) gives 2 + 2 = 4, so it is on the wall. (0, 0) gives 0 ≤ 4, inside. (4, 4) gives 12 > 4, outside.

Balls. The Euclidean ball with centre c and radius r is every point within distance r of c: {x : ‖x − c‖₂ ≤ r}. With c = (1, 1) and r = 2, the point (2, 2) is √2 ≈ 1.414 away from the centre, so it is inside.

Ellipsoids. Stretch a ball along some directions and you get an ellipsoid: {x : (x − c)^TP⁻¹(x − c) ≤ 1}. Here P is a symmetric positive definite matrix (all eigenvalues > 0). The semi-axis lengths are the square roots of P's eigenvalues. With P = diag(4, 1) and c = 0, the set is x₁²/4 + x₂² ≤ 1, with semi-axes 2 and 1. A ball is the special case P = r²I. Boyd & Vandenberghe 2004

Norm balls. You can measure length in other ways than the straight-line distance. The p-norm is

‖x‖_p = (|x₁|^p + … + |x_n|^p)^1/p

The absolute values matter: without them, negative entries would break the formula. The norm ball ‖x‖_p ≤ 1 is a diamond for p = 1, a disk for p = 2, and a square for p = ∞, where the norm is the largest |x_i|. All of these are convex.

The formula only gives a true norm, and a convex ball, when p ≥ 1. Try p = ½, where the "ball" is √|x₁| + √|x₂| ≤ 1. The points (1, 0) and (0, 1) are in it. Their midpoint (0.5, 0.5) scores √0.5 + √0.5 ≈ 1.414 > 1, so it is outside. The shape is a four-pointed star, and you can drag across its dent in the widget.

Polyhedra. A polyhedron is what remains after cutting space with finitely many flat walls: {x : Ax ⪯ b, Cx = d}. The sign ⪯ means "≤ in every entry". Each row of A ⪯ b is one halfspace, and each row of Cx = d is one hyperplane. A bounded polyhedron is a polytope.

For example, x₁ ≥ 0, x₂ ≥ 0 and x₁ + x₂ ≤ 1 cut out a triangle. In matrix form, the rows of A are (−1, 0), (0, −1) and (1, 1), and b = (0, 0, 1).

Simplices. A simplex is the convex hull of affinely independent points (no point is an affine combination of the others): a segment, a triangle, a tetrahedron. The probability simplex holds the vectors with entries ≥ 0 that add up to 1, such as (0.2, 0.3, 0.5): every possible list of probabilities or mixing weights.

5. Three cones you will meet again

The ice-cream cone. The second-order cone is {(x, t) : ‖x‖₂ ≤ t}, where x ∈ ℝⁿ and t is a number. Each horizontal slice at height t is a disk of radius t, so the slices grow as you go up. CO4A 2-2-2

Check three points. (3, 4, 5) is on the surface, since ‖(3, 4)‖₂ = 5. (3, 4, 6) is inside, and (3, 4, 4) is outside.

This cone describes friction. A box on a floor will not slide while the sideways force stays within μ times the pushing-down force, where μ is the friction coefficient. With μ = 0.5 and 10 N pressing down, any sideways force up to 5 N, in any direction, is safe. Legged robots and grasping hands keep their contact forces inside such cones. Day 5 builds second-order cone programs from them.

The PSD cone. A symmetric matrix X is positive semidefinite (PSD), written X ⪰ 0, when all its eigenvalues are ≥ 0. Equivalently, v^TXv ≥ 0 for every vector v. The PSD matrices of size n form the PSD cone Sⁿ₊.

It is a convex cone for a one-line reason. If v^TXv ≥ 0 and v^TYv ≥ 0, then v^T(θ₁X + θ₂Y)v ≥ 0 for any θ₁, θ₂ ≥ 0.

For 2×2 matrices there is a quick test:

[a b; b c] ⪰ 0 ⇔ a ≥ 0, c ≥ 0, ac ≥ b²

[2 1; 1 1] passes, since 2·1 = 2 ≥ 1. Its eigenvalues are 0.382 and 2.618. [2 1.5; 1.5 1] fails, since 2·1 = 2 < 2.25. One of its eigenvalues is −0.081.

Every covariance matrix is PSD, including the uncertainty a car's localisation keeps about its own position. Problems over PSD matrices are semidefinite programs, also on Day 5.

6. Building new convex sets: the safe moves

Checking the segment test by hand gets tiring for complicated sets. Instead, build them from the basic shapes with moves that keep convexity. Boyd's slides list three ways to prove a set convex: the definition, convex functions (Day 3), or safe moves from simple sets. Boyd EE364a slides CO4A 2-3

Move 1: intersection. If a segment stays inside each of two sets, it stays inside their overlap. So the intersection of convex sets is convex, even of infinitely many. The unit disk intersected with the halfspace x₂ ≤ 0.5 is a "cut disk", one of the widget's shapes. In fact every closed convex set is the intersection of all the halfspaces that contain it. Boyd & Vandenberghe 2004

Move 2: affine maps. An affine map f(x) = Ax + b scales, rotates, shears, shifts or projects. It sends segments to segments. So the image of a convex set is convex, and so is the preimage {x : f(x) ∈ D} of a convex set D.

For example, A = diag(2, 1) stretches the unit disk into the ellipse x₁²/4 + x₂² ≤ 1, with semi-axes 2 and 1. A preimage example: the set {x : x₁A₁ + … + x_kA_k ⪯ B}, called a linear matrix inequality, is convex. It is the set of x for which the affine map x ↦ B − Σx_iA_i lands in the PSD cone. CMU 10-725 2016

Move 3: the perspective map. A pin-hole camera divides by depth. The perspective function P(z, t) = z/t, defined for t > 0, does the same. Images and preimages of convex sets under it are convex. For example, P(2, 4, 2) = (2/2, 4/2) = (1, 2). A convex object in front of a camera therefore shows up as a convex blob in the picture.

Following an affine map by the perspective gives a linear-fractional function f(x) = (Ax + b)/(c^Tx + d), defined where c^Tx + d > 0. It also keeps convex sets convex.

Not a safe move: union. Two disks side by side are each convex, but together they are not. In the widget's "two disks" shape, the midpoint of the two centres lies outside both.

A driving example. A planner for a self-driving car has to describe the free space it may drive into. The simplest pieces are halfspaces, one per straight edge: the two road edges, the lane line, and a line some metres behind the car ahead. Intersect a few of them and you get a convex piece of free space.

Here region A is the road with x ≤ 20, which keeps 10 m behind the car ahead. Region B is the road with y ≥ 1.75, the left lane. Each is a polyhedron, so each is convex.

The free space as a whole is "A or B", a union, and it is not convex. The point p = (15, 0) is in A and q = (45, 3.5) is in B. At x = 25 their segment has y ≈ 1.17: still in the right lane, but past the 20 m line. It is in neither region.

Planners handle this by picking a convex piece, or a chain of them, and optimizing inside it. Picking is a yes-or-no decision (Day 15). Day 19 builds whole planners this way, following the IRIS paper and graphs of convex sets.

7. Separating and supporting hyperplanes

Picture two fields that do not overlap. If both are convex, you can always build a straight fence between them. If one has a dent that wraps around the other, a straight fence may be impossible.

That is the separating hyperplane theorem. If C and D are nonempty convex sets with no common point, then there are a ≠ 0 and b such that Boyd & Vandenberghe 2004 CO4A 2-5

a^Tx ≤ b for all x ∈ C, and a^Tx ≥ b for all x ∈ D

For example, take the disks of radius 1 centred at (−2, 0) and (2, 0). The wall x₁ = 0, meaning a = (1, 0) and b = 0, separates them. Every point of the left disk has x₁ ≤ −1, and every point of the right disk has x₁ ≥ 1.

Now lay a ruler against a convex pebble. At any point on its edge, the ruler can touch the pebble with the whole pebble on one side. That is a supporting hyperplane. If C is convex and x₀ is on its boundary, there is an a ≠ 0 with a^Tx ≤ a^Tx₀ for every x ∈ C.

For example, take the unit disk and x₀ = (1, 0). The vertical line x₁ = 1, with a = (1, 0), supports it: every point of the disk has x₁ ≤ 1.

Read the separating theorem in one direction only. The sets C = D = {0} on the number line share a point, yet the "wall" x = 0 satisfies both inequalities. And a wall with a gap on both sides is not guaranteed. Two disjoint closed convex sets can creep arbitrarily close to each other, leaving no room for a gap. Boyd & Vandenberghe 2004

These two theorems are the geometry behind duality (Day 9) and behind linear classifiers, which look for a wall between two classes of data points. CO4A 2-5

Robot planners use the same fact. A car's footprint and an obstacle, both convex polygons, are apart exactly when some straight line separates them with a gap. Day 19 and the OBCA paper turn this into constraints a solver can handle.

1. The chord test: what "convex" means for a function

Draw the graph of a function. Pick two points on it and stretch a string between them. If the string never passes below the graph, wherever you pick the points, the function is convex. A bowl-shaped valley passes. A road over a hill fails, because some string would cut through the hill.

The string is called a chord. To write it down, call the two inputs x and y. Every point between them is θx + (1 − θ)y for some weight θ ∈ [0, 1]: θ = 1 gives x, θ = 0 gives y, and θ = ½ gives the midpoint. The chord's height above that point is the same blend of the two heights, θf(x) + (1 − θ)f(y).

A function f is a convex function if its domain dom f (the inputs where it is defined) is a convex set and the curve never rises above the chord: CO4A 3-1-1 Boyd & Vandenberghe 2004 §3.1.1

f(θx + (1 − θ)y) ≤ θf(x) + (1 − θ)f(y)

This must hold for all x and y in dom f and every θ ∈ [0, 1]. A worked pair: f(x) = x², x = −1, y = 3 and θ = ½. The midpoint is 1, where the curve is at f(1) = 1. The chord there is at (1 + 9)/2 = 5. The gap, chord minus curve, is 4, so this pair passes.

One passing pair proves nothing, because the rule covers every pair. One failing pair, though, settles the question. Take f(x) = x⁴ − 2x² with x = −1 and y = 1. Both ends sit at height −1, so the chord is flat at −1, but the curve at the midpoint is f(0) = 0. The curve is 1 above the chord, so this f is not convex.

Turn the picture upside down for the mirror idea. A function is a concave function if −f is convex: the chord never passes above the graph, like a dome. A straight line f(x) = a^Tx + b (where a^Tx = a₁x₁ + … + aₙxₙ) turns the inequality into an equality, so affine functions are both convex and concave. CO4A 3-1-1

With θ = ½ the chord test reads "f of the average ≤ the average of f". Extended to many points and to expected values, this is Jensen's inequality, f(E[X]) ≤ E[f(X)]. If X is 0 or 2 with equal odds and f(x) = x², then f(E[X]) = f(1) = 1, while E[f(X)] = (0 + 4)/2 = 2. CO4A 3-1-3 Boyd & Vandenberghe 2004 §3.1.8

Day 2 gives a second picture of the same idea. Shade everything on or above the graph. That region is the epigraph, epi f = {(x, t) : f(x) ≤ t}, where t stands for a height. A function is convex exactly when its epigraph is a convex set. The chord test is Day 2's segment test applied to two points on the graph: a chord that dips below the graph leaves the shaded region. CO4A 3-1-3 Boyd & Vandenberghe 2004 §3.1.7

Try it: pick a function and drag the two points (or use the sliders). The θ slider moves the pair of dots that compare chord and curve. Press Hunt for a bad chord to scan 41 × 41 pairs of points with 11 weights each. The coloured strip under the plot shows the sign of f″, which section 3 explains.

Notice that the default, x² with a = −1 and b = 3, shows a gap of 4.000 at x = 1 (chord 5.000 − curve 1.000) and a lowest gap of 0.000, so the chord is green. Switch to x³: its default chord from 0 to 2 is green as well, but Hunt finds a chord from −2 to 1.1 that dips 3.98 below the graph.

2. The tangent test: one slope gives a floor everywhere

Lay a ruler against a bowl so that it touches at one spot. The ruler touches the bowl there and nowhere cuts into it. For a convex function the tangent line is that ruler: it touches the graph at one point and stays on or below the graph everywhere else.

In symbols, write ∇f(x) for the gradient, the vector of slopes ∂f/∂x₁, …, ∂f/∂xₙ, which points uphill. The tangent at x predicts f(y) ≈ f(x) + ∇f(x)^T(y − x). The first-order condition says that a differentiable f on a convex domain is convex exactly when this prediction never overshoots: CO4A 3-1-3 Boyd & Vandenberghe 2004 §3.1.3

f(y) ≥ f(x) + ∇f(x)^T(y − x) for all x, y

A number check. For f(x) = x² at x = 1 the slope is 2, so the tangent line is 1 + 2(y − 1) = 2y − 1. At y = 3 the tangent says 5 and the true value is 9, so the floor holds. Now take sin x at x = 0. Its tangent is the line through the origin with slope 1. At y = 1 the line says 1, but sin 1 ≈ 0.841 is lower. The tangent has poked above the graph, so sin is not convex.

Here is why this matters. Suppose the gradient is zero at some point x*. The floor then reads f(y) ≥ f(x*) for every y, so nothing anywhere beats x*. For a convex function, a flat spot is a global minimum, not just a local one. That is Day 1's "local equals global", seen through tangents. Boyd and Vandenberghe call the first-order condition "perhaps the most important property of convex functions". Boyd & Vandenberghe 2004 §3.1.3

For a planner this is a guarantee you can act on. If a car's cost is convex in its plan and there are no constraints, a zero gradient means the best plan, not just a nearby dip. With constraints, the zero-gradient test becomes the KKT conditions of Day 10. A kink, like |x| at 0, has no single tangent, but it still has many lines that act as floors; Day 7 calls their slopes subgradients.

3. The curvature test: second derivatives and the Hessian

Drive along the graph from left to right. On a convex function the steering wheel only ever turns left (the road bends upward) or stays straight. It never turns right. Bending upward means the slope keeps growing.

The rate at which the slope grows is the second derivative f″. So a function of one variable, twice differentiable on an open interval, is convex exactly when f″(x) ≥ 0 everywhere. Boyd & Vandenberghe 2004 §3.1.4 For −log x on x > 0, f′ = −1/x and f″ = 1/x² > 0, so −log x is convex. For x⁴ − 2x², f″ = 12x² − 4, which is negative for |x| < 0.577. There the graph bends down, and that is exactly where the chord in the first figure failed.

With several variables, curvature becomes a matrix. The Hessian ∇²f(x) is the table of all second derivatives; its entry (i, j) is ∂²f/∂xᵢ∂xⱼ. "Bends up in every direction" means the Hessian is positive semidefinite, written ∇²f(x) ⪰ 0: all its eigenvalues are ≥ 0, or equivalently v^T∇²f(x)v ≥ 0 for every direction v. The second-order condition, for f twice differentiable on a convex open domain, is the n-variable version of f″ ≥ 0: CO4A 3-1-3 Boyd & Vandenberghe 2004 §3.1.4

f convex ⇔ ∇²f(x) ⪰ 0 for every x in dom f

Two quadratics that look alike show the difference. Write a matrix row by row, with ";" between rows. f(x) = x₁² + x₁x₂ + x₂² has the constant Hessian [2 1; 1 2], with eigenvalues 1 and 3. Both are positive, so f is a bowl. g(x) = x₁² + 4x₁x₂ + x₂² also has only positive coefficients, but its Hessian [2 4; 4 2] has eigenvalues −2 and 6. One negative eigenvalue makes g a saddle: up in one direction, down in another.

You can see the saddle without eigenvalues. Walk through the origin along the direction v = (1, −1), so that x = (t, −t). Then g = t² − 4t² + t² = −2t², a parabola that opens downward. A convex function must be convex along every line, so g is not convex. This trick, restricting to a line, turns any n-variable question into a one-variable one. B&V slides 3.7 CMU 10-725 2016

Try it: start from the saddle, then pick the bowl, or slide q₁₂ toward 0 and watch the verdict turn convex once |q₁₂| ≤ 1. With q₁₁ = q₂₂ = 1 the Hessian eigenvalues are 2 − 2q₁₂ and 2 + 2q₁₂.

Notice that the default saddle, x₁² + 4x₁x₂ + x₂², has Hessian eigenvalues −2 and 6. Its contours are hyperbolas, and the slices read f = 6t² along (1, 1) but f = −2t² along (1, −1), the same −2t² as in the text.

How curved is curved enough? Convexity allows straight stretches: the chord of |x| from 1 to 3 lies on the graph, gap 0. A function is strictly convex if the chord lies strictly above the graph between any two different points (for 0 < θ < 1). Then it cannot have two lowest points, because the midpoint between them would be lower still. A quick test: if ∇²f(x) has all eigenvalues > 0 everywhere, f is strictly convex. The converse fails: x⁴ is strictly convex, yet f″(0) = 0. CO4A 3-1-1 Boyd & Vandenberghe 2004 §3.1.4

A function is strongly convex with parameter m > 0 if it is at least as curved as the bowl (m/2)‖x‖₂², where ‖x‖₂ = √(x₁² + … + xₙ²) is ordinary length. Formally, f(x) − (m/2)‖x‖₂² must still be convex. CO4A 3-1-1 The Hessian of that difference is ∇²f(x) − mI, where I is the identity matrix. So for twice-differentiable f, strong convexity means ∇²f(x) ⪰ mI: every eigenvalue of the Hessian is at least m, at every x.

Some numbers. x² is strongly convex with m = 2: x² − x² = 0 is still convex, while m = 3 would leave −½x². The bowl x₁² + x₁x₂ + x₂² has smallest Hessian eigenvalue 1, so m = 1. x⁴ is strictly but not strongly convex, because its curvature 12x² falls to 0 at x = 0. So strongly convex ⇒ strictly convex ⇒ convex, and neither arrow reverses. CMU 10-725 2016 Why care? Strong convexity puts a curved floor under f instead of a flat one, f(y) ≥ f(x) + ∇f(x)^T(y − x) + (m/2)‖y − x‖₂². Day 6 uses this floor to show that gradient descent then shrinks its error bound by a constant factor at every step.

4. A catalogue of convex functions

Most convexity proofs end with "and this piece is a function we already know". So it pays to recognise a dozen on sight, the way you read road signs without reading the words. The table lists the standard ones, each with a one-line reason. CO4A 3-1-2 Boyd & Vandenberghe 2004 §3.1.5 CMU 10-725 2016

A norm ‖x‖ measures the length of a vector. Three are common: ‖x‖₁ = |x₁| + … + |xₙ| (city-block length), ‖x‖₂ (straight-line length) and ‖x‖_∞ = maxᵢ|xᵢ| (the largest entry). For x = (3, −4) they give 7, 5 and 4. Every norm is convex for a one-line reason: scaling (‖θx‖ = θ‖x‖ for θ ≥ 0) and the triangle inequality (‖u + v‖ ≤ ‖u‖ + ‖v‖) together give ‖θx + (1 − θ)y‖ ≤ θ‖x‖ + (1 − θ)‖y‖, which is the chord test. In a planner, the distance ‖p − p_goal‖₂ from a car's position p to its goal is therefore convex in p.

A quadratic f(x) = ½x^TPx + q^Tx + r, with P a symmetric matrix, has the same Hessian P at every point. So it is convex exactly when P ⪰ 0. Least squares, ‖Ax − b‖₂², is the case P = 2A^TA. It is always convex, because v^TA^TAv = ‖Av‖₂² ≥ 0 for every v.

Log-sum-exp, lse(x) = log(e^x₁ + … + e^xₙ), is a smooth version of the maximum: max xᵢ ≤ lse(x) ≤ max xᵢ + log n. At x = (1, 2, 3) it is 3.408, between 3 and 3 + log 3 ≈ 4.099. Its Hessian is diag(z) − zz^T, where zᵢ = e^xᵢ/Σⱼe^xⱼ are the softmax weights and diag(z) is the diagonal matrix holding them. At x = (0, 0) that is [0.25 −0.25; −0.25 0.25], with eigenvalues 0 and 0.5, so it is positive semidefinite. Boyd & Vandenberghe 2004 §3.1.5 CMU 10-725 2016

Finally, the indicator function of a set C, I_C(x) = 0 for x in C and +∞ outside, turns a constraint "x ∈ C" into a cost. It is convex exactly when C is a convex set, which ties today back to Day 2. CO4A 3-1-2

5. Building convex functions from safe pieces

Checking a Hessian by hand gets painful fast. Engineers instead build costs like Lego, from catalogue pieces joined by moves known to keep convexity. Five moves cover almost everything. CO4A 3-2 Boyd & Vandenberghe 2004 §3.2

1. Nonnegative weighted sums. If f₁, …, fₘ are convex and every weight αᵢ ≥ 0, then α₁f₁ + … + αₘfₘ is convex: add up their chord inequalities. Example: 3x² + 2|x| + e^x. A negative weight can break it: x² − 2x⁴ has f″ = 2 − 24x², which is negative for |x| > 0.29.
2. Affine composition. If f is convex, so is f(Ax + b) for any matrix A and vector b, because a straight line in x maps to a straight line in Ax + b. Example: ‖Ax − b‖ is convex for every norm.
3. Pointwise maximum. The pointwise maximum f(x) = max{f₁(x), …, fₘ(x)} of convex functions is convex, and so is a supremum over infinitely many. Picture the upper edge of several bowls: still one bowl. Example: f(x) = max(2x + 1, −x + 4) has a kink at x = 1 with f(1) = 3, and the chord from x = 0 (f = 4) to x = 2 (f = 5) passes above it at 4.5. The same rule makes a worst-case cost convex: the largest of several convex costs, one for each guess about what the car ahead will do.
4. Partial minimization. If f(x, y) is convex in x and y jointly and C is a convex set, then g(x) = inf_{y ∈ C} f(x, y) is convex, as long as g stays above −∞. This is partial minimization. Example: the distance to a convex set, dist(x, C) = min_{y ∈ C} ‖x − y‖. For C = [0, 1] it equals max(0, −x, x − 1), so dist(3, [0, 1]) = 2. How far a car's position lies outside a convex lane corridor is therefore a convex penalty. Boyd & Vandenberghe 2004 §3.2.5
5. Composition. For f(x) = h(g(x)), convex pieces are not enough; the direction in which the outer function h slopes matters too. The composition rule comes from the chain rule f″ = h″(g)·g′² + h′(g)·g″: when both terms are ≥ 0, f is convex. The table lists the four cases. CMU 10-725 2016 Boyd & Vandenberghe 2004 §3.2.4

A failure shows why the slope direction matters. h(u) = e^−u is convex but decreasing, and g(x) = x² is convex. Their composition, the bell curve e^−x², is not convex: at x = 0 it is 1, above the chord from −1 to 1, which sits at e⁻¹ ≈ 0.368.

One move looks safe and is not: a plain minimum. min(x², (x − 2)²) is 0 at x = 0 and x = 2 but 1 at x = 1, above the chord at 0. Two bowls side by side make two valleys. Partial minimization is different: one jointly convex function, minimized over a convex set. CO4A 3-2 opens its minimization part by calling the minimum of convex functions convex; that sentence is wrong, though the formula after it is right. CO4A 3-2

Here is the toolkit on a real cost. A car's speed planner chooses accelerations a₁, …, a_N, one every Δt seconds. A comfort cost penalizes acceleration and jerk, the change in acceleration per second, (a_k+1 − a_k)/Δt:

J(a) = Σ a_k² + w Σ ((a_k+1 − a_k)/Δt)²

Each term is a square, which is convex, of an affine function of a (move 2), and the weights 1 and w are nonnegative (move 1). So J is convex for every w ≥ 0 and Δt > 0, with no Hessian in sight. Numbers, with N = 3, Δt = 1 and w = 1: plan A = (0, 2, 0) costs 4 + 8 = 12, plan B = (2, 0, 2) costs 8 + 8 = 16, and their average (1, 1, 1) costs 3 + 0 = 3, far below the average cost of 14. Blending two plans never costs more than blending their costs. That is the chord test again, and it is why a model predictive controller (Day 17) can find the single best plan reliably.

6. Two looser cousins: quasiconvex and log-convex functions

Some useful functions are not bowls but still have only one valley. Count the 10-tonne trucks needed to move x tonnes: 0 for nothing, 1 for up to 10 tonnes, 2 for up to 20. The graph is a staircase, ⌈x/10⌉. It fails the chord test, since 5 tonnes needs 1 truck, above the chord's 0.5 between 0 and 10 tonnes. Yet it never goes down, so it never makes a second valley.

The idea that captures this is the sublevel set {x : f(x) ≤ α}, all the inputs that score at most α. A function is quasiconvex if every sublevel set is convex. In one variable that means the graph goes down and then up, or only one way, with no second dip. There is a chord-style test too, with the higher endpoint in place of the chord: CO4A 3-4-1 CO4A 3-4-2

f(θx + (1 − θ)y) ≤ max{f(x), f(y)}

Every convex function is quasiconvex, because its sublevel sets are convex. CO4A 3-1-3 The reverse fails. f(x) = √|x| has the sublevel sets [−α², α²], all intervals, so it is quasiconvex. It is not convex: between 0 and 4 the chord at x = 2 is at 1, while f(2) ≈ 1.414. But 1.414 ≤ max{f(0), f(4)} = 2, so the weaker test passes. Boyd & Vandenberghe 2004 §3.4.1

Weaker means fewer guarantees. A zero slope no longer marks a minimum: x³ only rises, so it is quasiconvex, and its slope at 0 is 0, but 0 is not a minimizer. For differentiable f, the correct first-order test is: if f(y) ≤ f(x), then ∇f(x)^T(y − x) ≤ 0. Every point at least as good as x lies on the downhill side of the gradient. CO4A 3-4-3 garbles this statement, so follow Boyd and Vandenberghe here. Boyd & Vandenberghe 2004 §3.4.3

Two more warnings. A local minimum need not be global: at x = 5 tonnes nothing nearby needs fewer trucks, yet x = 0 needs none. Boyd & Vandenberghe 2004 §4.2.2 And sums can break quasiconvexity: x³ and −x each run one way, but x³ − x has the sublevel set {x : x³ − x ≤ 0} = (−∞, −1] ∪ [0, 1], two pieces.

The second cousin works on a log scale. A positive function f is log-concave if log f is concave, and log-convex if log f is convex. The bell curve e^−x²/2 is the standard log-concave example: it is not concave everywhere (its tails bend up), but its log, −x²/2, is a dome. The Gaussian cumulative distribution function and det X on positive definite matrices are log-concave too. On the other side, 1/x on x > 0 is log-convex, since its log, −log x, is convex. CO4A 3-5 Boyd & Vandenberghe 2004 §3.5.1

Four facts are worth keeping. A log-convex f is convex, since f = e^{log f} is a convex, nondecreasing function of the convex log f. Products of log-concave functions are log-concave, since their logs add. Sums of log-convex functions stay log-convex: log(f + g) is log-sum-exp applied to (log f, log g), which is convex and nondecreasing in each argument. (CO4A 3-5 argues this step with "the log of a convex function is convex", which is false.) Sums of log-concave functions need not stay log-concave: two bells at ±2, e^{−(x−2)²/2} + e^−(x+2)²/2, have (log f)″(0) ≈ +3 > 0. Boyd & Vandenberghe 2004 §3.5.2

7. How to check convexity in practice

Faced with a new cost, work down this list and stop at the first step that answers. It is the recipe from the Boyd and Vandenberghe slides, cheapest step first. B&V slides 3.18

1. Recognise it. Is it in the catalogue of section 4? Then it is convex.
2. Build it. Can you write it as catalogue pieces joined by the five moves of section 5? Then it is convex, and you never needed a derivative.
3. Differentiate it. For a function of one or two variables, compute f″ or the Hessian and check its sign or eigenvalues at every point of a convex domain.
4. Slice it. Otherwise restrict to lines: g(t) = f(x + tv) must be convex in t for every point x and direction v.
5. Attack it. If you suspect it is not convex, hunt for one bad chord, as the widget does, or one direction with a negative Hessian eigenvalue. One failure settles it. No number of passing samples proves convexity.

Steps 1 and 2 are automated in CVXPY, a Python modelling library, under the name disciplined convex programming (DCP). Every built-in function (an "atom", such as square, norm or log_sum_exp) carries its known curvature and monotonicity, and CVXPY applies the composition rule at every node of your expression (as of Oct 2026). CVXPY docs If the rules certify an expression, it is convex. If they cannot, CVXPY reports UNKNOWN, which may only mean the expression was written in an unlucky way: DCP is sufficient, not necessary. B&V slides 3.33

import cvxpy as cp
x = cp.Variable()
a = cp.Variable(3)                        # a 3-step acceleration plan

comfort = cp.sum_squares(a) + 4 * cp.sum_squares(cp.diff(a))  # w/Δt² = 4
print(comfort.curvature)                  # CONVEX
print(cp.log_sum_exp(a).curvature)        # CONVEX
print(cp.sqrt(1 + cp.square(x)).curvature)       # UNKNOWN, though convex
print(cp.norm(cp.hstack([1, x]), 2).curvature)   # CONVEX: same function
print(cp.minimum(cp.square(x), cp.square(x - 2)).curvature)  # UNKNOWN: truly not convex

The √(1 + x²) line is the textbook case. The function equals ‖(1, x)‖₂, a norm of an affine expression, which CVXPY certifies at once when it is written that way. Day 4 takes the next step and builds whole optimization problems from convex pieces.

1. Writing a problem down: standard form

Every optimization problem says the same thing. Pick some numbers to make one number as small as possible, while obeying some rules. A delivery robot picks a route to make the travel time small without leaving the sidewalk.

Standard form is the agreed way to write this down, so that theorems and solvers can read it. CO4A 4-1

minimize f(x)
subject to g_i(x) ≤ 0, i = 1, …, m
h_j(x) = 0, j = 1, …, r

Here x ∈ ℝⁿ is the optimization variable, the numbers you get to choose. f is the objective (CO4A also calls it the criterion), the number you want small. Each g_i(x) ≤ 0 is an inequality constraint and each h_j(x) = 0 is an equality constraint. All the rules are written as "something ≤ 0" or "something = 0", so a rule like x ≤ 3 becomes x − 3 ≤ 0.

A point that obeys every rule is a feasible point, and all of them together form the feasible set. The smallest objective value over the feasible set is the optimal value f*. A feasible point x* with f(x*) = f* is an optimal point, also called a solution or minimizer.

Two more names help. A feasible x with f(x) ≤ f* + ε is ε-suboptimal; solvers return such points, because they stop at a tolerance. An inequality with g_i(x) = 0 is active at x: you are touching that fence. CO4A 4-1

The problem is a convex problem when f is convex, every g_i is convex, and every h_j is affine, so that all the equalities together read Ax = b. Why must equalities be affine? The rule h(x) = 0 is the same as two rules, h(x) ≤ 0 and −h(x) ≤ 0. Both are convex only when h is both convex and concave, which means affine. Tibshirani 2016 To maximize a concave function, minimize its negative instead.

Worked example: maximize 3x + 2y subject to x + y ≤ 4, x ≤ 3, x ≥ 0 and y ≥ 0. In standard form it reads: minimize −3x − 2y subject to x + y − 4 ≤ 0, x − 3 ≤ 0, −x ≤ 0 and −y ≤ 0. Linear functions are both convex and affine, so this is a convex problem. The feasible set is a four-sided region with corners (0, 0), (3, 0), (3, 1) and (0, 4), where 3x + 2y is 0, 9, 11 and 8. The best point is (3, 1), with maximum 11, so f* = −11 in the minimize form. There the rules x + y ≤ 4 and x ≤ 3 are active, and x ≥ 0, y ≥ 0 are not. (Day 5 explains why a corner always wins here.)

A driving example in one dimension: a lane-keeping planner chooses a lateral offset d, in metres from the lane centre. It would like d = 0.8 to give a cyclist room, but the lane edge allows at most d = 0.5. In standard form: minimize (d − 0.8)² subject to d − 0.5 ≤ 0. The answer is d* = 0.5 with the rule active, and f* = 0.09. Keep this example in mind for section 3.

2. Why local means global

Imagine walking downhill in thick fog. You can only feel the slope under your feet, so you stop when every step goes up. In a single round valley you are now at the bottom. In the Alps you might be standing in a small hollow halfway up a mountain.

A feasible x is a local optimum if no feasible point within some distance R > 0 is better: f(x) ≤ f(y) for all feasible y with ‖y − x‖₂ ≤ R. It is a global optimum if no feasible point anywhere is better. For a convex problem, every local optimum is a global optimum. Tibshirani 2016 Boyd & Vandenberghe 2004

The one-line reason: suppose some feasible y were better than a local optimum x. The straight segment from x to y stays feasible, because the feasible set is convex. Along it, f stays below the chord, because f is convex. So points on the segment just next to x are also better than x, and x was not a local optimum after all.

Worked example. The function f(x) = x⁴ − 3x² + x is not convex. Its slope is zero at three points: x ≈ −1.301 (f ≈ −3.514, the global minimum), x ≈ 0.170 (a local maximum) and x ≈ 1.131 (f ≈ −1.070, a local minimum that is not global). Gradient descent from x = 2 with step 0.01 slides into the hollow at 1.131 and stops there. The convex function x² + x has exactly one point of zero slope, x = −0.5, and that point is the global minimum.

Two more facts. First, the set of all optimal points is convex: if x and y are both optimal, every point between them is too. For example, f(x) = max(|x| − 1, 0) is flat on [−1, 1], and that whole interval is optimal. Second, if f is strictly convex (no flat stretches), there is at most one optimal point. CO4A 4-2

3. The first-order optimality condition

Back to the lane-keeping planner. At d* = 0.5 the slope of the cost is f′(0.5) = 2(0.5 − 0.8) = −0.6, not zero. The cost would keep falling if d could grow, but the lane edge stops it. Every allowed move goes to some y ≤ 0.5, and the predicted change f′(d*)(y − d*) = −0.6(y − 0.5) is never negative. So every allowed move goes uphill, and d* is optimal.

The general rule is the first-order optimality condition. Take a convex problem with a differentiable f and feasible set C. Then x is optimal exactly when x ∈ C and CO4A 4-3

∇f(x)^T(y − x) ≥ 0 for every y ∈ C

Read ∇f(x)^T(y − x) as "the first-order change in f when you start moving from x toward y". The condition says no feasible direction goes downhill. Picture it: −∇f(x) is the steepest-downhill arrow, and at the optimum it points out of C. The line through x at right angles to ∇f(x) is a supporting hyperplane of C: it touches C at x, and all of C lies on its uphill side. Boyd & Vandenberghe 2004

With no constraints at all (C = ℝⁿ), y − x can point anywhere, including along −∇f(x). Then the condition becomes the familiar ∇f(x) = 0. Example: f(x) = ½x^TQx + b^Tx with Q = diag(2, 4) and b = (−2, −8). Setting ∇f(x) = Qx + b to zero gives x* = (1, 2). Tibshirani 2016

Now a two-dimensional example with constraints. Minimize f(x) = (x₁ − 3)² + (x₂ − 2)² over the triangle x₁ ≥ 0, x₂ ≥ 0, x₁ + x₂ ≤ 2. The level sets of f are circles around the target (3, 2), which lies outside the triangle. The best point is the triangle's closest point to the target, x* = (1.5, 0.5), with f* = 1.5² + 1.5² = 4.5.

Check the condition. The gradient is ∇f(x*) = (2(1.5 − 3), 2(0.5 − 2)) = (−3, −3). For any y in the triangle:

∇f(x*)^T(y − x*) = −3(y₁ + y₂ − 2) ≥ 0

It is non-negative because y₁ + y₂ ≤ 2. At the corner y = (0, 0) it equals 6, and anywhere on the edge x₁ + x₂ = 2 it equals 0. Compare a point that is not optimal, x = (1, 1), where f = 5 and ∇f = (−4, −2). Toward y = (2, 0) the product is −4·1 − 2·(−1) = −2 < 0. That direction goes downhill: the point (1.2, 0.8) has f = 4.68 < 5.

Try it: drag the orange target and watch x* and the arrow −∇f(x*). The pink half-plane holds the points y with ∇f(x)^T(y − x) < 0, which lie downhill to first order. Then pick "a candidate point x" and drag x around; the widget shows the corner that goes most downhill.

Notice that with the target at (3, 2), the best point is x* = (1.5, 0.5), with f* = 4.5 and ∇f(x*) = (−3, −3), which is not zero. The arrow −∇f = (3, 3) points straight out through the edge x₁ + x₂ = 2. The smallest value of ∇f(x*)^T(y − x*) over the triangle is 0, so the pink downhill side does not enter the triangle.

If the only constraints are equalities Ax = b, the condition becomes ∇f(x) + A^Tu = 0 for some vector u. This is the Lagrange condition from calculus, and its numbers u return as dual variables on Day 9 and in the KKT conditions on Day 10. Tibshirani 2016

Here is the triangle problem in CVXPY, a Python modelling tool you will see all course:

import cvxpy as cp
x = cp.Variable(2)
f = cp.sum_squares(x - [3, 2])          # bowl centred on the target (3, 2)
rules = [x >= 0, cp.sum(x) <= 2]        # the triangle
prob = cp.Problem(cp.Minimize(f), rules)
prob.solve()
print(x.value, prob.value)              # about [1.5 0.5] and 4.5
print(rules[1].dual_value)              # about 3: how hard the edge pushes back

That last number, 3, is the fence's push. At x*, ∇f + 3·(1, 1) = (−3, −3) + (3, 3) = 0, so the slope and the push cancel.

4. Equivalent problems and slack variables

Equivalent problems are two problems where solving one quickly gives you a solution of the other. They can look very different, with different sizes or different smoothness. Boyd & Vandenberghe 2004 The rest of today is a toolbox of such rewrites, each used to make convexity visible or to fit what a solver accepts.

The simplest rewrite applies an increasing function to the objective. Minimizing ‖Ax − b‖₂ and minimizing ‖Ax − b‖₂² give the same x, because squaring a non-negative number keeps its order. Only the squared version is differentiable everywhere, which is why least squares uses it. CO4A 4-5

A slack variable turns an inequality into an equality by naming the leftover room. Think of an airline's 23 kg suitcase limit: the slack is how many kilograms you still have spare. In symbols: CO4A 4-7

g_i(x) ≤ 0 ⇔ g_i(x) + s_i = 0, s_i ≥ 0

Worked example: x + y ≤ 4 becomes x + y + s = 4 with s ≥ 0. At (3, 1), s = 0, so the rule is active. At (1, 1), s = 2: there are 2 units of room. In driving terms, a following-distance rule "gap ≥ 20 m" becomes 20 − gap + s = 0, and the slack s is the number of spare metres to the car ahead. Slacks are how linear programs reach their standard form on Day 5.

One warning. The new equality g_i(x) + s_i = 0 is affine only if g_i is affine. Add a slack to a curved rule such as x² ≤ 1 and you get x² + s = 1, a curved equality, and the problem stops being convex in standard form. CO4A 4-7

5. Epigraph form: move the objective into the rules

Think of the graph of f as a tent roof. Its epigraph is everything on or above the roof, as on Day 3. Finding the lowest point of the roof is the same as finding the lowest height t at which a point (x, t) still fits in the epigraph.

That gives the epigraph form of a problem. Add one new variable t and minimize it: Boyd & Vandenberghe 2004 Boyd slides

minimize t over (x, t)
subject to f(x) − t ≤ 0, g_i(x) ≤ 0, Ax = b

The objective is now linear, and all the curvature has moved into the constraints. If f is convex, f(x) − t ≤ 0 is a convex rule, so the new problem is still convex.

Worked example: minimize max(2x + 1, 4 − x). The objective has a kink where the two lines cross, so it is not differentiable there. In epigraph form it becomes: minimize t subject to 2x + 1 ≤ t and 4 − x ≤ t. That is a linear program. The two lines meet where 2x + 1 = 4 − x, so x = 1 and t = 3. In planning, the same trick turns "make the worst of several errors small" into a linear program.

6. Substitutions: change of variables and removing equalities

Some problems are convex in disguise. A change of variables x = φ(z), with φ one-to-one and covering the feasible set, gives an equivalent problem in z that may be visibly convex. CO4A 4-5

The classic case is the geometric program (GP). A monomial is γx₁^a₁⋯x_n^a_n with γ > 0 and positive x. A posynomial is a sum of monomials, each with its own γ and exponents. A GP minimizes a posynomial subject to posynomials ≤ 1 and monomials = 1. With z_i = log x_i, each monomial becomes e^{aTz + log γ}, the exponential of an affine function. The log of a posynomial becomes a log-sum-exp, which is convex. Tibshirani 2016

Worked example: minimize x + y subject to xy = 1, with x, y > 0. As written it is not a convex problem, because xy = 1 is a curved equality. Put u = log x and v = log y. The problem becomes: minimize e^u + e^v subject to u + v = 0. Now the objective is convex and the rule is affine. Setting v = −u leaves e^u + e^−u, smallest at u = 0. So x = y = 1, with value 2.

A second substitution removes equality constraints. Find one point x₀ with Ax₀ = b and a matrix M whose columns span the null space of A. Then every solution of Ax = b is x = Mz + x₀, and you can minimize over a free z instead. CO4A 4-6

Worked example: minimize x₁² + x₂² subject to x₁ + x₂ = 2. Take x₀ = (2, 0) and M = (1, −1)^T, so (x₁, x₂) = (2 + z, −z). Minimize (2 + z)² + z²: the derivative 2(2 + z) + 2z is zero at z = −1. So x* = (1, 1) and f* = 2. In large problems this trick is often a poor idea in practice, because M can be costly to compute. Tibshirani 2016 CO4A adds that the new problem can lose sparsity.

7. Optimize in stages, or loosen the rules

Partial optimization means minimizing over some variables first, then over the rest. If f(x, y) is jointly convex and C is convex, then g(x) = min over y ∈ C of f(x, y) is convex in x. So the outer problem is still convex. CO4A 4-4

Worked example: minimize (x − 2)² + (y − x)² with y ∈ [0, 1]. For a fixed x, the best y is the point of [0, 1] closest to x. For x ≥ 1 that is y = 1, leaving (x − 2)² + (x − 1)², smallest at x = 1.5. For x ≤ 1 the cost is at least 1, which is worse. So the solution is (x, y) = (1.5, 1), with value 0.5.

CO4A's main example is the support vector machine. Its slack variables can be minimized out by hand, which leaves the hinge function a₊ = max{0, a} in the objective. A training point whose margin is 0.6, short of the target 1, is charged 1 − 0.6 = 0.4. CO4A 4-4

A relaxation loosens the rules: you solve over a bigger set C̃ ⊇ C. With more choices you can only do as well or better, so the relaxed optimal value is a lower bound on the true one. CO4A 4-8 Example: minimize (x − 0.3)² with x ∈ {0, 1}, which is not a convex set. The answer is x = 0 with value 0.09. Relax to x ∈ [0, 1] and you get x = 0.3 with value 0, a lower bound but not the same answer. Branch and bound on Day 15 (branch-and-bound card) is built on exactly such bounds. CO4A also relaxes a convex but non-affine equality h(x) = 0 to h(x) ≤ 0.

1. Why sort problems into families?

When you ship a parcel, you pick the smallest container that fits: an envelope, a box, a padded box, or a climate-controlled crate. Everything fits in the crate, but the crate costs more and fewer carriers take it.

Convex problems are sorted the same way, by the kind of objective and constraints they use. CO4A orders the families like this, where ⊆ means "is a special case of": CO4A 5

LP ⊆ QP ⊆ QCQP ⊆ SOCP ⊆ SDP ⊆ conic

Why does the family matter? Modelling tools such as CVXPY rewrite your problem into a standard form, "e.g., LP, QP, SOCP, or SDP", and hand it to a solver built for that form. Boyd slides Fewer solvers handle each step up. In CVXPY's solver table, 31 solvers accept LPs, 25 accept QPs, 17 accept SOCPs and 9 accept SDPs (as of Oct 2026). CVXPY 2026

The rule of thumb follows: write your problem in the lowest family that holds it. That keeps the most solvers open to you, including the specialised ones built for simpler forms.

2. Linear programs: flat cost, flat walls

Tilt a flat tray with a marble on it. The marble rolls until it hits a wall, then slides along the wall into a corner. A linear program behaves the same way.

A linear program (LP) has an affine objective and affine constraints: CO4A 5-1

minimize c^Tx
subject to Gx ⪯ h, Ax = b

Here c is the cost vector, and Gx ⪯ h means each entry of Gx is at most the matching entry of h. The feasible set is a polyhedron, a region cut out by finitely many flat walls. Its corners are vertices.

The key fact: if an LP has a solution and its polyhedron has at least one corner, then some vertex is optimal. The one-line reason: the level sets c^Tx = k are parallel lines (or planes). Slide one as far downhill as the region allows. The last part of the region it touches is a vertex, or a whole edge in a tie.

Worked example, the classic diet problem. Tibshirani 2016 Food 1 costs $2 per unit and gives 3 g of protein and 1 mg of a vitamin. Food 2 costs $3 per unit and gives 1 g of protein and 3 mg of the vitamin. You need at least 6 g of protein and 6 mg of the vitamin.

minimize 2x₁ + 3x₂
s.t. 3x₁ + x₂ ≥ 6, x₁ + 3x₂ ≥ 6, x ⪰ 0

The region has three corners. At (0, 6) the cost is $18, at (6, 0) it is $12, and at (1.5, 1.5) it is $7.50. The optimum is 1.5 units of each food, where both nutrient rules are active.

Solvers often want one fixed layout, the standard form LP: minimize c^Tx subject to Ax = b and x ⪰ 0. You reach it with two tricks. Add a slack (Day 4) to each inequality, so Gx ⪯ h becomes Gx + s = h with s ⪰ 0. Split each free variable into two non-negative parts, x = x⁺ − x⁻. Boyd & Vandenberghe 2004 The simplex and interior-point methods (Day 12) work on this form.

Try it: the widget's polygon is x ≥ 0, y ≥ 0, x ≤ 3, x + y ≤ 4 and x + 2y ≤ 6. It maximizes c^Tx for the direction c = (cos θ, sin θ), which is the same as minimizing (−c)^Tx. Slide θ and watch the best corner. Then switch to the QP mode and drag the target.

Notice that at θ = 30° the best corner is (3, 1), where c^Tx = 3.1. Slide θ past 45° and the optimum jumps to (2, 2); past 63.43° it jumps to (0, 3). Exactly at 45°, c is at right angles to the edge x + y = 4, and every point of that edge ties. The optimum jumps whenever c lines up with an edge's outward direction.

A hand check: maximize 3x + 2y over the same polygon. Its corners give 0, 9, 11, 10 and 6, so (3, 1) wins with 11. The direction (3, 2) has angle 33.7°, inside the 0° to 45° range where the widget picks (3, 1).

3. Quadratic programs: a bowl on the tray

Keep the flat walls but make the cost curved. A quadratic program (QP) is CO4A 5-2

minimize ½x^TPx + q^Tx
subject to Gx ⪯ h, Ax = b

The matrix P must be symmetric and positive semidefinite, written P ⪰ 0: u^TPu ≥ 0 for every vector u, or equally, every eigenvalue of P is at least 0. That makes the cost a convex bowl. If P has a negative eigenvalue, the cost bends down in some direction and the problem is no longer convex. Tibshirani 2016 Setting P = 0 gives back an LP, so LP ⊆ QP.

The bowl changes where the optimum can sit. In the widget's QP mode the cost is ½‖x − p‖₂², a round bowl centred on a target p. With p = (1, 1) the bowl's bottom is inside the polygon, so x* = (1, 1) and f* = 0. With p = (3.5, 2), the wall x + y ≤ 4 stops it at (2.75, 1.25), the middle of an edge, with f* = 0.5625. With p = (3, 3) it lands on the corner (2, 2), with f* = 1. So a QP optimum can sit inside, on an edge, or at a corner.

A driving example: smooth speed planning. A car is moving at 20 m/s, and a school zone 4 seconds ahead has a 10 m/s limit. The planner picks speeds v₁, …, v₄ at one-second steps. Harsh braking is uncomfortable, so it penalizes each squared speed change, which is the squared acceleration:

minimize (v₁ − 20)² + (v₂ − v₁)²
+ (v₃ − v₂)² + (v₄ − v₃)²
subject to v₄ ≤ 10

This is a QP: a sum of squares as the cost, one linear rule. The four speed changes must add up to at least −10 m/s, and squares punish big numbers, so the best plan splits the drop evenly. The answer is v = (17.5, 15, 12.5, 10), a steady 2.5 m/s² deceleration, with cost 4 × 2.5² = 25. Braking all at once, v = (20, 20, 20, 10), costs 10² = 100, four times as much. Real planners add speed tracking, jerk and many more steps, but the shape stays a QP. This is the heart of model predictive control on Day 17.

import cvxpy as cp
v = cp.Variable(4)                        # speeds at t = 1, 2, 3, 4 s
cost = cp.square(v[0] - 20) + cp.sum_squares(cp.diff(v))
prob = cp.Problem(cp.Minimize(cost), [v[3] <= 10, v >= 0])
prob.solve(solver=cp.OSQP)                # a QP solver is enough
print(v.value)                            # about [17.5 15. 12.5 10.]

Many familiar problems are QPs: least squares with linear constraints, the lasso, support vector machines, and Markowitz portfolios. Tibshirani 2016 OSQP, a solver built for exactly this form, is the subject of the OSQP paper card.

4. QCQP: curved walls

Now let the walls curve too. A robot vacuum must stay within 1 m of its dock; that rule is a disc, not a polygon. A QCQP (quadratically constrained quadratic program) allows convex quadratic rules: CO4A 5-3

minimize ½x^TP₀x + q₀^Tx + r₀
s.t. ½x^TP_ix + q_i^Tx + r_i ≤ 0, Ax = b

Every P_i must be positive semidefinite, so each rule describes an ellipse-shaped region (an ellipsoid in higher dimensions). Setting P_i = 0 in the rules gives back a QP, so QP ⊆ QCQP.

Worked example: minimize ‖x − (2, 2)‖₂² subject to ‖x‖₂² ≤ 1. The target (2, 2) lies outside the unit disc, so the answer pulls it straight back to the circle: x* = (2, 2)/‖(2, 2)‖₂ = (0.7071, 0.7071). The cost there is 2 × (2 − 0.7071)² = 3.343.

In practice a QCQP rarely gets its own solver. CVXPY's target forms are LP, QP, SOCP and SDP, so a convex QCQP is rewritten as an SOCP, the next family up. Boyd slides

5. Second-order cone programs: length limits

Picture an ice-cream cone standing on its tip. The higher you go, the wider it gets. That shape is the second-order cone, the set of points (x, t) with ‖x‖₂ ≤ t. Some books call it the ice-cream cone.

A second-order cone program (SOCP) has a linear cost and rules of the form "a length must stay below a linear budget": CO4A 5-4 Boyd & Vandenberghe 2004

minimize f^Tx
s.t. ‖A_ix + b_i‖₂ ≤ c_i^Tx + d_i, Fx = g

Each rule says the point (A_ix + b_i, c_i^Tx + d_i) lies inside a second-order cone. With A_i = 0 and b_i = 0, a rule becomes an ordinary linear inequality, so LP ⊆ SOCP.

A driving example: the friction circle. A simple tyre model says the total horizontal acceleration a = (a_lon, a_lat) is limited by grip, ‖a‖₂ ≤ μg. With friction μ = 0.8 and g = 9.81 m/s², the limit is 7.85 m/s². Braking at 6 m/s² while turning at 4 m/s² uses √(36 + 16) ≈ 7.2 m/s², which is allowed. A box rule |a_lon| ≤ 7.85, |a_lat| ≤ 7.85 is linear, but it permits the corner (7.85, 7.85), whose length is 11.1 m/s². That is about 41% more grip than the tyres have.

A second classic SOCP is the robust LP, an LP whose rules must hold for every value of an uncertain coefficient. Boyd & Vandenberghe 2004 Worked example: maximize x₁ + x₂ subject to a^Tx ≤ 2 for every a in the disc {(1, 1) + 0.5u : ‖u‖₂ ≤ 1}, with x ⪰ 0.

The worst a turns (1, 1)^Tx into (1, 1)^Tx + 0.5‖x‖₂, because u^Tx is largest when u points along x. So the rule becomes x₁ + x₂ + 0.5‖x‖₂ ≤ 2, an SOC rule. Without uncertainty the best value of x₁ + x₂ is 2. With it, the best point is x₁ = x₂ = 0.7388, so x₁ + x₂ = 1.4776. Being safe for every a in the disc costs about 26%.

Why is a QCQP also an SOCP? A convex quadratic rule can always be rewritten as a length rule. The simplest case: ‖x‖₂² ≤ 1 is the same set as ‖x‖₂ ≤ 1. The general identity is in the box below. SOCPs also run rockets: powered-descent guidance for landing is solved as an SOCP after a convexification trick, as the lossless convexification paper card explains.

6. Semidefinite programs: matrices that never bend the wrong way

An LP asks every number in a list to be non-negative. A semidefinite program asks a whole symmetric matrix to be positive semidefinite: all its eigenvalues non-negative, so u^TXu ≥ 0 in every direction u. Think of a bowl that curves up, or stays flat, in every direction and never curves down.

A linear matrix inequality (LMI) is a rule x₁F₁ + … + x_nF_n + G ⪯ 0, where the F_i and G are symmetric matrices and "⪯ 0" means negative semidefinite (all eigenvalues ≤ 0). A semidefinite program (SDP) is CO4A 5-5 Boyd & Vandenberghe 2004

minimize c^Tx
s.t. x₁F₁ + … + x_nF_n + G ⪯ 0, Ax = b

The standard form uses a matrix variable X: minimize tr(CX) subject to tr(A_iX) = b_i and X ⪰ 0, where tr is the trace, the sum of the diagonal. If every F_i and G is diagonal, the LMI is just n separate linear inequalities, one per diagonal entry. So LP ⊆ SDP.

Worked example: make the largest eigenvalue of a matrix as small as possible. Write matrices row by row, with rows separated by semicolons.

A(x) = [2 + x, 1; 1, −x]

The rule λ_max(A(x)) ≤ t is the same as A(x) − tI ⪯ 0, an LMI in (x, t). So "minimize t subject to A(x) − tI ⪯ 0" is an SDP. Here the eigenvalues are 1 ± √(1 + (x + 1)²), so the largest is smallest at x = −1, with t* = 2. At x = 0 it would be 1 + √2 ≈ 2.414.

Every SOCP is also an SDP, through a matrix test called the Schur complement: Tibshirani 2016

‖x‖₂ ≤ t ⇔ [tI, x; x^T, t] ⪰ 0

Check it with x = (3, 4), whose length is 5. At t = 5 the block matrix has eigenvalues 0, 5 and 10: positive semidefinite, right on the boundary. At t = 4.9 they are −0.1, 4.9 and 9.9: not positive semidefinite, matching 5 > 4.9.

SDPs reach robotics too. Robot mapping (pose-graph SLAM) has an SDP relaxation that can certify a map is globally optimal; see the SE-Sync paper card.

7. Conic programs, and choosing a solver

All these families fit one template. A cone K is a set that stays inside itself when you scale by any α ≥ 0. The generalized inequality u ⪯_K v means v − u lies in K. A conic program is CO4A 5-6

minimize c^Tx
s.t. Fx + g ⪯_K 0, Ax = b

Choose K = ℝ₊ⁿ, the non-negative numbers, and you get an LP. Choose a product of second-order cones and you get an SOCP. Choose the cone of positive semidefinite matrices and you get an SDP. Solvers such as SCS and Clarabel work in this form. They also accept an exponential cone, which handles log and exp terms like the geometric programs of Day 4. CVXPY 2026

To place a problem, read its pieces from the bottom of the ladder up, and stop at the first family that fits:

1. Objective and every rule affine → LP.
2. Convex quadratic objective (P ⪰ 0), affine rules → QP.
3. Some rules convex quadratic → QCQP, solved as an SOCP.
4. Rules like "‖Ax + b‖₂ ≤ c^Tx + d", or worst cases over a ball → SOCP.
5. Eigenvalues, matrix variables that must be ⪰ 0, matrix norms → SDP.

Here is what that choice means in practice. The table lists a few solvers from CVXPY's solver table, as of Oct 2026. CVXPY 2026

So the speed plan (a QP) can run on OSQP, the friction-circle plan (an SOCP) needs an SOCP solver such as ECOS or Clarabel, and the eigenvalue problem (an SDP) needs an SDP solver such as SCS, Clarabel or MOSEK. Writing a QP as an SDP would still give the right answer. It would just close the door on most of the table.

1. The update: walk downhill

Picture yourself on a hillside in thick fog. You cannot see the valley, but you can feel which way the ground tilts under your feet. A sensible plan is to step in the steepest downhill direction, stop, feel the ground again, and repeat. That loop is the whole algorithm.

On Day 3 you met the gradient ∇f(x), the vector of slopes of f at x. It points in the direction where f rises fastest, so −∇f(x) points where f falls fastest. Gradient descent starts at a point x₀ and repeats one line: CO4A 6-1

x_k+1 = x_k − t∇f(x_k)

Here x_k is the point after k steps. The number t > 0 is the step size; machine learning calls it the learning rate. Steps are long where the ground is steep and short where it is flat, because ∇f itself shrinks near the bottom. You stop when ‖∇f(x_k)‖₂ is tiny.

There is a second way to read the same rule. At x, gradient descent trusts the tangent plane of f, but it charges a fee for moving far from x. The next point minimizes that "tangent plane plus fee" model: Tibshirani 2016

x⁺ = argmin_y f(x) + ∇f(x)^T(y − x) + ‖y − x‖₂²/(2t)

Set the derivative in y to zero and you get x⁺ = x − t∇f(x) again. A small t means a big fee, so the steps are cautious. A big t means you trust the tangent plane far away, which is risky because f curves.

Here is today's running example. Take f(x, y) = ½(x² + 10y²), a bowl ten times steeper in y than in x. Its gradient is ∇f = (x, 10y). From the start (10, 1) the gradient is (10, 10). With t = 0.1 the next point is (10 − 1, 1 − 1) = (9, 0). The steep y-direction is solved in one step. After that each step multiplies x by 0.9: 8.1, 7.29, 6.56, and so on. Boyd & Vandenberghe 2004

2. Fixed step size and the 2/L cliff

Think of a shower knob. Tiny turns take forever to warm the water. Big turns overshoot to scalding, then to freezing, and each swing can be worse than the last. A fixed step, the same t at every iteration, has the same two failure modes.

The safe range depends on how sharply the bowl can curve. A function is L-smooth (its gradient is L-Lipschitz) if the gradient never changes faster than L: CO4A 6-3-1

‖∇f(x) − ∇f(y)‖₂ ≤ L‖x − y‖₂

For a twice-differentiable f this means ∇²f(x) ⪯ LI: no direction curves more than L. So L is the steepest curvature of the bowl. Our example curves by 1 in x and by 10 in y, so L = 10.

Why does trouble start at t = 2/L? On our bowl the two coordinates update separately, each multiplied by its own factor:

x_k+1 = (1 − t)·x_k
y_k+1 = (1 − 10t)·y_k

A coordinate shrinks only if its factor lies strictly between −1 and 1. For y that means 0 < t < 2/10. In general, on a quadratic every curvature λ between the smallest and L gives a factor 1 − tλ, and all of them are safe exactly when 0 < t < 2/L. CO4A 6-2-1

The sign of the factor tells the story, as the figure shows. Between 0 and 1, the coordinate creeps toward 0 from one side. Between −1 and 0, it jumps across the valley floor and back, but each jump is shorter. Below −1, each jump is longer, and the iterates fly away. At t = 2/L exactly the factor is −1, and y flips between +1 and −1 forever.

Here are the real counts from (10, 1), stopping when f ≤ 10⁻⁶. The best fixed step sits in the middle of the safe window, not at its edge.

The best fixed step is t = 2/(m + L), where m is the smallest curvature (here m = 1). It balances the two coordinates so that both shrink by the same factor, (κ − 1)/(κ + 1) = 9/11 per step, where κ = L/m. CO4A 6-3-4 The catch: on a real problem you rarely know m or L. So the next section lets the algorithm find t by itself.

3. Let the algorithm choose t: backtracking and exact line search

Think of reversing into a parking spot. You start with a bold move. If you are about to hit the curb, you ease off and try a shorter one. Backtracking line search does exactly this at every iteration.

Fix two numbers, α and β. At each iteration let g = ∇f(x), start with t = 1, and shrink t ← βt as long as Tibshirani 2016

f(x − tg) > f(x) − α·t·‖g‖₂²

As soon as the inequality is false, you take the step x⁺ = x − tg. The inequality you are waiting for, f(x − tg) ≤ f(x) − αt‖g‖², is called the sufficient-decrease or Armijo condition.

What do α and β mean? The tangent line predicts that a step of length t lowers f by t‖g‖². That is too hopeful, because f curves upward. Boyd and Vandenberghe read α as the fraction of that predicted drop you insist on getting. With α = 0.3 you say "give me at least 30% of what the tangent line promised." β is how hard you cut t after a failure; β = 0.5 halves it. Boyd & Vandenberghe 2004

They take 0 < α < 0.5 and 0 < β < 1, and report typical values of α between 0.01 and 0.3 and β between 0.1 and 0.8. Conventions differ a little: Tibshirani's slides and CO4A also allow α = ½, and this lesson's example uses α = β = ½. CO4A 6-2-2

Here is that first step by hand. At (10, 1) we have f = 55, g = (10, 10) and ‖g‖² = 200. With α = ½ the test asks for f(x − tg) ≤ 55 − 100t.

• t = 1: f = 405, but we need at most −45. Shrink.
• t = 0.5: f = 92.5, but we need at most 5. Shrink.
• t = 0.25: f = 39.4, but we need at most 30. Shrink.
• t = 0.125: f = 38.6 ≤ 42.5. Accept, and move to x₁ = (8.75, −0.25).

The very next step accepts t = 0.5, which is bigger than 2/L = 0.2, and nothing breaks. The test checks each step before it is taken, so an occasional bold step is safe. The full run needs 30 iterations plus 77 rejected trial points. It keeps the O(1/k) guarantee of Section 4, with t replaced by t_min = min{1, β/L}. CO4A 6-3-3

Exact line search goes one step further. Along the downhill ray, it jumps straight to the lowest point: CO4A 6-2-3

t = argmin_{s ≥ 0} f(x − s∇f(x))

For a quadratic with Hessian H this has a formula, t = g^Tg / g^THg. From (10, 1) it gives t = 200/(100 + 1000) = 2/11 at every step. The iterates follow x_k = 10(9/11)^k and y_k = (−9/11)^k, so f drops by (9/11)² ≈ 0.669 per step and needs 45 iterations. This is Boyd and Vandenberghe's own example. Boyd & Vandenberghe 2004

Exact search matched the best fixed step here, but each step required solving a small optimization problem. For a general f that inner problem is expensive, so exact search is rarely worth it. Backtracking gets most of the benefit for a few extra function evaluations. Tibshirani 2016 The whole method fits in a few lines (see the gradient descent card):

def gradient_descent(f, grad, x, alpha=0.5, beta=0.5, tol=1e-6):
    while True:
        g = grad(x)
        if norm(g) <= tol:            # flat enough: stop
            return x
        t = 1.0
        while f(x - t*g) > f(x) - alpha*t*(g @ g):
            t = beta * t               # backtrack: shorter step
        x = x - t*g                    # the gradient step

Try it: start with the default (κ = 10, t = 1/L), then switch the rule to each option and watch the path and the bars. Next, slide κ up to 100 and see every rule slow down. Tap the plot to start somewhere else.

Notice: from (10, 1) at κ = 10, t = 1/L needs 85 iterations, the best fixed step 2/(1+L) and exact line search 45 each, and backtracking (α = β = ½) 30 plus 77 rejected trials. At 2/L, f stalls at 5. At 2.1/L it diverges at iteration 137. On the log chart every converging rule is close to a straight line.

4. How fast? O(1/k), linear rates and the condition number

There are two speed guarantees. A weak one holds for every smooth convex bowl. A strong one holds for bowls that curve upward in every direction.

The weak guarantee rests on a picture. If f is L-smooth, then at any point x a steep bowl of curvature L sits on top of f and touches it at x. That is the quadratic upper bound: CO4A 6-3-2

f(y) ≤ f(x) + ∇f(x)^T(y − x) + (L/2)‖y − x‖₂²

With t ≤ 1/L the step goes downhill on that upper bowl, so f falls by at least (t/2)‖∇f(x)‖² at every step. Adding up these drops gives the theorem for convex, L-smooth f: CO4A 6-3-1

f(x_k) − f* ≤ ‖x₀ − x*‖₂² / (2tk)

Here x* is a minimizer and f* = f(x*). The error falls like a constant over k, written O(1/k) and called a sublinear rate. To guarantee an error ε you need about 1/ε iterations, so one more correct digit costs ten times more steps.

On our bowl, t = 0.1 and ‖x₀ − x*‖² = 101, so the bound is 505/k. Taken at face value it promises 10⁻⁶ only after 5.05 × 10⁸ steps. The actual run needs 85. The bound covers the worst smooth convex function; our bowl is much kinder, because it is also strongly convex.

A function is strongly convex with parameter m > 0 if a shallow bowl of curvature m fits underneath it at every point: f(y) ≥ f(x) + ∇f(x)^T(y − x) + (m/2)‖y − x‖₂². For twice-differentiable f this is ∇²f ⪰ mI. Then every step removes a fixed share of the remaining error: CO4A 6-3-4

f(x_k) − f* ≤ c^k·(L/2)‖x₀ − x*‖₂², 0 < c < 1

This is linear convergence: on a log plot the error is a straight line. Reaching ε needs O(log(1/ε)) iterations, so each extra digit costs the same number of steps. For t = 1/L, CO4A's proof gives c = 1 − m/L. On our bowl c = 0.9, and the bound reaches 10⁻⁶ by iteration 191. The real error is 50·0.81^k, which gets there at 85.

The condition number κ = L/m is the ratio of the steepest curvature to the flattest. It measures how stretched the bowl is. In a round bowl (κ = 1) the downhill arrow points straight at the bottom, and one step of t = 1/L lands on it. In a long narrow canyon the arrow points mostly at the nearest wall, not along the floor. So gradient descent bounces from wall to wall while it creeps along the floor. That bouncing is the zig-zag. Boyd & Vandenberghe 2004

The widget's κ slider shows the cost. From (10, 1), iterations to f ≤ 10⁻⁶:

The count grows roughly in proportion to κ, as Boyd and Vandenberghe's analysis predicts. Boyd & Vandenberghe 2004 Large κ often comes from badly scaled variables. Say a path planner's cost is ½(e² + 100θ²), with lateral error e in metres and heading error θ in radians. Its curvatures are 1 and 100, so κ = 100 even though nothing is wrong with the model. Measuring heading in a new unit, θ′ = 10θ, gives ½(e² + θ′²), a round bowl that one step solves. Day 11's Newton method does this kind of rescaling automatically, using the Hessian.

When should you stop? You cannot measure f − f* directly, but you can measure the gradient. For m-strongly convex f, ‖∇f(x)‖₂ ≤ √(2mε) guarantees f(x) − f* ≤ ε. With m = 1 and ε = 10⁻⁶, stop once ‖∇f‖ ≤ 1.41 × 10⁻³. The catch is that you usually do not know m. Boyd & Vandenberghe slides

5. Can we do better?

Picture a speed-limit sign for a whole family of methods. A first-order method learns about f only through gradients. Each new point lies in x₀ plus the span of the gradients seen so far. Gradient descent is one such method, but there are cleverer ones.

Nesterov proved a lower bound for this whole family. For any k ≤ (n − 1)/2 and any start, there is an L-smooth convex function on ℝⁿ where every first-order method has CO4A 6-3-6 Tibshirani 2016

f(x_k) − f* ≥ 3L‖x₀ − x*‖₂² / (32(k + 1)²)

So no first-order method can beat a rate of 1/k² on the worst functions. Gradient descent only guarantees 1/k, which leaves a gap. Put in our numbers, L = 10 and ‖x₀ − x*‖² = 101. At k = 10 the unavoidable error is 0.78, while gradient descent's guarantee is 50.5. The lower bound is about hard functions in many dimensions, not about our 2-D bowl.

The gap is real. Accelerated gradient methods reach O(1/k²), that is O(1/√ε) iterations. For ε = 10⁻⁶ that is on the order of a thousand steps instead of a million. That is Day 8 (see the FISTA card). Two other limits of plain gradient descent point to other days. It needs a gradient, so kinks like |x| need Day 7's subgradients. And it ignores curvature, which Day 11's Newton method uses.

6. Gradient boosting: gradient descent on predictions

Picture a line of tutors. Each one looks only at the mistakes the earlier tutors left behind and tries to fix those. Boosting builds simple models one after another, each fitting what the current team still gets wrong. Bagging, by contrast, trains many models independently on resampled data and averages them. CO4A 6-4

Boosting is gradient descent in disguise. Treat the vector of predictions u = (u₁, …, u_n) on the n training points as the variable, with total loss Σ_i L(y_i, u_i). A gradient step on u would move along d_i = −∂L(y_i, u_i)/∂u_i. But the model must also predict new points, so you fit a small tree T_k to d by least squares and step along the tree instead: Tibshirani 2016

u^(k) = u^(k−1) + α_k·T_k

Here α_k is the boosting step size (not the backtracking α). The final model is a weighted sum of trees, which is what gradient boosting is. Take two points with labels y = (3, 5), a first guess u⁽⁰⁾ = (4, 4), and squared loss ½(y − u)². Then d = y − u = (−1, 1), which is just the residuals. A one-split tree that separates the two points matches d exactly. With α = 1 you get u⁽¹⁾ = (3, 5), a perfect fit. For squared loss, "fit the residuals" and "follow the negative gradient" are the same move. Other losses, such as the logistic loss, only change d.

7. Stochastic gradient descent: why machine learning runs on it

Imagine asking for directions. You could survey the whole town before every step, or you could ask one random passer-by. One answer is noisy, but you get a thousand of them for the price of one survey.

Most learning objectives are sums with one term per training example: f(x) = Σ_i=1..m f_i(x). A batch method such as plain gradient descent needs all m gradients for every step. Stochastic gradient descent (SGD) uses just one: Tibshirani 2016 CO4A 6-5

x_k+1 = x_k − t_k∇f_ik(x_k)

The index i_k follows the cyclic rule (1, 2, …, m, 1, 2, …) or the randomized rule (uniform at random), and the random one is more common. In cost, m stochastic steps equal about one batch step. When there are 500 million examples, that difference decides whether training is possible at all. The rule of thumb is that SGD thrives far from the optimum and struggles close to it. Tibshirani 2016 (SG)

Why does it struggle near the end? Near x* the individual gradients ∇f_i disagree, even though their sum is zero. With a fixed step each one keeps kicking the iterate, so it hovers in a cloud around x*. The theory matches: a fixed step only promises to get within a band around f* whose width grows with t. Steps that shrink, with Σt_k = ∞ but Σt_k² < ∞ (for example t_k = 1/k), close the band. These are called diminishing step sizes or step-size decay. Tibshirani 2016 (SG)

A tiny example shows both behaviours. Let f(x) = Σ ½(x − a_i)² with a = (1, 2, 3, 6), so x* = 3, the mean. Batch gradient descent with t = 0.25 from 0 lands on 3 in one step. Cyclic SGD with a fixed t = 0.25 visits 0.25, 0.69, 1.27, 2.45 in its first pass. After 200 passes it is still circling: 2.94, 2.70, 2.78, 3.58. With t_k = 1/k, SGD computes a running mean (1, 1.5, 2, 3) and lands on 3 after one pass.

In practice you use a minibatch: average the gradients of b random examples per step. Averaging cuts the noise, and a GPU computes the b gradients in parallel. One epoch is one pass over the data, so each epoch takes m/b steps. Adam is a stochastic gradient method that also adapts a separate step size for each coordinate. Kingma & Ba 2014 As of Oct 2026 it is a common choice for training neural networks. See the SGD card for a summary.

Try it: compare the blue full-gradient path with the orange SGD scribble. Then switch to b = 10, and then to the decaying step. Watch how low the orange band sits at the end, and when blue passes it.

Notice (defaults: b = 1, fixed t = 0.1): after one epoch full GD has taken a single step and still has f − f* = 0.84. In the same work SGD took 1,000 cheap steps and is already in its noise band; at that instant it happens to dip to 6.6×10⁻⁵. It never leaves the band, which has a median of 8.5×10⁻³ in the last epoch. GD keeps falling in a straight line, passes the band at epoch 32 and ends at 5.6×10⁻⁴. With b = 10 or the decaying step, the band drops below 10⁻³.

1. Why kinks need a new tool

Many objectives worth minimizing have corners. The absolute value |x| has one at 0. The L1 norm ‖x‖₁ = |x₁| + … + |x_n| has one wherever a coordinate is 0. A maximum of several functions has one wherever two of them tie.

These corners are not an accident. A path planner that cares about the worst of several obstacle clearances is working with a maximum. A lane-line fit that should shrug off one bad GPS fix uses absolute errors, because a squared error lets a single outlier dominate. A model that should use only a few of many sensors adds an L1 penalty, and Day 8 builds a fast solver for exactly that.

On Day 6 every function was smooth. For a smooth convex function, the tangent line at any point never rises above the graph Tibshirani 2016:

f(y) ≥ f(x) + ∇f(x)^T(y − x) for all y

Check it on f(x) = x² at x = 1. The tangent is 1 + 2(y − 1) = 2y − 1. The gap is y² − (2y − 1) = (y − 1)², which is never negative. So the tangent stays under the bowl everywhere.

Now look at the tip of a V. Coming from the left the slope is −1, coming from the right it is +1, so there is no single tangent. But many lines still pass through the tip and stay under the V. That is today's whole idea: keep the "line under the graph" property, and drop the demand that the line be unique.

2. Subgradients: every line that fits under the graph

The right-hand picture above is the definition. Each line in the fan touches the V at its tip and never rises above it. The slope of any such line is called a subgradient.

In symbols: a vector g is a subgradient of f at x if the line (or plane) with slope g through the point (x, f(x)) stays under the graph everywhere CO4A 7-1:

f(y) ≥ f(x) + g^T(y − x) for all y

It is the tangent inequality from Section 1, with ∇f(x) replaced by g. A kink is a point where f is continuous but has no derivative, like the tip of the V.

Try it on f(x) = |x| at x = 0. With g = 0.5 the line is 0.5y, and |y| ≥ 0.5y for every y, so 0.5 is a subgradient. With g = 1.2 the line at y = 1 gives 1.2, but |1| = 1. The line pokes above the graph, so 1.2 fails. Any g from −1 to 1 works.

Three facts make subgradients usable. A convex function has at least one subgradient at every point inside its domain. Where a convex f is differentiable, its only subgradient is ∇f(x), so nothing changes for smooth functions. A nonconvex function may have none at all, because no line fits under a dip CO4A 7-2.

Here are the subgradients you will use most often Tibshirani 2016:

Try it: the widget draws the fan at a point x₀. Drag "Your slope g" to 1.2 and watch the line turn red where it rises above f. Then switch to function (c) and look at the last tile. Move x₀ off the kink to see the fan collapse to one slope.

Notice: at the kink x₀ = 0 of |x|, the fan fills every slope from −1 to 1. The default line g = 0.5 sits inside the fan and stays under the V. And because 0 is in ∂f(0) = [−1, 1], the widget reports that x₀ = 0 is a minimizer.

3. The subdifferential, and the rules for building it

Collect every subgradient at x into one set. That set is the subdifferential ∂f(x), the whole fan of allowed slopes CO4A 7-2:

∂f(x) = { g : f(y) ≥ f(x) + g^T(y − x) for all y }

The set is always closed and convex, so in one dimension it is an interval. For convex f it is never empty. It shrinks to the single point {∇f(x)} exactly where f is differentiable.

Take f(x) = max(1 − x, 2x − 2). The two pieces cross at x = 1, where f = 0. The left piece has slope −1 and the right piece has slope 2, so ∂f(1) = [−1, 2]. At x = 0 only 1 − x is active, so ∂f(0) = {−1}: one slope, as at any smooth point.

You rarely check the definition by hand. As with derivatives, a few rules build ∂f from simple pieces CO4A 7-2-2:

• Scaling: ∂(af) = a·∂f for a > 0.
• Sum: ∂(f₁ + f₂) = ∂f₁ + ∂f₂. Add every member of one set to every member of the other.
• Affine inside: if h(x) = f(Ax + b), then ∂h(x) = A^T∂f(Ax + b).
• Maximum: for f = max_i f_i, ∂f(x) is the convex hull of the subgradients of the active pieces, meaning all weighted averages of them.

The L1 norm is a sum of |x_i| terms, so its subgradient is built one coordinate at a time. At x = (2, 0, −1) the first and third coordinates are fixed at sign(2) = 1 and sign(−1) = −1. The middle one is free in [−1, 1]. So ∂‖x‖₁ = {(1, s, −1) : s ∈ [−1, 1]}. That makes (1, 0.3, −1) a subgradient and (1, 1.5, −1) not one.

The sum rule also explains flat-bottomed valleys. Take f(x) = |x| + |x − 2|. At x = 1 both terms are smooth, and ∂f(1) = {1} + {−1} = {0}. At x = 0, ∂f(0) = [−1, 1] + {−1} = [−2, 0], which also contains 0. In fact every x in [0, 2] has 0 in its fan, and f = 2 on that whole stretch.

Constraints fit into the same language. The indicator function I_C of a set C is 0 on C and +∞ outside it. Adding I_C to an objective is the same as demanding x ∈ C, because any point outside costs infinity.

Picture standing with your back against a wall. You can push in any direction that does not point back into the room. That set of pushes is the normal cone, and it is exactly the subdifferential of the indicator at a point x of C CO4A 7-2-1:

∂I_C(x) = N_C(x) = { g : g^T(y − x) ≤ 0 for all y ∈ C }

Take C = [0, 1] on the number line. At x = 0.5, inside, no push is allowed except g = 0, so N_C(0.5) = {0}. At the right end any push to the right is fine: N_C(1) = [0, ∞). At the left end, N_C(0) = (−∞, 0].

One more function you will meet is the distance to a convex set, dist(x, C) = ‖x − P_C(x)‖₂. Here the projection P_C(x) is the point of C closest to x. Outside C the distance is smooth, and its gradient is the unit arrow pointing away from the nearest point CO4A 7-3-5:

∇dist(x, C) = (x − P_C(x)) / ‖x − P_C(x)‖₂

Take C = the unit disk and x = (3, 4). The nearest point is P_C(x) = (0.6, 0.8), the distance is 4, and the gradient is (0.6, 0.8). The formula needs x outside C, where the distance is positive. CO4A does not say so, but the slides do Tibshirani 2016.

4. The optimality rule: 0 ∈ ∂f(x*)

At a minimum you cannot go lower in any direction. In the picture, the flat line through (x*, f(x*)) fits under the whole graph. That is the same as saying slope 0 is a subgradient. This is the subgradient optimality condition CO4A 7-3-1:

f(x*) = min_x f(x) ⇔ 0 ∈ ∂f(x*)

The proof is one line. Put g = 0 into the definition and you get f(y) ≥ f(x*) for all y, which is what "minimum" means. The rule holds for any function. Convexity is what makes ∂f easy to compute, and so makes the rule easy to use.

For max(1 − x, 2x − 2), ∂f(1) = [−1, 2] contains 0, so x* = 1 and f* = 0. Now take max(0.5x, 2x − 1.5). It also has a kink at x = 1, with ∂f(1) = [0.5, 2]. Zero is not in that interval, so this kink is not a minimum. Both slopes are positive, so moving left lowers f. This is function (c) in the widget above.

Constraints work the same way. To minimize a smooth f over C, minimize f + I_C. The sum rule turns 0 ∈ ∂(f + I_C) into a simple test: the downhill direction must point out of C CO4A 7-3-2.

−∇f(x*) ∈ N_C(x*)

Example: minimize (x − 2)² over C = [0, 1]. At x = 1 the gradient is 2(1 − 2) = −2, so −∇f = 2. That lies in N_C(1) = [0, ∞). The ball wants to roll right, the wall at 1 stops it, and x* = 1. Day 10 grows this rule into the KKT conditions.

The rule's most famous use is the lasso: least squares plus an L1 penalty.

minimize ½‖y − Xβ‖₂² + λ‖β‖₁

Here X is the data matrix (one column per feature), y holds the measurements, β holds the coefficients, and λ > 0 sets how hard the penalty pushes. Write X_i for column i. Applying 0 ∈ ∂f one coordinate at a time gives two cases CO4A 7-3-3:

X_i^T(y − Xβ) = λ·sign(β_i) if β_i ≠ 0
|X_i^T(y − Xβ)| ≤ λ if β_i = 0

Read X_i^T(y − Xβ) as the pull the data exert on coefficient i. A coefficient stays at exactly zero while that pull is at most λ, like a block held by static friction. Only a pull stronger than λ gets it moving. This is why the lasso keeps a few features and sets the rest to zero.

With X = I the conditions solve in closed form. Each β_i is y_i moved toward 0 by λ, and set to 0 if |y_i| ≤ λ. This map is soft-thresholding S_λ CO4A 7-3-4:

S_λ(y) = y − λ if y > λ
S_λ(y) = 0 if |y| ≤ λ
S_λ(y) = y + λ if y < −λ

Take y = (3, 0.4, −1.5) and λ = 1, so S₁(y) = (2, 0, −0.5). Check coordinate 1: 3 − 2 = 1 = λ·(+1). Coordinate 2: |0.4 − 0| = 0.4 ≤ 1. Coordinate 3: −1.5 − (−0.5) = −1 = λ·(−1). All three conditions hold. Day 8 turns soft-thresholding into a fast solver for any X.

5. The subgradient method

Now the algorithm. Take gradient descent and replace the gradient by any subgradient. That is the whole subgradient method CO4A 8-1:

x_k+1 = x_k − t_kg_k, g_k ∈ ∂f(x_k)

Here t_k > 0 is the step size at step k. Any member of the fan will do, so you never have to compute the whole set.

There is one catch. A negative subgradient need not point downhill. Near a kink a step can jump to the far side and raise f. So this is not a descent method, and you keep the best point found so far. Its value is f_best:

f_best(k) = min{ f(x₀), f(x₁), …, f(x_k) }

Think of walking down a faceted rock in fog. Each step follows the slope of the face under your feet. Sometimes you step over a ridge onto a face that tilts the other way.

A tiny example shows the catch. Take f(x) = |x|, start at x₀ = 1, and use a fixed step t = 0.3. While x > 0 the subgradient is 1, so the iterates are 0.7, 0.4, 0.1. The next step jumps past 0 to −0.2, where the subgradient is −1. Then it goes back to 0.1, then −0.2 again, forever. The best value stays at f_best = 0.1.

In code the method is a few lines. The only additions to gradient descent are the subgradient and the bookkeeping for the best point.

import numpy as np

def subgradient_method(f, subgrad, x0, step, iters=1000):
    x = np.array(x0, float)
    best_x, best_f = x.copy(), f(x)
    for k in range(1, iters + 1):
        g = subgrad(x)                    # any member of ∂f(x)
        x = x - step(k, x, g) * g
        if f(x) < best_f:                 # not a descent method,
            best_x, best_f = x.copy(), f(x)   # so keep the best point
    return best_x, best_f

f  = lambda x: abs(x[0]) + 2 * abs(x[1])
sg = lambda x: np.array([np.sign(x[0]), 2 * np.sign(x[1])])
print(subgradient_method(f, sg, [2, -1.5], lambda k, x, g: 1 / k)[1])  # ≈ 1.0e-3

Using np.sign(0) = 0 at a kink is fine, because 0 is in [−1, 1]. The subgradient method card collects the update, the rates and the pitfalls in one place.

6. Step sizes, and how fast it can go

The step-size rule decides everything. Unlike on Day 6, the steps are fixed before the run starts. No line search is used CO4A 8-1-1.

• Fixed step: t_k = t for every k.
• Diminishing steps: t_k → 0 with Σt_k = ∞ and Σt_k² < ∞. The standard example is t_k = 1/k.

Why those two sums? Σt_k = ∞ means the steps add up to enough distance to reach any point. Σt_k² < ∞ means the bouncing around the minimum shrinks to nothing. For t_k = 1/k, 1 + 1/2 + 1/3 + … grows without bound, while 1 + 1/4 + 1/9 + … = π²/6 ≈ 1.645.

The guarantees assume f is convex with Lipschitz constant G: |f(x) − f(y)| ≤ G‖x − y‖₂. In words, f never changes faster than G per unit of distance. Then every subgradient has length at most G CO4A 8-1-3.

One master bound covers every rule. Let R = ‖x₀ − x*‖₂ be the distance from the start to a minimizer. After k steps the basic inequality says CO4A 8-1-2:

f_best(k) − f* ≤ (R² + G²Σt_i²) / (2Σt_i)

With a fixed t the bound becomes R²/(2kt) + G²t/2. The first term dies away, but the second does not. So the method can stall up to about G²t/2 above f*. With diminishing steps both sums grow, their ratio goes to 0, and f_best → f*.

How many steps for accuracy ε? Choose t = ε/G², and the bound drops below ε once k ≥ R²G²/ε². With R = G = 1 and ε = 0.01 that is 10,000 iterations CO4A 8-1-4. Gradient descent on an L-smooth function needs O(1/ε) iterations, roughly 100 here up to constants. That gap, O(1/ε²) against O(1/ε), is the price of kinks.

If you happen to know the optimal value f*, there is a smarter step. The Polyak step picks the t that makes the one-step bound on the distance to x* as small as possible CO4A 8-1-6:

t_k = (f(x_k) − f*) / ‖g_k‖₂²

Example: f(x) = |x₁| + 2|x₂| at x₀ = (2, −1.5), where f* = 0. Here f = 5, g = (1, −2) and ‖g‖² = 5, so t = 1. The step lands at (1, 0.5), where f = 2. In the worst case the Polyak step is still O(1/ε²); it just uses information the other rules lack.

Try it: the widget runs all three rules on f = |x₁| + 2|x₂| from (2, −1.5), where G = √5 ≈ 2.236. Switch the fixed step to 0.5 and to 0.01. A bigger step stalls higher, and a smaller one starts slower.

Notice: the fixed step t = 0.1 stalls at f_best = 0.2 from about k = 20 on, under its ceiling G²t/2 = 0.25. The rule t_k = 1/k keeps improving and reaches about 1.0e−3 at k = 1000. The Polyak step is below 1e−10 by k = 50, but that comes from this sharp test function, which grows linearly away from its minimum. It is not the general rate.

7. Constraints: project after each step

Often x must stay in a set C: a box of actuator limits, the nonnegative orthant, a ball. The projected subgradient method takes the usual step, then snaps back to the nearest point of C CO4A 8-1-8:

x_k+1 = P_C(x_k − t_kg_k)

It keeps the same guarantees, provided the projection is cheap. For C = ℝ²₊, the nonnegative quadrant, a raw step that lands at (1, −2) projects to (1, 0): just clip the negative entry to 0. Clipping also projects onto a box, such as limits on steering angle and acceleration.

Easy sets to project onto include {x : Ax = b}, affine images {Ax + b}, the nonnegative orthant, L1, L2 and L∞ balls, and simple cones. A general polyhedron {x : Ax ≤ b} is hard: projecting onto it is itself a quadratic program Tibshirani 2016.

A neat special case finds a point in the intersection of convex sets C₁, …, C_m. Minimize f(x) = max_i dist(x, C_i). If the sets meet, f* = 0, so the Polyak step is available. Working it out, each Polyak step simply projects onto the farthest set. For two sets this is alternating projections CO4A 8-1-7. CO4A's closing sentence says "closest"; that is a typo, and the slides say farthest Tibshirani 2016.

Here are real numbers. Let C₁ be the line x₂ = 0 and C₂ the disk of radius 1.5 centred at (3, 1). Start at (0, 3). C₁ is 3 away and C₂ is 2.106 away, so project onto C₁: (0, 0). Then C₂ is farther, so project onto it: (1.577, 0.526). Then back onto C₁: (1.577, 0). The worst distance falls 1.662, 0.526, 0.239, 0.138, 0.083, … toward the meeting point (≈ 1.882, 0).

1. Smooth plus simple: composite objectives

Yesterday's subgradient method treats the whole objective as kinky, and it pays for that with O(1/ε²) iterations. But most kinky objectives are kinky only in a simple, known way. The lasso is a smooth least-squares fit plus a penalty λ‖β‖₁, whose kinks all sit where some coordinate is 0.

A composite objective is a sum of two such parts CO4A 9-1:

f(x) = g(x) + h(x)

Here g is convex and differentiable, and h is convex, possibly not differentiable, but simple. The lasso splits as g(β) = ½‖y − Xβ‖₂² and h(β) = λ‖β‖₁.

A courier drives the highway to the right block, then walks the last bit on foot. Nobody tries to drive through the pedestrian zone. Today's method works the same way: a gradient step handles the smooth part, and a short exact "walk" handles the kinks.

Sparse sensor selection is a lasso. Say a perception stack has 40 candidate input channels, such as radar range bins or lidar features. You want a linear model that uses only the few that matter. The L1 penalty drives useless weights to exactly 0, so the pattern of zeros tells you which channels to keep. The widget in Section 7 solves a small version: 20 measurements, 40 correlated channels, and only 3 that truly matter.

2. Gradient descent as "jump to the bottom of a simple model"

Here is a new way to see a gradient step. At the current point x, replace f by a simple model: its tangent plane plus a round bowl. Then jump to the bottom of the bowl Tibshirani 2016.

x⁺ = argmin_z f(x) + ∇f(x)^T(z − x) + (1/2t)‖z − x‖₂²

The bowl term (1/2t)‖z − x‖² has Hessian (1/t)I. It pretends the curvature is 1/t in every direction. CO4A writes (1/2t)I here, but the ½ already sits in the bowl term, so the Hessian is (1/t)I. Setting the derivative in z to zero gives ∇f(x) + (z − x)/t = 0, so z = x − t∇f(x). That is exactly the Day 6 step.

A number check: f(x) = x², x = 1, t = 0.25. The model is 1 + 2(z − 1) + 2(z − 1)². Its bottom is at z = 1 − 0.25·2 = 0.5, the same point a gradient step reaches.

Now the key move. Use the simple model only for g, and keep h exactly as it is:

x⁺ = argmin_z g(x) + ∇g(x)^T(z − x)
+ (1/2t)‖z − x‖² + h(z)

Completing the square tidies this up CO4A 9-1:

x⁺ = argmin_z (1/2t)‖z − (x − t∇g(x))‖² + h(z)

Read it as two moves. First take a gradient step on g. Then look for the point that stays close to where that step landed while keeping h small.

3. The proximal operator

That second move has a name. The proximal operator of h with step t is CO4A 9-1:

prox_t(x) = argmin_z (1/2t)‖x − z‖₂² + h(z)

It is a compromise. The first term says "stay near x", the second says "make h small", and a larger t gives h more say. Three examples show the range:

• h = 0: nothing to fix, so prox_t(x) = x.
• h = I_C, the indicator from Day 7: stay as close to x as possible while staying in C. That is the projection P_C(x).
• h = λ‖x‖₁: soft-thresholding at level λt, the subject of Section 4.

Proximal gradient descent repeats the two moves CO4A 9-1:

x_k+1 = prox_t(x_k − t∇g(x_k))

The prox depends only on h, not on g. So g can be as messy as you like, as long as you can compute its gradient. When h's prox is cheap, each iteration costs about as much as a gradient step Tibshirani 2016.

It helps to write the update as x_k+1 = x_k − tG_t(x_k), where the generalized gradient is

G_t(x) = (x − prox_t(x − t∇g(x))) / t

When h = 0 it is just ∇g. In general it plays the gradient's role: G_t(x) = 0 exactly when x minimizes f.

4. ISTA: proximal gradient for the lasso

For h = λ‖β‖₁ the prox splits into one tiny problem per coordinate. Each is solved by Day 7's soft-thresholding, at level λt instead of λ:

[S_λt(β)]_i = β_i − λt if β_i > λt
[S_λt(β)]_i = 0 if −λt ≤ β_i ≤ λt
[S_λt(β)]_i = β_i + λt if β_i < −λt

CO4A drops the minus sign in the middle case and writes λt ≤ β_i ≤ λt. The band is symmetric: −λt ≤ β_i ≤ λt Tibshirani 2016.

Think of a flat tax of λt on every coefficient's size. Big coefficients pay it and survive, a little smaller. Anything smaller than the tax is wiped out to exactly zero.

The gradient of g(β) = ½‖y − Xβ‖² is −X^T(y − Xβ). Put the gradient step and the prox together and you get ISTA, the iterative soft-thresholding algorithm CO4A 9-1 Beck & Teboulle 2009:

β⁺ = S_λt(β + tX^T(y − Xβ))

Try one step on a lasso with a single coefficient: g(β) = ½(4 − 2β)² and λ = 2. Here g′(β) = −2(4 − 2β) and L = 4, so take t = 1/4. Start at β = 0.

• Gradient step: 0 − 0.25·(−8) = 2.
• Prox: S_0.5(2) = 1.5, since λt = 2 × 0.25 = 0.5.
• Check Day 7's rule: g′(1.5) = −2(4 − 3) = −2, and −2 + λ·sign(1.5) = 0. So β = 1.5 is optimal.

It lands in one step only because, in one dimension with t = 1/L, the model of g is exact. In many dimensions ISTA needs many steps. Here is one prox with λt = 1: a gradient step that lands at (3, 0.4, −1.5) is soft-thresholded to (2, 0, −0.5). The middle coefficient falls inside the band and is zeroed.

Try it: drag λ and watch the flat band widen. Then drag the input y through the band and watch the output stick at 0.

Notice: at λ = 1 the input y = 2 comes out as S(2) = 1. Every input in the shaded band [−1, 1], such as 0.5, comes out as exactly 0, and −2 comes out as −1. Outside the band the curve has slope 1, shifted toward 0 by λ.

The whole solver fits in a dozen lines. The same loop gives FISTA (Section 7) when the momentum weight is switched on.

import numpy as np

def soft(z, a):                         # prox of a·‖·‖₁
    return np.sign(z) * np.maximum(np.abs(z) - a, 0)

def lasso(X, y, lam, iters=1000, fista=True):
    t = 1 / np.linalg.norm(X, 2) ** 2   # t = 1/L, L = λmax(XᵀX)
    b = b_old = np.zeros(X.shape[1])
    for k in range(1, iters + 1):
        w = (k - 2) / (k + 1) if fista else 0.0
        v = b + w * (b - b_old)          # momentum (FISTA only)
        b_old = b
        b = soft(v + t * X.T @ (y - X @ v), lam * t)
    return b

The ISTA card and the FISTA card summarize both methods.

5. How fast is ISTA, and how to pick t

The rate matches gradient descent. Assume g has a Lipschitz gradient with constant L, also called L-smooth: ‖∇g(x) − ∇g(y)‖₂ ≤ L‖x − y‖₂. For the lasso, L is the largest eigenvalue of X^TX. With a fixed step t ≤ 1/L CO4A 9-2:

f(x_k) − f* ≤ ‖x₀ − x*‖₂² / (2tk)

That is O(1/k) error, or O(1/ε) iterations, against O(1/ε²) for the subgradient method. Take ‖x₀ − x*‖² = 1 and t = 1. Reaching ε = 10⁻⁴ needs k ≥ 1/(2·10⁻⁴) = 5,000 iterations by this bound. The subgradient bound with R = G = 1 would ask for R²G²/ε² = 10⁸.

Two cautions. The count is iterations, not work. Each iteration costs one prox, which is cheap for soft-thresholding and expensive when it needs a matrix factorization (Section 6). Also, ISTA's objective never goes up: it is a descent method Beck & Teboulle 2009.

You rarely know L. Backtracking line search works much as on Day 6. Start each iteration at t = 1 and shrink t ← βt, with 0 < β < 1, until the smooth part's quadratic upper model holds at the new point CO4A 9-2:

g(x − tG_t(x)) ≤ g(x) − t∇g(x)^TG_t(x)
+ (t/2)‖G_t(x)‖²

Only g is tested, never h. The rate becomes ‖x₀ − x*‖²/(2t_mink) with t_min = min{1, β/L}. Every trial step costs one prox call.

On the seeded lasso from Section 7, computed offline with β = 0.5, 200 backtracking iterations used 1,319 prox calls, about 6.6 per iteration. They accepted steps as large as 0.25, where 1/L = 0.00395, and reached an error of 1.8×10⁻⁴. Fixed-step ISTA was at 9.1×10⁻² after 200 iterations and at 4.9×10⁻⁴ after 1,319. So here the bigger steps were worth their extra prox calls.

6. Special cases: one method, many names

Proximal gradient is also called generalized gradient descent, because familiar methods fall out of it when you choose g or h CO4A 9-4:

• h = 0: the prox is the identity, and you get plain gradient descent.
• h = I_C: the prox is the projection, and you get projected gradient descent, x⁺ = P_C(x − t∇g(x)).
• g = 0: you get the proximal minimization algorithm, x⁺ = prox_t(x). It beats the subgradient method, but only when the prox has a closed form.

A projected-gradient example with numbers: minimize (x − 3)² over C = [0, 1], with L = 2 and t = 0.25, from x = 0. The gradient step goes to 0 − 0.25·(−6) = 1.5, and the projection brings it back to P_C(1.5) = 1. From x = 1 the gradient step goes to 1 − 0.25·(−4) = 2, which projects back to 1 again. The point has stopped moving, and a point the step leaves unchanged is the optimum.

Robots use projected gradient all the time. If the variables are steering and acceleration commands with hard limits, C is a box and P_C just clips each command to its range. Box limits on inputs come back in Day 17's model predictive control.

One heavier example shows why the cost of a prox matters. In matrix completion you see some entries of a matrix Y, such as a few ratings, and want a low-rank B that matches them. The trace norm ‖B‖_tr, the sum of the singular values, is the L1 norm of the singular values. It is a convex stand-in for rank CO4A 9-3.

Its prox soft-thresholds the singular values. Singular values (5, 2, 0.5) with λt = 1 become (4, 1, 0), and the rank drops from 3 to 2. With t = 1 the method is called Soft-Impute: fill the missing entries with your current guess, then shrink the singular values. Each step needs a singular value decomposition (SVD), which is far more work than a soft-threshold.

7. Acceleration: FISTA

ISTA forgets its past. Each step starts from where it stands, so along a long, gently sloping valley it inches forward. Momentum fixes that: before each step, keep going a little in the direction you were already moving.

Think of pushing a swing. Small pushes timed with the motion build speed, because each push adds to the motion already there.

FISTA (fast ISTA) adds one line to ISTA. Set x₋₁ = x₀. Then for k = 1, 2, … CO4A 9-5-1:

v = x_k−1 + ((k − 2)/(k + 1))(x_k−1 − x_k−2)
x_k = prox_t(v − t∇g(v))

The momentum weight (k − 2)/(k + 1) starts small and creeps toward 1. It is 0 at k = 2, 0.25 at k = 3, 0.727 at k = 10 and 0.970 at k = 100. At k = 1 it is −0.5, which does no harm because x₀ − x₋₁ = 0.

This is the form CO4A and the CMU slides use. Beck and Teboulle's original paper builds the weight from a sequence t_k, which there is a momentum counter, not a step size Beck & Teboulle 2009. Start with t₁ = 1 and set

t_k+1 = (1 + √(1 + 4t_k²)) / 2, weight = (t_k − 1)/t_k+1

Both forms carry the same guarantee. On the widget's example they differ only slightly: at k = 100 the error is 6.2×10⁻⁴ with CO4A's weights and 5.6×10⁻⁴ with Beck and Teboulle's.

With t ≤ 1/L, FISTA guarantees CO4A 9-5-2:

f(x_k) − f* ≤ 2‖x₀ − x*‖₂² / (t(k + 1)²)

That is O(1/k²) error, or O(1/√ε) iterations. With ‖x₀ − x*‖² = 1 and t = 1, ε = 10⁻⁴ needs (k + 1)² ≥ 20,000, so k = 141, against 5,000 for ISTA. In the worst case over smooth convex problems, no method that only uses gradients can do better.

Why does looking back help? On a valley floor, consecutive steps point the same way, so the momentum adds up and the effective step grows. Across a valley, steps alternate in direction, and the momentum mostly cancels.

FISTA is not a descent method. Momentum can carry it past the bottom, so f rises for a while. These bumps are called Nesterov ripples CO4A 9-5-1. In the widget's example, FISTA's objective rose on 430 of its first 1,000 steps, the first time at k = 35. ISTA never rose. If you need a curve that only goes down, report the best value so far.

Try it: the chart runs ISTA and FISTA on the seeded lasso (20 measurements, 40 correlated channels, λ = 1, t = 1/L with L = 253.1). Switch the FISTA weights between the two forms, and press "Show bounds" to draw both guarantees.

Notice: on this example, at k = 100 ISTA is still 0.465 above f* while FISTA is at 6.2e−4, about 750 times closer. To get within 10⁻³ of f*, ISTA needs 1,146 iterations and FISTA needs 98.

Acceleration is not always worth it CO4A 9-5-4. A warm start begins a solve from the answer to a nearby problem. Along a path of λ values, λ₁ > λ₂ > …, each solve starts so close to its answer that plain ISTA often does just as well. In matrix completion, backtracking costs several SVDs per step, and the extrapolated point v need not be low rank, which makes each SVD slower.

Prox steps return later. Day 13's ADMM is built from prox steps on each part of a split problem, and Day 14 shows coordinate descent, another fast lasso solver.

1. A referee who proves "you can't do better"

Picture a factory that must ship at least 2 tons of product. Plant X makes it for $1 a ton and plant Y for $3 a ton. Write x and y for the tons each plant makes. The cheapest plan solves this linear program (LP):

minimize x + 3y
s.t. x + y ≥ 2, x ≥ 0, y ≥ 0

Here the answer is obvious: use only the cheap plant, (x, y) = (2, 0), for a cost of 2. With 10,000 variables nothing is obvious. A solver hands you a plan and says it is the best. How could anyone check that claim without trying every plan?

Bring in a referee. She never looks at plans. She only combines the rules that every plan must obey. Take the rule x + y ≥ 2 once and the rule y ≥ 0 twice, and add them: (x + y) + 2y ≥ 2 + 0, which reads x + 3y ≥ 2. The left side is exactly the cost. So no plan can cost less than 2, and the plan (2, 0) is the best one. CO4A 10-01

That is duality in miniature. A floor that matches a plan you hold is a certificate: a short proof of optimality that anyone can check with a few additions. The referee's weights (1, 0 and 2 for the three rules) are the certificate.

Now make the trick systematic. Give each rule a weight: u₁ for x + y ≥ 2, u₂ for x ≥ 0, and u₃ for y ≥ 0. Each weight must be ≥ 0, because multiplying an inequality by a negative number flips it. Adding the weighted rules gives

(u₁ + u₂)x + (u₁ + u₃)y ≥ 2u₁

If the left side equals the cost, that is u₁ + u₂ = 1 and u₁ + u₃ = 3, then 2u₁ is a floor under every plan. The referee wants the highest floor, so she solves a problem of her own: CO4A 10-01

maximize 2u₁
s.t. u₁ + u₂ = 1, u₁ + u₃ = 3, u ≥ 0

Because u₂ and u₃ cannot be negative, these constraints say u₁ ≤ 1 and u₁ ≤ 3. The best choice is u₁ = 1, with floor 2: exactly the cost of the best plan. The referee's problem is the dual problem. The original one is the primal problem.

The weight u₁ also has a price tag. Raise the demand from 2 tons to 2.1 tons, and the cheapest plan costs 2.1 instead of 2. Each extra ton of demand costs u₁ = 1 dollar. A multiplier is the price of its rule: what one more unit of tightening would cost you. Boyd & Vandenberghe 2004 Hold on to that reading. It is the story behind everything below.

2. Every LP has a twin

The referee's recipe works for any LP. Take the form "pay as little as possible while meeting every requirement":

primal: minimize c^Tx
s.t. Ax ≥ b, x ≥ 0

Here x ∈ ℝⁿ lists the amounts, c lists their unit costs, and each row of Ax ≥ b is one requirement. Give each row a weight, and collect the weights in a vector u ≥ 0. Every feasible x then obeys u^TAx ≥ u^Tb. If also A^Tu ≤ c entry by entry, then x ≥ 0 gives c^Tx ≥ (A^Tu)^Tx = u^TAx ≥ b^Tu. The best floor of this kind is the dual:

dual: maximize b^Tu
s.t. A^Tu ≤ c, u ≥ 0

The factory is the case c = (1, 3), A = [1 1] and b = 2. Its dual is "maximize 2u₁ s.t. u₁ ≤ 1, u₁ ≤ 3, u₁ ≥ 0". In Section 1 the extra weights u₂ and u₃ did the same job: u₁ + u₂ = 1 with u₂ ≥ 0 just says u₁ ≤ 1. Tibshirani 2016

Notice the swap. Each primal requirement becomes one dual variable, and each primal variable becomes one dual constraint. The costs c and the requirements b trade places. The table collects the sign rules. CO4A 10-02

A "≤" requirement is multiplied by −1 first. Why may an equality's weight have any sign? You can scale an equation by a negative number without breaking it: a = b also gives −a = −b. CO4A 10-02 states the same rules for the form "Ax = b, Gx ≤ h". Watch the letters on that page: it uses u for the equalities and v for the inequalities, the reverse of chapter 11 and of this lesson.

A dual you can see: max flow and min cut

Picture water pipes running from a source s to a sink t, each pipe with a capacity. The max flow is the most water per second you can push from s to t. A cut splits the nodes into two groups, one holding s and the other holding t. Its capacity is the total capacity of the pipes that run from the s group to the t group. The min cut is the cut with the smallest capacity.

Every drop of water must cross every cut, so any cut's capacity is a ceiling on the flow. The primal here is a maximum, so the referee's job flips. Her certificate now says "you can't push more than this".

Max flow is an LP. Its variables are the flows f_ij in the pipes, with 0 ≤ f_ij ≤ c_ij on each pipe (CO4A prints < here; it should be ≤) and "flow in = flow out" at every node except s and t. Its dual turns out to be the LP relaxation of min cut: the same problem with the 0-or-1 labels "which side is this node on" allowed to take fractional values. So max flow ≤ relaxed min cut ≤ min cut, and the max-flow min-cut theorem says all three are equal. CO4A 10-03 Tibshirani 2016

The figure shows a small case. The best flow sends 4 units. The cut that separates {s, a} from {b, t} crosses three pipes with capacities 2, 1 and 1, so it also totals 4. The other cuts are larger: {s} alone gives 5, {s, b} gives 7 and {s, a, b} gives 5. A flow of 4 and a cut of 4 certify each other, just as the factory plan and the referee's floor did.

3. Prices instead of walls: the Lagrangian

For an LP the referee can combine rules by hand. Curved problems need a recipe that always works. Here it is: turn each hard rule into a price. Instead of a wall that forbids breaking a rule, charge a fine of u for each unit of violation, and pay a rebate of u for each unit of room to spare.

First the letters, because the books disagree. This lesson follows CO4A chapter 11 and Tibshirani's slides.

So in Boyd & Vandenberghe the letter h means an equality, while here it means an inequality. Day 4 followed CO4A chapter 4 and wrote the rules as g_i(x) ≤ 0 and h_j(x) = 0. Today the letter g belongs to the dual function, so the inequalities become h_i and the equalities ℓ_j. Day 10 keeps h_i and ℓ_j but names the multipliers λ and ν, as CO4A chapter 12 does.

Day 4's standard form, in today's letters, is

minimize f(x)
s.t. h_i(x) ≤ 0, i = 1, …, m
ℓ_j(x) = 0, j = 1, …, r

Write p* for its best value. The Lagrangian adds every rule to the cost, each with its own price:

L(x, u, v) = f(x) + Σ_i u_ih_i(x) + Σ_j v_jℓ_j(x)

The prices u_i ≥ 0 and v_j (any sign) are the Lagrange multipliers, also called dual variables. There is one for each rule, just like the referee's weights. CO4A 11-1

One fact makes everything work: on the feasible set, L never exceeds f. If x obeys the rules, each h_i(x) ≤ 0, so u_ih_i(x) ≤ 0. That term is a rebate, or nothing. Each ℓ_j(x) = 0, so the v terms vanish. CO4A 11-1

L(x, u, v) ≤ f(x)
for every feasible x and every u ≥ 0

A worked example. Take min x² s.t. x ≥ 1, whose answer is x* = 1 with p* = 1. In "≤ 0" form the rule is h(x) = 1 − x ≤ 0, so

L(x, u) = x² + u(1 − x)

At the feasible point x = 1.5 with price u = 1, L = 2.25 + (1 − 1.5) = 1.75. That is below f(1.5) = 2.25: the point has 0.5 units of room and earns a rebate of 0.5. At the infeasible point x = 0.5, L = 0.25 + 0.5 = 0.75. That point breaks the rule by 0.5 and pays a fine of 0.5.

4. The dual function g(u, v), and why it is always concave

Now set x free. Fix the prices, ignore the rules, and find the lowest point of L over all x. That lowest value is the Lagrange dual function:

g(u, v) = min over all x of L(x, u, v)

(Strictly, "min" should be "inf", the greatest lower bound, because the lowest value may never be reached.) Why is g a floor? Freedom plus a fair fine can only make things cheaper. Formally, let x* be the best feasible point. Then g(u, v) ≤ L(x*, u, v) ≤ f(x*) = p*. The first step holds because g takes the minimum over all x. The second holds because x* is feasible. CO4A 11-2

Back to min x² s.t. x ≥ 1. For a fixed u, L(x, u) = x² + u(1 − x) is a bowl in x. Its slope 2x − u is zero at x = u/2. Plug that in:

g(u) = u²/4 + u(1 − u/2) = u − u²/4

Every floor is at most p* = 1, as promised. Look at u = 1: the lowest point x = 0.5 breaks the rule. That is fine, because g(u) is a floor, not a plan. The best price, u = 2, gives a floor of exactly 1.

The value of g can also be −∞, which is a useless floor. This happens all the time in LPs. For the factory, L = x + 3y + u₁(2 − x − y) − u₂x − u₃y. Collect terms:

L = (1 − u₁ − u₂)x + (3 − u₁ − u₃)y + 2u₁

This is linear in x and y. Unless both brackets are zero, x or y can run off and make L as negative as you like. So g = 2u₁ when u₁ + u₂ = 1 and u₁ + u₃ = 3, and −∞ otherwise. The general recipe has rebuilt the referee's LP from Section 1. CO4A 10-04 Prices with u ≥ 0 and g(u, v) > −∞ are called dual feasible.

Now the surprising part. Fix one x and watch L(x, u) = f(x) + u·h(x) as the price u changes. It is a straight line in u that starts at height f(x) with slope h(x). Each x gives one line, and g(u) is the lowest line at each u: their lower envelope. The figure draws five of these lines for our example.

The lowest of many straight lines always bends down, never up. It is a concave function: a cap shape, the negative of a convex function. Check it with numbers. Halfway between u = 0 and u = 2, the chord sits at (0 + 1)/2 = 0.5, while g(1) = 0.75 sits above it. Nothing in this argument used the shape of f or h. So g is concave for every problem, however nonconvex the primal is. CO4A 11-3

The dual problem asks for the best floor. Write d* for its best value:

d* = maximize g(u, v) s.t. u ≥ 0

Maximizing a concave function over the convex set u ≥ 0 is a convex problem, because maximizing a concave g is the same as minimizing the convex −g. So the dual is always a convex problem, even when the primal is not. There is one catch. Each value of g needs a minimization over x, and that can be hard when f is nonconvex. CO4A 11-3

5. Weak duality, strong duality, and Slater's condition

Back to the referee. Her best floor can never pass the true answer. This is weak duality: CO4A 11-3

d* ≤ p* (always)

The proof is the one line from Section 4: every g(u, v) is at most p*, so the best one is too. It holds for every problem, convex or not. The difference p* − d* ≥ 0 is the optimal duality gap. When it is zero we have strong duality, p* = d*: the referee's best floor touches the answer.

In the x² example, d* = g(2) = 1 = p*, so there is no gap. The best price u* = 2 also carries the meaning from Section 1. Tighten the rule to x ≥ 1.1, and the best cost rises from 1 to 1.21. That is an increase of 0.21, close to u* × 0.1 = 0.2. The multiplier is the price of the rule. Boyd & Vandenberghe 2004

When is the gap zero? The standard test is Slater's condition. Suppose the problem is convex: f and every h_i are convex, and every ℓ_j is affine. Suppose also that some point is strictly feasible: h_i(x) < 0 for every i, and ℓ_j(x) = 0 for every j. Then p* = d*. Affine inequalities may even hold with equality; only the curved rules need room to spare. CO4A 11-4 Tibshirani 2016

In words: a convex problem with a little wiggle room inside every curved rule has no gap. For x² with x ≥ 1, the point x = 2 has h(2) = 1 − 2 = −1 < 0. Slater holds, and the gap is indeed zero.

An LP has only affine rules, so for LPs feasibility is enough. Strong duality holds whenever the primal or the dual is feasible, and it fails only if both are infeasible. Tibshirani 2016 That is why the max flow equals the relaxed min cut in Section 2.

Without convexity, a real gap can open. Take f(x) = x⁴ − 50x² + 100x s.t. x ≥ −4.5. The cost has two valleys, near x = −5.4 and x = 4.4. The rule cuts off the deeper one, so the best feasible point sits on the boundary: p* = f(−4.5) = −1052.4375. CO4A 11-3

The Lagrangian is L = x⁴ − 50x² + (100 − u)x − 4.5u. The best price is u = 100. It cancels the linear term and leaves two equal valleys at x = ±5, so g(100) = 625 − 1250 − 450 = −1075. No other price does better, so d* = −1075 and the gap is 22.5625. Note that strict feasibility holds here: x = 0 is strictly inside the rule. What fails is convexity. CO4A 11-3

The gap is also useful while a solver runs. For any feasible x and any dual-feasible (u, v), the duality gap is f(x) − g(u, v). Since g(u, v) ≤ p* ≤ f(x), it bounds how far x is from the best value: CO4A 11-5

f(x) − p* ≤ f(x) − g(u, v)

For example, x = 1.2 costs f = 1.44, and the price u = 1.5 gives g = 1.5 − 0.5625 = 0.9375. The gap is 0.5025, so x is at most 0.5025 above the best value. (The truth is 0.44.) A gap of zero proves that both are optimal. Solvers stop once the gap falls below a tolerance ε, because then f(x) − p* ≤ ε is guaranteed. CO4A 11-5

Try it: drag the price u and watch the dashed priced cost L(x, u) tilt. Its lowest point, the orange dot, gives the floor g(u), the purple dot below. Then switch to the nonconvex quartic and hunt for the best price.

Notice: at the default price u = 1 the floor is g(1) = 0.75, a gap of 0.25, and L's lowest point x = 0.5 breaks the rule x ≥ 1. The purple curve stays under p* = 1 everywhere and touches it only at u = 2 (Slater holds: x = 2 is strictly feasible). On the quartic, the best floor is −1075 at u = 100, where L has two lowest points, x = −5 and x = 5. That floor stops 22.56 short of p* = −1052.44.

6. Worked example: the dual of a least-squares problem

A small robot has two thrusters. Thruster 1 pushes with force x₁. Thruster 2 sits on a lever, so its force x₂ counts double. Together they must deliver a push of exactly 5, so x₁ + 2x₂ = 5. Energy grows with the square of each force, so we want the least total effort:

minimize x₁² + x₂²
s.t. x₁ + 2x₂ = 5

This asks for the shortest solution of a linear equation, a least-squares problem with a constraint. In matrix form it is min ‖x‖₂² s.t. Ax = b, with A = [1 2] and b = 5. There is one equality, so there is one multiplier v, of any sign.

Step 1, the Lagrangian:

L(x, v) = x₁² + x₂² + v(x₁ + 2x₂ − 5)

Step 2, minimize over x. L is a bowl in x. Its gradient is (2x₁ + v, 2x₂ + 2v), which is zero at x₁ = −v/2 and x₂ = −v.

Step 3, plug back in. Then x₁² + x₂² = 1.25v² and x₁ + 2x₂ − 5 = −2.5v − 5, so

g(v) = 1.25v² + v(−2.5v − 5)
= −1.25v² − 5v

This is a downward parabola, concave as promised. Any v gives a floor. For example, v = −1 gives g = −1.25 + 5 = 3.75.

Step 4, maximize. The slope −2.5v − 5 is zero at v* = −2, so d* = g(−2) = −5 + 10 = 5.

Step 5, compare. Put v* = −2 into Step 2 to get x = (1, 2). It obeys 1 + 2·2 = 5, and its effort is 1 + 4 = 5 = d*. The gap is zero, as Slater predicts: the problem is convex and its only rule is affine.

Two things to notice. The multiplier v* = −2 is negative, which is allowed for an equality. And the dual's answer handed back the primal's answer through Step 2. Day 10 turns that into a rule.

The same five steps work for any A and b. The lowest point of L is x = −A^Tv/2, and

g(v) = −¼ v^TAA^Tv − b^Tv

Boyd & Vandenberghe work this exact problem in chapter 5, and CO4A 11-2 runs the same steps on quadratic programs. Boyd & Vandenberghe 2004 CO4A 11-2 The few lines of Python below check every number.

import numpy as np
A = np.array([[1.0, 2.0]]); b = np.array([5.0])      # x1 + 2*x2 = 5
AAt = A @ A.T
x_star = A.T @ np.linalg.solve(AAt, b)               # least-norm solution
g = lambda v: -0.25 * v @ AAt @ v - b @ v            # dual function
v_star = -2 * np.linalg.solve(AAt, b)                # where g peaks
print(x_star, x_star @ x_star)                       # [1. 2.] 5.0
print(v_star, g(v_star), g(np.array([-1.0])))        # [-2.] 5.0 3.75
print(-A.T @ v_star / 2)                             # [1. 2.]: primal from dual

7. Where duality goes next

Today's ideas carry much of the rest of the course. Three later days use them directly.

• Day 10, KKT conditions. When the gap is zero, the best x* and the best prices u* are locked together by four conditions you can check. You have already met one. In the x² example, u* = 2 is positive and the rule is tight at x* = 1. A rule with room to spare would get a price of zero.
• Day 13, dual methods. Instead of solving the primal, climb the concave hill g. In our example the slope of g is g′(u) = 1 − u/2. That equals h(x̂) = 1 − x̂ at L's lowest point x̂ = u/2. So "raise the price while the rule is broken, lower it while there is room" walks uphill on g. That is dual ascent (dual ascent card), and ADMM grows from it (ADMM card).
• Day 19, collision avoidance. A self-driving car must keep a safe distance from every obstacle. The distance between two convex shapes is itself a minimization, which is awkward to place inside a planner's constraints. OBCA writes that distance as a dual problem instead. By weak duality, any dual-feasible prices certify "the distance is at least this much". So the planner carries the prices as extra variables and asks the certified floor to exceed the safe distance. Strong duality means nothing is lost. See the OBCA paper.

1. Four checks at the bottom of the bowl

Picture a marble in a round bowl, with a straight fence across part of it. If the lowest point of the bowl is on the allowed side, the marble rolls there and the fence does nothing. If the lowest point is behind the fence, the marble ends up pressed against it. It stays put because the fence pushes back exactly as hard as the slope pulls.

Three facts hide in that picture. A fence can push but never pull. A fence the marble is not touching pushes with zero force. And at rest, the pull and all the pushes add up to zero. Add "the marble is on the allowed side", and you have the whole of today's lesson.

Now with numbers. Minimize f(x) = (x₁ − 2)² + (x₂ − 1)², a round bowl with its lowest point at (2, 1). Add one fence: x₁ + x₂ ≤ 2. The bowl's bottom breaks this rule, since 2 + 1 = 3, so the answer must lie on the fence line x₁ + x₂ = 2. The figure shows three versions of the problem.

Now the general version. A constrained problem asks to minimize f(x) subject to h_i(x) ≤ 0 for i = 1, …, m and l_j(x) = 0 for j = 1, …, r. As on Day 9, each inequality gets a multiplier and each equality gets a multiplier. Day 9 called them u_i and v_j, like CO4A's chapter 11 and the CMU slides. From today we switch to the more common names λ_i and ν_j. The Lagrangian is

L(x, λ, ν) = f(x) + Σ λ_ih_i(x) + Σ ν_jl_j(x)

Think of λ_i as the force of fence i. The gradient ∇h_i(x) points straight out of the allowed region, at right angles to the fence. A constraint with h_i(x) = 0 is an active constraint: the marble touches that fence.

The KKT conditions are four checks on a candidate (x, λ, ν). They are named after Karush, who wrote them in a 1939 master's thesis, and Kuhn and Tucker, who published them in 1951. CO4A 12

Read stationarity as −∇f = Σ λ_i∇h_i + Σ ν_j∇l_j. The left side is the downhill direction. The right side mixes the outward normals of the walls you touch, with weights λ_i ≥ 0. So at the optimum, every downhill direction leads straight through a wall. If f or h_i has kinks, use the subgradients of Day 7 in place of gradients. CO4A 12-01

Check the middle panel. Here h(x) = x₁ + x₂ − 2, so ∇h = (1, 1), and ∇f(x) = (2(x₁ − 2), 2(x₂ − 1)). Try x⋆ = (1.5, 0.5) with λ = 1. (The star ⋆ marks an optimum; the plain * is saved for section 7.)

• Stationarity: ∇f(x⋆) + λ∇h = (−1, −1) + (1, 1) = (0, 0). ✓
• Primal feasibility: 1.5 + 0.5 − 2 = 0 ≤ 0. ✓
• Dual feasibility: λ = 1 ≥ 0. ✓
• Complementary slackness: λ · h = 1 · 0 = 0. ✓

The cost is f(x⋆) = 0.5. Day 9's dual function here is g(λ) = λ − λ²/2, which peaks at λ = 1 with value 0.5: the duality gap is zero. In the corner panel, ∇f(1.8, 0.2) = (−0.4, −1.6), and the wall normals (1, 1) and (0, 1) cancel it with λ₁ = 0.4 and λ₂ = 1.2. Both walls are doing work.

To find λ by hand, guess which constraints are active, solve stationarity with those constraints as equalities, and check the signs. Here the guess "active" gives x₁ = 2 − λ/2 and x₂ = 1 − λ/2, and x₁ + x₂ = 2 forces λ = 1. Had λ come out negative, the guess would be wrong, and that fence should be released.

The multiplier is also a price. Move the fence out by a small amount ε, and the best cost drops by about λε: with x₁ + x₂ ≤ 2.01 the optimum falls from 0.5 to 0.49005. Boyd & Vandenberghe 2004 For a car planner that keeps a margin from the lane edge, λ is the cost saved per centimetre of extra room.

2. When the checks are enough, and when they are required

Two different questions. If a point passes all four checks, is it optimal? That asks whether KKT is sufficient. If a point is optimal, does it pass? That asks whether KKT is necessary.

Sufficient, for convex problems. Suppose f and every h_i are convex and every l_j is affine (flat, like a line or plane). Then any (x, λ, ν) that passes all four checks is optimal, with zero duality gap. The one-line reason: stationarity says x minimizes the convex function L(·, λ, ν), so g(λ, ν) = L(x, λ, ν). Slackness and feasibility make that equal to f(x). A feasible point whose cost equals a lower bound cannot be beaten. CO4A 12-01 Boyd & Vandenberghe 2004

Necessary, when strong duality holds. Strong duality means the primal and dual optimal values are equal. Then every optimal x⋆, paired with any dual optimal (λ⋆, ν⋆), passes the four checks. Day 9 gave an easy test for strong duality, Slater's condition: the problem is convex and some point is strictly feasible, with every h_i(x) < 0. CO4A 12-01

Put together: for a convex problem that meets Slater's condition, "x is optimal" and "x passes KKT with some λ, ν" mean the same thing. So solving by hand means writing down the four checks and solving them.

Solvers use the same checks, but a computer never hits zero exactly. So it measures how far each check is from holding. These are the KKT residuals, such as the size of the stationarity left side or the amount of constraint violation. For convex problems, small residuals mean "nearly optimal".

The duality gap f(x) − g(λ, ν) is a second certificate, since it is never smaller than the true error f(x) − f⋆. CO4A 13-01 In our example, the feasible guess x = (1, 1) has f = 1, and λ = 0.5 gives g = 0.375. The gap promises an error of at most 0.625; the true error is 0.5.

This is how a model predictive controller (MPC) decides that it is done. A car's MPC solves a fresh quadratic program many times a second (Day 17). The OSQP solver stops when its primal residual ‖Ax − z‖_∞ and dual residual ‖Px + q + A^Ty‖_∞ both drop below a tolerance: 10⁻³, absolute and relative, by default as of Oct 2026. OSQP docs OSQP settings Its iterates meet complementary slackness by construction, so these two residuals are the KKT checks still open. The OSQP paper and the ADMM card show how it gets there. Interior-point solvers (Day 12) instead aim at a relaxed KKT system with λ_ih_i(x) = −1/t and tighten it as t grows. Boyd & Vandenberghe 2004

3. Only equality constraints: KKT is one linear system

Now replace the fences with rails: the point must stay on a line or a flat plane. With no inequalities, checks 3 and 4 have nothing to say. Two remain: the pull is balanced by the rails, and the point is on the rails. For a quadratic cost both are linear equations, so the whole problem is one linear solve.

The simplest picture: find the point on the line x₁ + x₂ = 1 closest to the origin. Your eye says (0.5, 0.5), where the line meets the direction from the origin at a right angle. KKT gets there mechanically.

The general problem is to minimize ½x^TPx + q^Tx + r subject to Ax = b, with P ⪰ 0. Stationarity reads Px + q + A^Tν = 0 and feasibility reads Ax = b. Stacked, they form one square system: CO4A 12-02 Boyd & Vandenberghe 2004

P x + A^Tν = −q
A x = b

[ P A^T ; A 0 ] [ x ; ν ] = [ −q ; b ]

Numbers: minimize ½(x₁² + x₂²) subject to x₁ + x₂ = 1. Here P is the 2×2 identity, q = 0, A = [1 1] and b = 1. The three equations are x₁ + ν = 0, x₂ + ν = 0 and x₁ + x₂ = 1. So x⋆ = (0.5, 0.5) and ν⋆ = −0.5. The multiplier is negative, and that is fine. A rail can push either way, so only inequality multipliers carry a sign.

import numpy as np
P = np.eye(2); q = np.zeros(2)
A = np.array([[1.0, 1.0]]); b = np.array([1.0])
K = np.block([[P, A.T], [A, np.zeros((1, 1))]])   # the KKT matrix
sol = np.linalg.solve(K, np.concatenate([-q, b]))
x, nu = sol[:2], sol[2:]
print(x, nu)                                       # [0.5 0.5] [-0.5]

Two links forward. On Day 11, Newton's method for equality-constrained problems solves this same system at every step, with P = ∇²f(x) and q = ∇f(x). CO4A 12-02 And a linear MPC with a quadratic cost, whose only constraints are the vehicle's linear dynamics, is exactly this kind of problem. One linear solve gives the whole plan. Add a steering limit or a lane bound and you get inequalities, which is why the solvers on Day 17 have to iterate.

4. Water-filling: pour the power into the quiet channels

A radio has four channels, each with its own noise level α_i, and a fixed power budget P to split across them. Quiet channels turn power into data well, but every channel has diminishing returns. How should you split it?

Draw each channel as a vessel whose floor sits at height α_i. Pour the whole budget in as water and let it settle. The water finds one common level, so quiet channels (low floors) get deep water and noisy ones get little or none. The depth in vessel i is its power x_i. KKT says this picture gives the exact optimum.

The problem is Boyd and Vandenberghe's Example 5.2: CO4A 12-03 Boyd & Vandenberghe 2004

min −Σ log(α_i + x_i)
s.t. x ≥ 0, x₁ + … + x_n = P

Give λ_i to the constraint −x_i ≤ 0 and ν to the budget. The four checks become

−1/(α_i + x_i) − λ_i + ν = 0
λ_ix_i = 0, λ_i ≥ 0, x_i ≥ 0, Σ x_i = P

Split into two cases. If channel i gets power (x_i > 0), slackness forces λ_i = 0, so α_i + x_i = 1/ν: its water reaches the level 1/ν. If it gets none (x_i = 0), then λ_i = ν − 1/α_i ≥ 0, which means α_i ≥ 1/ν: its floor is above the level. One formula covers both cases:

x_i = max{0, 1/ν − α_i}

Only one number is still unknown, the water level 1/ν. The budget fixes it through Σ max{0, 1/ν − α_i} = P. The left side rises steadily with the level, so exactly one level works. This is water-filling. CO4A 12-03

Numbers: α = (0.1, 0.4, 0.6, 1.2) and P = 1. Guess that three channels are wet and write L = 1/ν. Then (L − 0.1) + (L − 0.4) + (L − 0.6) = 1, so 3L = 2.1 and L = 0.7. The guess holds up: the third floor (0.6) is below 0.7 and the fourth (1.2) is above. So x = (0.6, 0.3, 0.1, 0) and ν = 1/0.7 ≈ 1.4286.

The dry channel's multiplier is λ₄ = ν − 1/α₄ = 1.4286 − 0.8333 ≈ 0.595. Read it as a price again. Force a small amount of power ε into channel 4, taking it from the wet channels, and Σ log(α_i + x_i) drops by about 0.595ε. At the optimum Σ log(α_i + x_i) = −0.888, against −1.272 for an equal split of 0.25 each.

Try it: drag the budget and watch the channels switch on one by one. Then drag the level by hand. The formula keeps stationarity, dual feasibility and slackness at zero for any level, so only the budget residual moves. "Match budget" snaps it back.

Notice that at budget 1 the level is 0.700 (ν = 1.4286), the powers are 0.600, 0.300, 0.100 and 0, and the dry channel 4 has λ₄ = 0.595. All four KKT residuals read 0.0000. The channels switch on at budgets 0, 0.30, 0.70 and 2.50, as the level reaches each floor.

5. SVM: the points that hold the line

A support vector machine separates two classes of points with the widest possible street. Only the points on the edge of the street, or inside it, matter. Delete a point that sits comfortably on its own side, and the street does not move. Complementary slackness is the reason.

The soft-margin SVM picks a direction β, an offset β₀, and a slack ξ_i for each point, which measures how far point i intrudes into the street: CO4A 12-04

min ½‖β‖₂² + C Σ ξ_i
s.t. y_i(x_i^Tβ + β₀) ≥ 1 − ξ_i, ξ_i ≥ 0

Here x_i is a data point, y_i = ±1 is its class, and C > 0 is the cost of an intrusion. Give the multiplier w_i to the margin constraint and v_i to ξ_i ≥ 0. Stationarity gives

β = Σ w_iy_ix_i, Σ w_iy_i = 0, w_i = C − v_i

Since v_i ≥ 0, every w_i lies between 0 and C. Complementary slackness, w_i(1 − ξ_i − y_i(x_i^Tβ + β₀)) = 0 and v_iξ_i = 0, then sorts the points into three groups:

• Clear of the street (margin y_i(x_i^Tβ + β₀) > 1): the constraint is slack, so w_i = 0, and the point drops out of β.
• On the edge (ξ_i = 0 and margin exactly 1): w_i can be anywhere from 0 to C.
• Inside the street or on the wrong side (ξ_i > 0): then v_i = 0, so w_i = C.

The points with w_i > 0 are the support points, also called support vectors. They alone build β.

Numbers, in one dimension with C = 1: x = −1 (class −1), and x = 1 and x = 3 (class +1). The solution is β = 1 and β₀ = 0, with no intrusions, and w = (0.5, 0.5, 0). Check stationarity: β = 0.5·(−1)·(−1) + 0.5·1·1 + 0·1·3 = 1, and Σ w_iy_i = −0.5 + 0.5 + 0 = 0. The point at 3 has margin 3 > 1, so slackness forces its w to 0. The primal cost ½β² = 0.5 equals the dual value, so the gap is zero.

6. Budget or price: constrained and penalty forms

You can keep a model simple with a budget ("total size at most t") or with a price ("pay λ per unit of size"). A traveller told "spend at most $200" and one fined for every dollar spent can make the same choices, if the fine is set right. KKT says which price matches which budget.

(C) min f(x) s.t. h(x) ≤ t
(L) min f(x) + λ h(x)

The second line is the Lagrange form, also called the penalized form. Going from (C) to (L): if (C) has a strictly feasible point, its KKT multiplier λ⋆ makes x⋆ a solution of (L). The reason is that the stationarity check of (C) is, word for word, the optimality condition of (L). Going back: a solution x⋆ of (L) solves (C) with t = h(x⋆). So as t and λ sweep their ranges, both forms trace out much the same family of solutions. CO4A 12-05 CMU 10-725 2016

Numbers, with a one-dimensional lasso and data y = 3. The penalty form min ½(3 − β)² + 1·|β| gives β = 2: shrink 3 toward zero by λ = 1, the soft-thresholding of Day 8. The budget form min ½(3 − β)² s.t. |β| ≤ 2 also gives β = 2. Its stationarity check at β = 2 reads −(3 − 2) + λ·1 = 0, so its multiplier is λ = 1, the same price.

One more KKT fact, about the lasso min ½‖y − Xβ‖₂² + λ‖β‖₁ with n data points and p features. When p > n it could in principle have many solutions, but KKT shows that they all share the fitted values Xβ̂ and, for λ > 0, the value ‖β̂‖₁. If the entries of X come from a continuous distribution, the solution is unique with probability 1. It also has at most min{n, p} nonzeros, even when p is far larger than n. CO4A 12-06 So a lasso fit on 50 samples of continuous data picks one answer that uses at most 50 features.

7. Writing duals fast: dual norms and conjugates

CO4A chapter 13 collects tools that turn "derive the dual" from a page of algebra into a lookup. Two ideas do most of the work. CO4A 13

Dual norms. How big can z^Tx get if x must stay inside the unit ball of some norm? Picture a line at right angles to z, sliding outward until it last touches the ball. The point it last touches is the best x.

That largest value is the dual norm of z:

‖z‖_* = max { z^Tx : ‖x‖ ≤ 1 }

So |z^Tx| ≤ ‖x‖ · ‖z‖_* for every x and z, a general form of the Cauchy–Schwarz inequality. The pairs to know are ℓ₁ ↔ ℓ_∞, ℓ₂ ↔ ℓ₂, and ℓ_p ↔ ℓ_q with 1/p + 1/q = 1. For matrices, the trace norm (sum of singular values) pairs with the operator norm (largest singular value). The dual of the dual is the original norm. CO4A 13-03 Boyd & Vandenberghe 2004

With z = (3, −4), the diamond's best corner (0, −1) gives 4 = ‖z‖_∞. The square's best corner (1, −1) gives 3 + 4 = 7 = ‖z‖₁. The circle's best point z/5 gives 5 = ‖z‖₂.

A robotics use: suppose a car's position error w is only known to satisfy |w_i| ≤ 0.2 m in each coordinate. Then w can add at most 0.2‖c‖₁ to c^T(x + w), since ℓ₁ is the dual of the box's ℓ_∞. So a robust planner enforces c^Tx ≤ d − 0.2‖c‖₁ in place of c^Tx ≤ d.

Conjugate functions. Picture selling x units at price y each, when making them costs f(x). Your best profit at price y is the conjugate of f. On a graph, yx is a line through the origin with slope y. The best profit is the biggest vertical gap between that line and the curve, and it sits where the curve's slope equals y.

In symbols, the conjugate of f is

f*(y) = max_x [ y^Tx − f(x) ]

f* is always convex, even when f is not, because it is a maximum of functions that are linear in y. Fenchel's inequality, f(x) + f*(y) ≥ x^Ty, holds for every x and y, with equality when y is a slope (subgradient) of f at x. If f is closed and convex, then f** = f. CO4A 13-04 Boyd & Vandenberghe 2004

Numbers with f(x) = x²: the gap yx − x² peaks at x = y/2, so f*(y) = y²/4. Fenchel at x = 1 and y = 2 gives 1 + 1 = 2 = xy, an equality because 2 is the slope of x² at 1. At y = 3 it gives 1 + 2.25 = 3.25 ≥ 3.

The indicator function I_C is 0 on a set C and +∞ off it, which turns a constraint into a function. Its conjugate is the support function of C, max over x ∈ C of y^Tx. Four conjugates cover most needs:

For the quadratic, y^Tx − ½x^TQx is strictly concave in x, and its maximizer is x = Q⁻¹y. CO4A 13-04 copies a slip from the CMU slide and says "concave in y". CO4A 13-04 CMU 10-725 2016

The shifting trick. For min f(x) + g(Ax), name z = Ax and give that constraint a multiplier u. Minimizing the Lagrangian over x and over z separately produces two conjugates: CO4A 13-04-03

dual: max_u −f*(A^Tu) − g*(−u)

With A = I, this splits min f(x) + g(x). When g is a norm, the table turns g*(−u) into the constraint ‖u‖_* ≤ 1. CO4A 13-04-02 After solving the dual, you can read off x⋆ as the minimizer of L(x, u⋆), when that minimizer is unique. CO4A 13-02

The lasso shows the trick at work. Write z = Xβ in min ½‖y − z‖₂² + λ‖β‖₁. Minimizing over z gives z − y + u = 0, so z = y − u. (CO4A 13-04-01 and the CMU slide print "z − y + β = 0"; it should be u.) Minimizing over β leaves the constraint ‖X^Tu‖_∞ ≤ λ, because ℓ_∞ is the dual of ℓ₁. The dual is CO4A 13-04-01

max ½(‖y‖₂² − ‖y − u‖₂²)
s.t. ‖X^Tu‖_∞ ≤ λ

It finds the point of the region {u : ‖X^Tu‖_∞ ≤ λ} closest to y, and that point is the lasso's residual: Xβ⋆ = y − u⋆. With X = 1, y = 3 and λ = 1, the region is |u| ≤ 1, so u⋆ = 1, which is indeed 3 − 2. The dual value ½(9 − 4) = 2.5 matches the primal ½·1² + 1·2 = 2.5.

One subtlety: rewriting that dual as min ‖y − u‖₂² keeps the solution u⋆ = 1 but changes the optimal value to 4. A rewritten dual keeps its solution, not necessarily its value. CO4A 13-06

1. The Newton step: jump to the bottom of a fitted bowl

Gradient descent from Day 6 is like walking downhill in fog: you feel the slope under your feet and take a stride of a set length. Newton's method also feels how the slope is changing, the curvature. Slope plus curvature is enough to fit a bowl to the ground around you, and Newton jumps to the bowl's lowest point.

Start in one dimension. At the current point x_k, measure the slope f′(x_k) and the curvature f″(x_k). The parabola with the same height, slope and curvature is the local quadratic model q(x):

q(x) = f(x_k) + f′(x_k)(x − x_k) + ½f″(x_k)(x − x_k)²

The bottom of this bowl is where q′(x) = 0. Solving that one linear equation gives the 1-D Newton update:

x_k+1 = x_k − f′(x_k) / f″(x_k)

A worked example: f(x) = log(e^x + e^−x), whose minimum is at 0. At x_k = 1 the slope is 0.76 and the curvature is 0.42, so the jump goes to 1 − 0.76/0.42 = −0.81. That overshoots 0 a little. Repeat the recipe and the points go 0.41, −0.047, 0.00007, closing in very fast. Boyd & Vandenberghe 2004

In n dimensions the slope becomes the gradient ∇f(x) and the curvature becomes the Hessian ∇²f(x), the n × n matrix of second derivatives. The quadratic model of f around x is

f(y) ≈ f(x) + ∇f(x)^T(y − x) + ½(y − x)^T∇²f(x)(y − x)

Setting the model's gradient to zero gives the Newton step v, the jump to the bottom of the bowl. CO4A 14-01

∇²f(x) v = −∇f(x), x⁺ = x + v

Written with an inverse, v = −(∇²f(x))⁻¹∇f(x). In code you never form the inverse: you solve the linear system, usually by Cholesky factorisation.

Gradient descent is the same recipe with a lazier bowl: (1/t)I in place of ∇²f, a round bowl whose bottom is x − t∇f(x). Newton's bowl has the right shape, long where the valley is long and steep where it is steep. So gradient descent zig-zags across a narrow valley while Newton heads straight in. CO4A 14-01-01

The cleanest example is a quadratic, f = ½(10x₁² + x₂²) from (1, 1). Its quadratic model is the function itself, so one Newton step lands exactly on (0, 0). Gradient descent needs 106 steps with t = 0.1 to bring f below 10⁻¹⁰, and 62 with the best fixed step, t = 2/11.

A second view: minimising f means solving ∇f(x) = 0, and Newton's method began as a root finder. To solve F(x) = 0, replace F by its tangent and solve that: x⁺ = x − F′(x)⁻¹F(x). With F = ∇f, F′ is the Hessian, so this is the same step. CMU 10-725 2016 CO4A 14-02-01

2. Damped Newton: shorten the jump when the bowl is a poor guide

The bowl is fitted at one point. Near that point it matches f well; far away it can be a poor guide, and a full jump can land somewhere worse than where you started.

The function from section 1 shows this. Its full Newton step is x⁺ = x − ½sinh(2x). From x₀ = 1.1 the points go −1.13, 1.23, −1.70, 5.71, −23,000, each jump overshooting more. Always taking the full step t = 1 is called pure Newton, and here it diverges from every |x₀| > 1.0887. Boyd & Vandenberghe 2004

The fix keeps the direction and shortens the jump. Damped Newton picks the step length t by backtracking line search, as on Day 6: try t = 1, and while the drop in f is not "enough", shrink t.

t = 1; while f(x + tv) > f(x) + αt∇f(x)^Tv: t = βt

Here α ∈ (0, ½] and β ∈ (0, 1). "Enough" means f falls by at least a fraction α of what the slope along v promises. CO4A 14-04 Boyd and Vandenberghe suggest a small α such as 0.01 and β around 0.5. Boyd & Vandenberghe 2004

Back to x₀ = 1.1, now with α = 0.1 and β = 0.7. The first step is shortened once, to t = 0.7, and lands at −0.46. Then every full step passes the test: 0.068, −0.00021, 6·10⁻¹². Four steps, only one damped.

That pattern is general. Near the optimum the full step t = 1 always passes the test (for α < ½), so damped Newton turns into pure Newton by itself and keeps the fast finish. CO4A 14-05

3. The Newton decrement: a built-in progress meter

Before Newton jumps, it already knows how deep its bowl is below the current point. That depth forecasts how much f will drop, and it is almost free once v is computed. The Newton decrement λ(x) measures it:

λ(x)² = ∇f(x)^T (∇²f(x))⁻¹ ∇f(x) = −∇f(x)^Tv

It has three readings. CO4A 14-03

• f(x) minus the bottom of the bowl equals λ²/2. This is the predicted drop.
• The slope of f along the Newton step is ∇f(x)^Tv = −λ². That is the term the backtracking test uses.
• λ is the length of the step v measured with the Hessian as the ruler: λ² = v^T∇²f(x)v.

The first reading gives a stopping rule. You cannot compute f(x) − p⋆ during a run, because the optimal value p⋆ is unknown, but you can compute λ²/2. Stop when it falls below a tolerance ε. Put together, this is Boyd and Vandenberghe's Algorithm 9.5. Boyd & Vandenberghe 2004

given a start x, tolerance eps, alpha in (0, 0.5], beta in (0, 1)
repeat:
    g = grad f(x);  H = hess f(x)
    v = solve(H, -g)                  # Newton step (Cholesky)
    lam2 = -g @ v                     # Newton decrement squared
    if lam2 / 2 <= eps: stop          # stopping rule
    t = 1
    while f(x + t*v) > f(x) - alpha*t*lam2:
        t = beta * t                  # backtracking
    x = x + t*v

How good is the forecast? On the widget's function (section 5), λ²/2 goes 4.45, 0.86, 0.14, 4.7·10⁻³, 5.6·10⁻⁶, 7.9·10⁻¹², against true errors f − p⋆ of 6.60, 1.17, 0.157, 4.8·10⁻³, 5.6·10⁻⁶, 7.9·10⁻¹². Rough at the start, almost exact at the end, which is where a stopping rule needs it.

4. Why it finishes so fast: two phases

Near the answer, Newton's next error is roughly the square of its current error. Square 10⁻³ and you get 10⁻⁶. Square that and you get 10⁻¹². The number of correct digits doubles every step. This is quadratic convergence.

Compare linear convergence, where the error shrinks by a fixed factor each step. Gradient descent on the widget's function shrinks it by about 0.42 per step, one digit every two or three steps. Boyd & Vandenberghe 2004

In symbols, near a minimiser x⋆ there is a constant K with

‖x⁺ − x⋆‖ ≤ K ‖x − x⋆‖²

The one-line reason: the bowl gets the height, slope and curvature right, so it misses only through the change in curvature, and that miss shrinks like the distance squared. CO4A 14-02-03

This fast finish is local. Far from x⋆, Newton is in the damped phase: backtracking shortens some steps, and the theory promises only a fixed drop γ in f per step. Once the gradient is small enough, it enters the quadratic phase, takes t = 1 every time and never leaves. Counting both phases bounds the total number of steps: CO4A 14-05

iterations ≤ (f(x₀) − p⋆)/γ + log₂ log₂(ε₀/ε)

The second term is tiny: six quadratic-phase steps shrink the error bound by 2⁻⁶⁴ ≈ 5·10⁻²⁰. In practice that means some damped steps, then five or six full steps to machine precision. Boyd & Vandenberghe 2004

On the widget's function every Newton step is a full step, and the log₁₀ error goes 0.82, 0.07, −0.80, −2.32, −5.25, −11.10: the last three exponents roughly double. The widget's "log of that sum" option shows both phases, with three shortened steps (t = 0.058, 0.49, 0.49) before the full steps.

5. Units don't matter: affine invariance

Draw the same hillside on two maps, one in metres and one with the north–south axis in kilometres. Newton walks the same path on the ground, in the same number of steps, on both maps. Gradient descent does not: its bowl is round in the map's units, so stretching one axis changes its direction and its step count.

The reason in plain words: Newton's bowl is fitted to the function itself. Stretch the map and the bowl stretches with it, so its bottom stays at the same spot on the ground. This is affine invariance: an invertible linear change of coordinates leaves Newton's path unchanged.

In symbols, let x = Au with A invertible, and minimise g(u) = f(Au). The chain rule gives

∇g(u) = A^T∇f(x), ∇²g(u) = A^T∇²f(x) A

So Newton's step in the new units is Δu = −A⁻¹(∇²f)⁻¹∇f, and AΔu is exactly the old Newton step. The gradient step in the new units maps back to

AΔu = −t AA^T∇f(x)

That is a different direction unless AA^T is a multiple of I. (CO4A 14-02-02 writes the product as A∇f(x)A^T, the wrong order.) The Newton decrement is unchanged too, so the stopping rule ignores units as well. CO4A 14-02-02 Boyd & Vandenberghe 2004

The numbers make the point. On Boyd and Vandenberghe's test function (their eq. 9.20, used in the widget below), start at (−1, 1) and measure y in new units, y = s·u. Newton needs 5 steps for every s. Gradient descent needs 29 at s = 1, 53 at s = 0.2, 150 at s = 0.1, 596 at s = 5 and 2,345 at s = 10.

A car's planner mixes metres, radians and metres per second. Gradient descent needs someone to pick good scalings; Newton-type methods are far less sensitive, one reason robotics solvers lean on them.

Try it: read the step counts. Pick "Log of that sum" to see Newton's damped first steps, set "Rescale y by s" to 10, or tap the plot to move the start. (BFGS, the third path, comes in section 7.)

Notice: from (−1, 1), gradient descent needs 29 steps to reach an error of 10⁻¹⁰, BFGS 7 and Newton 5. Newton's error falls 6.6, 1.2, 0.16, 4.8·10⁻³, 5.6·10⁻⁶, 7.9·10⁻¹², all with full steps, while gradient descent zig-zags across the valley.

6. The price of a Newton step

Newton buys few iterations with expensive ones. Each step needs the Hessian (n² entries) and a linear solve with it. A dense Cholesky solve costs about n³/3 multiply-adds, so a step is O(n³). Gradient descent needs O(n) per step beyond the gradient. CO4A 14-07

With 8-byte numbers, at n = 1,000 the Hessian takes 8 MB and a solve about 3.3·10⁸ operations, a fraction of a second. At n = 10⁶ it would take 8 TB. Measured in seconds, Newton's lead shrinks, and for huge n it vanishes. CMU 10-725 2016

Structure rescues Newton. If the Hessian is banded or sparse, a step can cost O(n). CO4A 14-08 Trajectory problems are the standard case: each time step couples only to its neighbours, so the Hessian is block-banded and the cost grows linearly with the horizon. That is how model predictive control (Day 17) and trajectory optimisation (Day 18) afford Newton-type steps.

Equality-constrained Newton, briefly: to minimise f(x) subject to Ax = b, minimise the bowl while staying on the plane Ax = b. The step v and a multiplier estimate w solve the KKT system from Day 10:

∇²f(x) v + A^Tw = −∇f(x)
A v = −(Ax − b)

Example: minimise ½(10x₁² + x₂²) subject to x₁ + x₂ = 1 from (0.5, 0.5). The model is exact, so one step lands on the answer (1/11, 10/11) ≈ (0.091, 0.909), with w = −10/11, the optimal multiplier ν⋆. CO4A 14-08

Where Newton runs in practice: interior-point solvers (Day 12) take Newton steps on a smoothed KKT system. IPOPT, an open-source interior-point solver for nonlinear programs, uses exact second derivatives by default (as of Oct 2026) IPOPT docs and is a common back end for robot trajectory optimisation and nonlinear MPC (Day 16). SLAM and other least-squares problems use Newton's cousins Gauss–Newton and Levenberg–Marquardt (LM card). See also the Newton card.

7. Quasi-Newton: learn the curvature from the gradients

Quasi-Newton methods skip the Hessian and measure curvature from gradients they already have. After each step, compare the gradient before and after. The change per unit of step is a curvature reading, and it costs nothing extra.

In n dimensions, call the step s = x⁺ − x and the change in gradient y = ∇f(x⁺) − ∇f(x). The new curvature estimate B⁺ must explain the change you just saw. That requirement is the secant equation: CO4A 18-01

B⁺s = y, with s^Ty > 0

The second part is the curvature condition: along the step the slope went up, which a positive definite B⁺ needs. Strictly convex f gives it for free; otherwise the Wolfe conditions in the line search enforce it.

One vector equation cannot pin down a whole matrix, so each quasi-Newton method picks, in its own sense, the B⁺ closest to the old B. Three updates matter. CMU 10-725 2016

• SR1 (symmetric rank one): B⁺ = B + uu^T/(u^Ts) with u = y − Bs. Simplest, but it can lose positive definiteness, and the denominator can be near zero. CO4A 18-02
• DFP (Davidon–Fletcher–Powell): a rank-2 update of the inverse estimate H = B⁻¹ that stays positive definite when s^Ty > 0. CO4A 18-03
• BFGS (Broyden–Fletcher–Goldfarb–Shanno): DFP with the roles of B and H swapped. It stays positive definite, recovers from bad early estimates better than DFP, and is the standard. CO4A 18-04

BFGS is usually written for the inverse estimate H, so the step is a matrix–vector product, −H∇f, with no linear solve. With ρ = 1/(y^Ts):

H⁺ = (I − ρsy^T) H (I − ρys^T) + ρss^T

Both the update and the step cost O(n²), and the start can be H = I, so no Hessian is ever computed. The idea goes back to William Davidon, whose mid-1950s computer runs kept crashing before they finished, so he needed a method that converged faster. CO4A 18

A worked example: f = ½x^TQx with Q = diag(10, 1), B = H = I and step s = (1, 1), so y = Qs = (10, 1) and s^Ty = 11 > 0. SR1 recovers the exact Hessian diag(10, 1) in one update. BFGS gives H⁺ ≈ [[0.107, −0.074], [−0.074, 1.744]]: it satisfies H⁺y = s, stays positive definite (eigenvalues 0.104 and 1.747), and its corner 0.107 is already near the true 0.1.

The payoff is superlinear convergence: the ratio between successive errors goes to zero, slower than quadratic but far faster than linear. In the widget, BFGS's first step is exactly gradient descent's, because it starts from H = I; after that it learns the curvature and finishes in 7 steps to Newton's 5. A rough rule of thumb from the slides: n quasi-Newton steps do the work of one Newton step. CO4A 14-09 CO4A 18-06

BFGS still stores an n × n matrix, 8 TB at n = 10⁶. L-BFGS (limited-memory BFGS) keeps only the last m pairs (s, y), typically m = 3 to 20, and rebuilds H∇f from them with two loops of dot products, never forming a matrix. Memory is 2mn numbers (160 MB for m = 10, n = 10⁶) and the work about 4mn multiplications per step. Chen 2025 (UW CS726) CO4A 18-07

# L-BFGS direction p = -H g from the last m pairs (s[i], y[i])
q = -g
for i in newest ... oldest:          # i = k-1 down to max(k-m, 0)
    a[i] = (s[i] @ q) / (y[i] @ s[i])
    q = q - a[i] * y[i]
p = gamma * q                        # gamma = s@y / y@y of the newest pair
for i in oldest ... newest:          # i = max(k-m, 0) up to k-1
    b = (y[i] @ p) / (y[i] @ s[i])
    p = p + (a[i] - b) * s[i]
return p                             # about 4mn multiplications

CO4A's figure of this recursion, copied from the slides, prints min(k − m, 0) where max(k − m, 0) is meant. The price of the diet is a weaker guarantee, only linear convergence in general, yet L-BFGS is the default for large smooth problems. As of Oct 2026, SciPy's minimize uses BFGS when unconstrained and L-BFGS-B with bounds, SciPy 1.18 docs scikit-learn's logistic regression defaults to lbfgs, scikit-learn 1.9 docs and IPOPT's limited-memory option keeps 6 pairs. IPOPT docs See the BFGS card.

Sources: CO4A 14-07, the CMU slides and Chen's UW notes. BFGS's O(n²) is standard arithmetic, not stated in those sources.

1. Inequalities are walls: soften them with a log barrier

Picture a robot vacuum heading for the dirtiest spot in a room. If the spot is in the middle, the walls do not matter. If it is in a corner, two walls decide the answer. You do not know in advance which walls you will end up touching.

That is the hard part of every problem with "≤" rules h_i(x) ≤ 0. A rule that holds with equality at the answer is active, and the set of active rules is the active set. With m rules there are many possible active sets, and trying them one by one gets slow. CO4A 15-1-1

One exact trick moves the rules into the cost. The indicator function I_C(x) is 0 when x is in the feasible set C and +∞ outside it. Minimizing f(x) + I_C(x) is the same problem as before. But I_C is a cliff: flat, then infinitely high. Newton's method (Day 11) learns from slope and curvature, and a cliff has neither.

The fix is a smooth wall that starts rising before you reach it. For one rule x ≤ b, use the penalty −log(b − x). The gap b − x is the slack, the room you have left. The penalty is 0 at slack 1. It is 2.30 at slack 0.1, 4.61 at slack 0.01 and 6.91 at slack 0.001. Every 10× step closer adds 2.30, and at slack 0 the penalty is infinite.

For general rules h_i(x) ≤ 0, the slack is −h_i(x), and the log barrier puts one such wall on every rule: CO4A 15-1-2

φ(x) = −Σ_i log(−h_i(x))

φ is defined only where every h_i(x) < 0. Such points are strictly feasible: they keep some room from every wall. If every h_i is convex, so is φ, and it is smooth there.

To imitate the cliff, scale the wall by 1/t. The barrier parameter t > 0 sets how sharp the wall is. With t = 2, the penalty at slack 0.01 drops from 4.61 to 2.30, yet the wall itself stays infinite. So a larger t hugs the indicator more closely.

Newton also needs the slope. For a linear rule a_i^Tx ≤ b_i, write the slack as s_i = b_i − a_i^Tx. The chain rule gives: CO4A 15-1-3

∇φ(x) = Σ_i a_i / s_i
∇²φ(x) = Σ_i a_ia_i^T / s_i²

Read each term −a_i/s_i as a push away from wall i with strength 1/slack. At slack 1 the push is 1, and at slack 0.01 it is 100. Boyd and Vandenberghe call this a force field. Boyd & Vandenberghe 2004, §11.2

The point where all the pushes cancel minimizes φ by itself. It is the analytic center, the point that is farthest from all the walls in this log sense. Worked example: the pentagon x₁ ≥ 0, x₂ ≥ 0, x₁ ≤ 3, x₂ ≤ 3, x₁ + x₂ ≤ 4. At (1, 1) the five slacks are 1, 1, 2, 2 and 2, and the pushes add up to zero. So (1, 1) is the analytic center. This pentagon is today's running example.

2. The barrier problem and the central path

Put a ball on a table with springy walls, and tilt the table toward one corner. A gentle tilt leaves the ball near the middle. A steep tilt rolls it close to the low corner, but the springs keep it off the walls. The tilt is t, and the ball's resting place is x*(t).

In symbols, multiply the cost by t and add the barrier. The barrier problem is: CO4A 15-1-2

minimize t·f(x) + φ(x)
subject to Ax = b

It has the same minimizer as f + (1/t)·φ, and no inequality rules are left, so Newton's method (Day 11) can solve it. Call its minimizer x*(t).

As t grows from 0 to ∞, x*(t) traces a smooth curve called the central path. Small t means "mostly stay away from the walls", so the path starts near the analytic center. Large t means "mostly care about cost", so the path ends at the optimum x*. CO4A 15-2

The smallest example: minimize x subject to x ≥ 0. The barrier problem is minimize t·x − log x. Setting the slope t − 1/x to zero gives x*(t) = 1/t. So t = 1 gives 1, t = 10 gives 0.1, and t = 100 gives 0.01. The true answer is 0, and the error is exactly 1/t.

Now the pentagon. Minimize c^Tx with c = (−1, −2), that is, push up and to the right, with x₂ counting twice. The optimum is the corner (1, 3), with cost p* = −7. Setting the gradient of t·c^Tx + φ(x) to zero gives the centering condition:

t·c + Σ_i a_i / s_i = 0

In words: at x*(t), the pull of the cost, −t·c, balances the total push of the walls. Newton's method gives these points (m = 5 rules; section 3 explains the last column):

The path bends. x₁ first rises to about 1.14, then returns to 1 as the top wall and the diagonal wall close in. In the widget below, the cost line through x*(t) just touches the barrier level curve there, because the two gradients balance. CO4A 15-2

Try it: drag the slider from t = 0.01 to t = 1000, or press the t buttons. The orange dots are the Newton iterates of the last move. A warm start from the previous t usually needs only a few steps.

Notice: at the default t = 1, x*(t) = (1.1384, 2.1161), reached in 6 Newton steps from the analytic center (1, 1). The primal cost −5.3705 minus the dual value −10.3705 is exactly 5 = m/t. Press t = 10 and t = 100: the point moves to (1.0062, 2.8985) and (1.0001, 2.9900), and the guaranteed gap falls to 0.5 and 0.05. The dot never touches a wall.

3. The duality gap m/t: a certificate for free

A delivery app says "at most 5 minutes away". You do not know the exact arrival time, but you have a guaranteed bound. Each central point gives the same kind of promise about cost: you are at most m/t above the best value, even though you do not know that value.

Where does it come from? Day 9 showed that any dual feasible point gives a lower bound on p*, and a central point hands you one. The barrier problem's optimality condition is t∇f(x) + ∇φ(x) + A^Tw = 0, where w is the multiplier of Ax = b. Set: CO4A 15-3-1

u_i(t) = −1 / (t·h_i(x*(t))), v(t) = w / t

Each u_i is positive, because h_i < 0. For an LP it is simply u_i = 1/(t·s_i). Divide the centering condition by t, and it becomes ∇f + Σ u_i∇h_i + A^Tv = 0. That is the stationarity condition of Day 10, so x*(t) minimizes the Lagrangian at (u, v), and (u, v) is dual feasible.

Compare the whole list with the KKT conditions. Only one line really changes: complementary slackness becomes "u_ih_i = −1/t" instead of "= 0". These are the perturbed KKT conditions. As t → ∞ they turn back into KKT.

Now the gap. The dual function at (u, v) is the Lagrangian at x*(t). In it, each product u_ih_i equals −1/t, and the term v^T(Ax − b) is zero. Adding the m products gives: CO4A 15-3-2

f(x*(t)) − g(u, v) = m / t
so f(x*(t)) − p* ≤ m / t

The duality gap, primal value minus dual value, is exactly m/t, where m is the number of inequality rules. Weak duality says g(u, v) ≤ p*, which turns the gap into a bound on the true error.

Check it on the pentagon at t = 1. The slacks at x*(1) give u = (0.8784, 0.4726, 0.5372, 1.1313, 1.3413). This LP's dual value is −b^Tu = −10.3705, and the primal value is −5.3705. The gap is exactly 5.0000 = 5/1. The true error is only 1.63, so the bound is safe but not tight.

The multipliers also tell you which walls matter. As t → ∞, u approaches (0, 0, 0, 1, 1), the true KKT multipliers u*. They are nonzero only for x₂ ≤ 3 and x₁ + x₂ ≤ 4, the two walls that meet at (1, 3). Check: A^Tu* + c = (1, 2) + (−1, −2) = 0. In the 1-D example, u(t) = 1/(t·x*(t)) = 1 for every t, which is already the true multiplier.

This also gives a stopping rule. To be within ε of optimal, stop once m/t ≤ ε, that is, once t ≥ m/ε. For the pentagon with ε = 10⁻⁶, you need t ≥ 5 × 10⁶.

4. The barrier method: raise t, then re-center with Newton

If t ≥ m/ε is all we need, why not set t = 5 × 10⁶ at once? CO4A calls this barrier method v.0, and it is a bad idea. With a huge t, the barrier problem is badly conditioned: its valley is a steep, narrow canyon pressed against the walls. Newton's method started far away needs many small damped steps to get in. CO4A 15-4

Climb stairs instead. Each riser multiplies t by a factor μ > 1, which divides the guaranteed gap by μ. Each tread is a short Newton run that starts from the previous center. That start is already close to the new center, so Newton is in its fast, quadratic zone (Day 11).

Each re-centering is a centering step. Raising t is an outer iteration, and each Newton step inside a centering is an inner iteration. The factor μ is the update factor. In pseudocode: CMU 10-725 2016

given strictly feasible x, t := t0 > 0, mu > 1, eps > 0
repeat
    # centering (inner loop): Newton's method on t*f + phi,
    # subject to Ax = b, started from the current x
    x := argmin  t*f(x) + phi(x)
    if m / t <= eps: stop          # certificate from section 3
    t := mu * t                    # outer iteration

The number of outer iterations is known in advance. After k raises the guaranteed gap is m/(μ^kt⁽⁰⁾), so you need this many raises, plus the first centering: CO4A 15-5

⌈ log(m / (t⁽⁰⁾ε)) / log μ ⌉

For the pentagon with t⁽⁰⁾ = 1 and ε = 10⁻⁶, log(5 × 10⁶) = 15.4. With μ = 2 that is ⌈22.2⌉ = 23 raises, so 24 centerings. With μ = 10 there are 8 centerings, and with μ = 50 or μ = 150 there are 5.

So how should you pick μ? A small μ gives many outer iterations, but each needs only 1–3 Newton steps, because the new center is close. A large μ gives few outer iterations, each needing more Newton steps. The totals hardly change. Boyd and Vandenberghe found the total nearly flat for μ from about 3 to 200, and suggest 10 ≤ μ ≤ 100. Boyd & Vandenberghe 2004, p. 573

On the pentagon (start (1, 1), t⁽⁰⁾ = 1, ε = 10⁻⁶), μ = 2 takes 35 Newton steps in total, μ = 10 takes 39, μ = 50 takes 30 and μ = 150 takes 27. With μ = 2, the centerings take 5, 3, 2, 2, 2, 2, 2 steps and then 1 step each. Centering need not be exact either; practical codes stop each Newton run early (CO4A's barrier method v.2). CO4A 15-6

Try it: read each staircase as the gap m/t against total Newton steps. Then switch to a random LP with m = 100 rules: the problem is 20 times larger, but the totals grow far less.

Notice: on the default pentagon, μ = 2, 10, 50 and 150 need 35, 39, 30 and 27 Newton steps (24, 8, 5 and 5 centerings), and primal-dual reaches η ≤ 10⁻⁸ in 13 steps. With m = 100, the barrier totals are 73, 58, 45 and 41, and primal-dual needs 21.

5. Phase I: finding a point strictly inside

The barrier method must start strictly inside every wall, because φ is undefined anywhere else. Often you have no such point. A driving planner may start from a guess that cuts into an obstacle margin. Is there any trajectory that obeys every rule?

Ask a different question whose start is easy: how far must every rule be loosened before my guess obeys it? Add one extra variable s and solve the phase I problem: CO4A 15-7

minimize s over (x, s)
subject to h_i(x) ≤ s, Ax = b

Any x with Ax = b, together with any s larger than max_i h_i(x), is strictly feasible for this problem. So the barrier method can start on it right away.

Then read the sign of the optimal s*. If s* < 0, every rule holds with room to spare: a strictly feasible start. If s* > 0, no point obeys all the rules, so the problem is infeasible. In practice you stop phase I as soon as s drops below 0.

Pentagon example: start from the bad guess x = (5, 5). The rule values a_i^Tx − b_i are (−5, −5, 2, 2, 6), so s⁽⁰⁾ = 7 is strictly feasible for phase I. The phase I optimum is s* = −4/3 at x = (4/3, 4/3). There the slacks of x₁ ≥ 0, x₂ ≥ 0 and x₁ + x₂ ≤ 4 all equal 4/3. Since s* < 0, a strictly feasible start exists.

A variant gives each rule its own slack: minimize 1^Ts subject to h_i(x) ≤ s_i, s ≥ 0 and Ax = b. When the problem is infeasible, the nonzero s_i point to the rules that cannot be met, for example a speed limit that clashes with an arrival deadline. That makes it a good debugging tool. CMU 10-725 2016

Modern solvers rarely run a separate phase I. Primal-dual methods (next section) can start from infeasible points. MOSEK's interior-point optimizer solves a homogeneous self-dual model, which returns either a solution or a certificate that none exists. Clarabel also uses a homogeneous embedding. MOSEK docs Clarabel docs

6. Duality revisited, and the primal-dual method

CO4A's chapter 16 looks at an LP from both sides at once. Take an LP in standard form and its dual: CO4A 16-1 CMU 10-725 2016

primal: minimize c^Tx
s.t. Ax = b, x ≥ 0
dual: maximize b^Ty
s.t. A^Ty + s = c, s ≥ 0

Here s is the dual slack. The KKT conditions are Ax = b, A^Ty + s = c, x ≥ 0, s ≥ 0 and x_is_i = 0 for every i. The central path relaxes only the last condition, to x_is_i = τ with τ = 1/t > 0. CO4A 16-2

This view separates two families of LP methods. The simplex method keeps the equations and x_is_i = 0 true at every step, and fights for x, s ≥ 0. It walks from corner to corner: on the pentagon, from (0, 0) to (0, 3) and then to (1, 3). Interior-point methods keep x, s > 0 true and drive x_is_i toward 0. They cut through the middle.

The symmetry goes further. The dual of the LP barrier problem is the barrier problem of the dual LP. Semidefinite programs (Day 5) follow the same pattern. The primal "minimize tr(CX) s.t. tr(A_iX) = b_i, X ⪰ 0" has the dual "maximize b^Ty s.t. Σ y_iA_i ⪯ C". The barrier is −log det X, and the central path is XS = τI. CO4A 16-3 CO4A 17-5

Now the method. The barrier method walks to each station on the central path before heading for the next one. The primal-dual interior-point method cuts the corner: it never re-centers fully, and aims each step at a moving target on the path. CO4A 17-2

Treat the perturbed KKT conditions as equations r(x, u, v) = 0 in the primal x and the duals u and v together. Each of the three blocks is a residual: how much one group of conditions is violated. CO4A 17-2-1

r_dual = ∇f(x) + Σ u_i∇h_i(x) + A^Tv
r_cent,i = u_ih_i(x) + τ
r_prim = Ax − b

Day 11's Newton method also solves equations: linearize r around the current point, solve one linear system for the step (Δx, Δu, Δv), and move. The primal-dual method takes one such step, lowers τ, and repeats. There is one loop instead of two.

The iterates need not be feasible: Ax = b and dual feasibility hold only in the limit. So the m/t certificate does not apply directly. Instead the method tracks the surrogate duality gap: CO4A 17-2-2

η = −h(x)^Tu = Σ_i u_is_i

It equals the true duality gap once r_dual = 0 and r_prim = 0. Each iteration sets τ = σ·η/m with σ = 1/μ, a fraction of the current average u_is_i. A careful line search keeps the iterates interior. It caps the step so that u stays positive, starts at 0.99 of that cap, halves until h(x) < 0, and halves again until the residual has shrunk enough. CO4A 17-2-3

Pentagon numbers. Start at x = (1, 1) with u_i = 1/s_i = (1, 1, 0.5, 0.5, 0.5), so η = 5. The dual residual A^Tu + c has length 2.24, so this u is not dual feasible: an infeasible start. With μ = 10, η goes 5 → 2.31 → 1.64 → 0.75 → 0.40 → 0.22 → 0.020 → 0.0022, and from then on it shrinks about 9× per step. After 13 iterations, η ≤ 10⁻⁸ and ‖r_dual‖ ≤ 10⁻⁸, with x = (1, 3) and u = (0, 0, 0, 1, 1). The barrier method needed 27–39 Newton steps for less accuracy.

7. Why a few tens of Newton steps, and who uses interior-point methods

A Newton step is expensive, because it solves a linear system. A gradient step (Day 6) costs about one matrix–vector product. Interior-point methods bet on a few expensive steps; first-order methods such as ADMM (Day 13) bet on many cheap ones.

How few is "a few"? Boyd and Vandenberghe solved 2,000 random standard-form LPs, with m from 10 to 1,000 equality constraints and n = 2m variables, using μ = 100. Cutting the duality gap by a factor of 10⁴ took on average about 21 Newton steps at the small end and 27 at the large end, with a standard deviation near 2. Their summary: the count grows very slowly with size and is almost always around a few tens. Boyd & Vandenberghe 2004, p. 575

The primal-dual method does even better. In the same chapter, it reached a gap of 10⁻⁸ in 15 to 35 iterations as the size grew, roughly like log m. On their small inequality-form LP (100 rules, 50 variables), the residual fell to zero, up to rounding, in 24 iterations. Boyd & Vandenberghe 2004, pp. 613–615

Our own quick check agrees: 5 random LPs per size, m rules, m/2 variables, and absolute targets (10⁻⁶ for barrier, 10⁻⁸ for primal-dual). The size grew 100×; the step count grew less than 2×.

The worst-case theory is far more pessimistic. With a self-concordant barrier and μ = 1 + 1/√m, the provable bound grows like √m·log(1/ε) Newton steps. For one small LP in the book, that bound is about 8,000 steps, while practice needs a few tens. Nobody runs the theory's tiny μ. Boyd & Vandenberghe 2004, p. 590

What does grow with size is the cost of each step: a dense Newton system in n variables costs about n³ operations. So real solvers exploit sparsity, and MPC solvers (Day 17) exploit the stage-by-stage structure of their Newton systems along the time horizon.

Which solvers use interior-point methods? Here are five common ones, described from their own documentation (as of Oct 2026). Day 21 compares solvers in practice.

IPOPT carries the barrier idea beyond convexity, to the nonlinear problems of Day 16. For contrast, OSQP, a popular QP solver for MPC, is not an interior-point method; it uses ADMM (Day 13, OSQP paper). The algorithm cards summarize both methods of today: barrier card and primal-dual card.

A short history. Dantzig's simplex method (1940s) is fast in practice, but Klee and Minty (1972) built LPs on which it takes about 2ⁿ steps. Khachiyan's ellipsoid method (1979) was polynomial but slow. Karmarkar's interior-point method (1984) was polynomial and practical. Modern LP solvers offer both simplex and interior-point methods. CO4A 17-3 CMU 10-725 2016

In CVXPY you never write the barrier yourself. You state the problem and pick an interior-point solver:

import cvxpy as cp, numpy as np
A = np.array([[-1, 0], [0, -1], [1, 0], [0, 1], [1, 1]])
b = np.array([0, 0, 3, 3, 4])
x = cp.Variable(2)
rules = [A @ x <= b]
prob = cp.Problem(cp.Minimize(-x[0] - 2 * x[1]), rules)
prob.solve(solver=cp.CLARABEL)    # an interior-point solver
print(prob.value, x.value)        # about -7 and [1, 3]
print(rules[0].dual_value)        # about [0, 0, 0, 1, 1]
print(prob.solver_stats.num_iters)  # interior-point iterations

1. Solve through the dual: set a price, not a rule

Two delivery robots share a charging station with 4 kWh left. Robot A would like 4 kWh and robot B would like 2, so together they want 6. The station could solve one joint problem for both of them. Instead it posts a price per kWh and lets each robot decide for itself.

On Day 9 a dual variable was the price of breaking a rule. Today the price becomes a control knob. Take the problem: minimize f(x) subject to Ax = b, where x ∈ ℝⁿ and A is a matrix. Attach a price vector u to the rule, one entry per row of A. (This chapter of CO4A calls the equality price u; Day 9 wrote v.) CO4A 20-1

L(x, u) = f(x) + u^T(Ax − b)
g(u) = min_x L(x, u)

The dual function g(u) is the lowest priced cost when x may roam with no rules. Day 9 showed it is a lower bound on the optimum. Today's new fact is about its slope. Let x(u) be the x that minimizes L(x, u). Then the gradient of g at u is the constraint residual Ax(u) − b, how far the rule is from holding. The one-line reason: L is linear in u, with slope Ax − b.

So you improve a price by stepping along the residual. This is dual (sub)gradient ascent: CO4A 20-1

x^(k) ∈ argmin_x L(x, u^(k−1))
u^(k) = u^(k−1) + t(Ax^(k) − b)

Here k counts iterations and t > 0 is the step size. Too much demand (Ax > b) raises the price, and too little lowers it. If f is strictly convex, x(u) is unique, g is differentiable, and this is plain gradient ascent on g. Otherwise Ax − b is only a subgradient, as on Day 7.

Worked example. The robots' total cost is f(x) = ½(x₁ − 4)² + ½(x₂ − 2)², and the rule is x₁ + x₂ = 4. At price u, robot A minimizes ½(x₁ − 4)² + u·x₁, so it picks x₁ = 4 − u. Robot B picks x₂ = 2 − u. Start at u = 0 with t = 0.25:

The price climbs toward u* = 1, where the robots pick x* = (3, 1) and use exactly 4 kWh. The gap to u* halves at every step. Both robots cut back by the same price, but robot A, which wanted more, still gets more.

The step needs care. With t = 1 the price jumps 0, 2, 0, 2, … and never settles. Here any 0 < t < 1 converges. In general, suppose f is strongly convex with parameter d, meaning it curves upward at least like (d/2)‖x‖₂². Then the dual is smooth, and dual ascent converges like gradient descent. The box below has the details. CO4A 20-1-1

2. Dual decomposition: broadcast a price, gather the usage

Notice what the station did not do. It never looked at the robots' cost functions. It posted one number, each robot solved its own tiny problem, and the station added up the answers. That pattern scales to thousands of pieces.

It works whenever the cost is a sum of separate pieces, f(x) = f₁(x₁) + … + f_B(x_B), and the shared rule is linear, A₁x₁ + … + A_Bx_B = b. For a fixed price the Lagrangian splits into B separate terms. So the x-step splits into B independent problems. This is dual decomposition: CO4A 20-2-1

x_i^(k) ∈ argmin f_i(x_i) + (u^(k−1))^TA_ix_i
u^(k) = u^(k−1) + t(Σ_i A_ix_i^(k) − b)

CO4A calls the two halves broadcast and gather. Broadcast sends the current price to every worker, which can be a separate machine or robot. Gather collects their usage A_ix_i, and one sum updates the price. Only prices and totals travel, and no worker sees another's cost.

Shared resources are usually limits, not exact amounts: Σ_i A_ix_i ≤ b. Then the price must stay non-negative, and the update clips it at zero: CO4A 20-2-2

u^(k) = (u^(k−1) + t(Σ_i A_ix_i^(k) − b))₊

Here (v)₊ = max(v, 0), entry by entry. Now the price story is complete: raise the price of an over-used resource, lower it when there is spare capacity, and never go below zero. Example: give the station 8 kWh instead of 4. At u = 0 the robots ask for 6. The update gives (0 + 0.25(6 − 8))₊ = 0, so the price stays at zero. A resource with room to spare is free, which is complementary slackness from Day 10.

3. The augmented Lagrangian: add a spring toward the rule

Dual ascent is fragile in two ways. You saw the first: the step size must suit the problem, or the price bounces. The second is worse. If f is not strictly convex, the x-step can fail outright. Take f(x) = x with the rule x = 1. For any price u ≠ −1, L(x, u) = x + u(x − 1) is a sloped straight line in x, with no minimum.

The fix is to add a spring. The augmented Lagrangian adds a quadratic penalty on the residual: CO4A 20-3

L_ρ(x, u) = f(x) + u^T(Ax − b) + (ρ/2)‖Ax − b‖₂²

The penalty parameter ρ > 0 sets how stiff the spring is. The new term is zero at every feasible point, so the solution does not change. But now the x-step always has a minimum. In the example, x + u(x − 1) + (ρ/2)(x − 1)² is a parabola for every u.

The method of multipliers is dual ascent on L_ρ, with the step size set to ρ: CO4A 20-3

x^(k) = argmin_x L_ρ(x, u^(k−1))
u^(k) = u^(k−1) + ρ(Ax^(k) − b)

Why a step of exactly ρ? The x-step's optimality condition reads 0 ∈ ∂f(x^(k)) + A^T(u^(k−1) + ρ(Ax^(k) − b)). The bracket is exactly the new price u^(k). So each new pair (x^(k), u^(k)) meets the stationarity condition from Day 10, and only feasibility is left to fix.

Worked example: the robots again, with ρ = 1. The x-step now gives a total use of (14 − 2u)/3. From u = 0 the total is 4.667, the residual 0.667, and the new price 0.667. Then the price goes 0.889, 0.963, 0.988. The gap to u* = 1 shrinks by a factor 1/(1 + 2ρ) = 1/3 per step, and by 1/21 for ρ = 10. Any ρ > 0 converges here, and nothing bounces. In the f(x) = x example with ρ = 1, the first step gives x = 0 and u = −1, and the second lands on x = 1.

The spring has a cost of its own. Expand (ρ/2)(x₁ + x₂ − 4)² and you find a cross term ρx₁x₂. Robot A's best choice now depends on robot B's, so the x-step no longer splits. The method of multipliers converges far more reliably, but it gives up dual decomposition. CO4A 20-3 CMU 10-725 2016 The dual ascent and method of multipliers card sums up both.

4. ADMM: take turns instead of deciding together

Can we keep the spring and still split the work? Yes, if the pieces take turns. That is the idea behind the alternating direction method of multipliers (ADMM). It works on problems with two blocks of variables, each with its own cost, tied by a linear rule: CO4A 20-4-1

minimize f(x) + g(z)
subject to Ax + Bz = c

(Here g names the second cost, not the dual function.) The lasso from Day 8 fits this shape once you make a copy. It minimizes ½‖y − Xβ‖₂² + λ‖β‖₁ for data X, y and a weight λ > 0. Call the copy α, set f(β) = ½‖y − Xβ‖₂² and g(α) = λ‖α‖₁, and require β − α = 0. Each piece alone is easy: f is least squares, and g is handled by soft-thresholding. CO4A 21-3

The augmented Lagrangian of the two-block problem is

L_ρ(x, z, u) = f(x) + g(z)
+ u^T(Ax + Bz − c) + (ρ/2)‖Ax + Bz − c‖₂²

The method of multipliers would minimize it over x and z together. ADMM minimizes over x, then over z using the new x, then updates the price: CO4A 20-4-1

x^(k) = argmin_x L_ρ(x, z^(k−1), u^(k−1))
z^(k) = argmin_z L_ρ(x^(k), z, u^(k−1))
u^(k) = u^(k−1) + ρ(Ax^(k) + Bz^(k) − c)

One pass of x and then z replaces the joint minimization. The x-step sees only f and the spring, and the z-step only g and the spring. So the work splits again, while the spring keeps the robustness.

The formulas get shorter with the scaled dual w = u/ρ. Completing the square folds the price term into the spring: CO4A 20-4-3

x ← argmin_x f(x) + (ρ/2)‖Ax + Bz − c + w‖₂²
z ← argmin_z g(z) + (ρ/2)‖Ax + Bz − c + w‖₂²
w ← w + (Ax + Bz − c)

Now w is simply the running sum of all residuals so far. It works like the integral term of a PI controller: a disagreement that lasts keeps pushing harder until it is gone.

When the rule is x = z (A = I, B = −I, c = 0), each step is a proximal operator. Recall from Day 8 that prox_f,λ(v) = argmin_x f(x) + ‖x − v‖₂²/(2λ): "move toward v, but keep f small". CO4A 21-2

x ← prox_f,1/ρ(z − w)
z ← prox_g,1/ρ(x + w)
w ← w + x − z

For the lasso this reads: ridge regression, then soft-threshold, then add the disagreement. Soft-thresholding S_t shrinks each entry toward 0 by t, and sets it to 0 if its size is at most t. CO4A 21-3

β ← (X^TX + ρI)⁻¹(X^Ty + ρ(α − w))
α ← S_λ/ρ(β + w)
w ← w + β − α

Worked example with a single number. Minimize ½(β − 3)² + |β|, whose answer is β* = S₁(3) = 2. Take ρ = 1 and start at β = α = w = 0, so the β-step is β = (3 + α − w)/2. Step 1: β = 1.5, α = S₁(1.5) = 0.5, and w = 1. Step 2: β = 1.25 and α = S₁(2.25) = 1.25, so w stays at 1. After that β = α = 1.625, 1.8125, 1.906, …, and step 8 reaches 1.988. The gap to 2 halves each step. The price u = ρw = 1 is the slope of |β| at the answer, just as KKT predicts.

The same three lines in Python, with the stopping rule from the next section:

import numpy as np
from scipy.linalg import cho_factor, cho_solve

def admm_lasso(X, y, lam, rho=1.0, tol=1e-4, max_iter=2000):
    p = X.shape[1]
    F = cho_factor(X.T @ X + rho * np.eye(p))   # factor once
    Xty, alpha, w = X.T @ y, np.zeros(p), np.zeros(p)
    for k in range(1, max_iter + 1):
        beta = cho_solve(F, Xty + rho * (alpha - w))          # ridge step
        old, v = alpha, beta + w
        alpha = np.sign(v) * np.maximum(np.abs(v) - lam / rho, 0)  # soft-threshold
        w = w + beta - alpha                                   # add the disagreement
        r = np.linalg.norm(beta - alpha)                       # primal residual
        s = rho * np.linalg.norm(alpha - old)                  # dual residual
        if r <= tol and s <= tol:
            return alpha, k
    return alpha, max_iter

5. When to stop: two residuals, two questions

ADMM's stopping rule asks two questions. Do the two copies agree yet? And have they stopped moving, so that the prices have settled? Each question has its own number: Boyd et al. 2011

r^(k) = Ax^(k) + Bz^(k) − c
s^(k) = ρA^TB(z^(k) − z^(k−1))

The primal residual r measures how badly the rule is broken; for the lasso it is β − α. The dual residual s measures how much z moved in this step, times ρ; for the lasso its size is ρ‖α^(k) − α^(k−1)‖₂. After the z-step, z's own optimality condition holds exactly, and x's is off by s. So when both are small, every KKT condition nearly holds, and f + g is near the optimal value p*. Boyd et al. 2011

Look back at the scalar example. At step 2, β = α = 1.25, so r = 0. Yet the answer is 2, not 1.25. The dual residual gives it away: s = 1 × |1.25 − 0.5| = 0.75. Stop on r alone and you would quit far too early.

Boyd et al. stop when ‖r‖₂ ≤ ε_pri and ‖s‖₂ ≤ ε_dual. Each tolerance mixes an absolute and a relative part:

ε_pri = √p·ε_abs + ε_rel·max{‖Ax‖, ‖Bz‖, ‖c‖}
ε_dual = √n·ε_abs + ε_rel·‖A^Tu‖

Here p is the length of r and n the length of x, and ε_rel = 10⁻³ or 10⁻⁴ is typical. The widget below uses a simpler rule: both norms at most 10⁻⁴. Boyd et al. 2011

Does ADMM always get there? The guarantee is broad. Suppose f and g are closed, proper and convex, and the plain Lagrangian has a saddle point. Then for every ρ > 0 the residual goes to 0, f + g goes to p*, and the price goes to an optimal u*. A and B need not have full rank. CO4A 20-4-2

Speed is another matter. Boyd et al. report that ADMM often reaches "modest accuracy" "within a few tens of iterations", but can be very slow to reach high accuracy. Boyd et al. 2011 In the widget's lasso at ρ = 0.1, the tolerances 10⁻³, 10⁻⁴ and 10⁻⁶ take 81, 130 and 229 iterations. Each extra digit costs about 50 more iterations, the pattern of a first-order method like gradient descent on Day 6. If you need many digits, an interior-point method (Day 12) is usually the better fit.

6. Choosing ρ, and letting it adapt

Think of ρ as how loudly the referee shouts "agree!". Shout too softly, and each copy chases its own cost while agreement comes slowly, so r stays large. Shout too loudly, and the copies agree at once but barely move toward the optimum, so s stays large. CO4A puts it this way: a large ρ puts too little weight on minimizing f + g, and a small ρ too little on feasibility. CO4A 21-3 CMU 10-725 2016

The widget's lasso shows both patterns. At ρ = 0.1 the first three steps give ‖r‖ = 2.568, 1.597, 0.581 and ‖s‖ = 0.000, 0.110, 0.129, so r is on top. At ρ = 10 they give ‖r‖ = 0.065, 0.002, 0.003 and ‖s‖ = 2.946, 2.637, 2.304. Now s is on top, and it is still 0.92 at step 10.

Since the two residuals pull ρ in opposite directions, a simple rule can balance them. Residual balancing raises ρ when r is much larger than s, and lowers it in the opposite case: Boyd et al. 2011

if ‖r‖ > μ‖s‖: ρ ← τρ and w ← w/τ
if ‖s‖ > μ‖r‖: ρ ← ρ/τ and w ← τw

Typical values are μ = 10 and τ = 2. Rescaling w keeps the real price u = ρw unchanged. If the x-step reuses a factorization, it must be refactored whenever ρ changes. One caution: the convergence proof assumes a fixed ρ, so it covers adaptive ρ only if ρ eventually stops changing. Boyd et al. 2011

Try it: the widget runs scaled ADMM on a seeded lasso with 30 data points and 10 features (λ = 0.2516). The answer is β* ≈ (1.825, −1.336, 0, 0, 0.925, 0, …, 0), which recovers the three true features 1, 2 and 5. Slide ρ, or tap a preset, and watch which residual stays on top. Then switch on adaptive ρ and try the ends of the slider, 0.01 and 100.

Notice that at the default ρ = 1 both residuals fall below 10⁻⁴ after 16 iterations. The lower chart is a U: 130 iterations at ρ = 0.1 and 123 at ρ = 10, about 8× more, and over 1,000 at the two ends of the slider. With adaptive ρ the purple line is nearly flat: 18, 17 and 19 iterations from starts of 0.1, 1 and 10, and 16 to 25 from any start.

7. ADMM at work: lasso, consensus, projections, and OSQP

The lasso at scale. The β-step solves with the same matrix X^TX + ρI every time, and that matrix is always invertible. So you factor it once, with a Cholesky factorization costing O(p³) for p features. After that each iteration costs only O(p²). CO4A 21-3 Boyd et al. report a lasso with 1,500 examples and 5,000 features, also with one factorization up front. With ρ = 1 it met the stopping rule in 15 iterations, and the whole solve took about 1.5 s, of which about 1 s was the first factorization. Boyd et al. 2011 The group lasso, which zeroes whole groups of features at once, swaps in group soft-thresholding, R_t(v) = (1 − t/‖v‖₂)₊v. CO4A 21-3

Consensus: many machines, one answer. Suppose the cost is a sum f₁(x) + … + f_N(x), and each f_i lives on a different machine, for example because each holds part of the data. Give machine i its own copy x_i, and require every copy to equal one shared value x̄. This is consensus ADMM: CO4A 21-5

x_i ← argmin f_i(x_i) + (ρ/2)‖x_i − x̄ + w_i‖₂²
x̄ ← average of the x_i
w_i ← w_i + x_i − x̄

The x_i-steps run in parallel, as in dual decomposition, but each has a spring, as in the method of multipliers. (Strictly, the z-step averages x_i + w_i. The w_i start at zero and always sum to zero, so this is a plain average.)

Worked example: three machines hold the numbers a = 1, 2 and 6, with f_i(x) = (x − a_i)² and ρ = 1. The shared value goes 2.0, 2.667, 2.889, 2.963, 2.988, 2.996, toward the mean 3, which is the true minimizer. The primal residual shrinks by about 2/3 per step: 2.49, 1.66, 1.11, 0.74, 0.49, 0.33. The prices w_i approach (−4, −2, 6). Each one records how hard, and in which direction, that machine's data pulls away from the shared answer.

Alternating projections. To find a point in the overlap of two convex sets C and D, minimize I_C(x) + I_D(z) subject to x = z. The indicator function I_C is 0 on C and +∞ off it, so its prox step is the projection P_C, the closest point of C: CO4A 20-4-4 CMU 10-725 2016

x ← P_C(z − w), z ← P_D(x + w)
w ← w + x − z

Without w this would be the classical back-and-forth between the two sets. Boyd et al. note that the ADMM version is like Dykstra's alternating projections, which is far more efficient than the classical method. Boyd et al. 2011

More splits. Robust PCA writes a video as a low-rank background L plus a sparse foreground S. ADMM alternates singular-value soft-thresholding for L with entrywise soft-thresholding for S. CO4A 21-4 And how you split matters. For a 2-D fused lasso on an image, splitting into row and column problems that are each solved exactly converges much faster than the generic split, at a similar cost per cycle. CO4A 21-6

OSQP and model predictive control. OSQP is a free solver for quadratic programs: minimize ½x^TPx + q^Tx subject to l ≤ Ax ≤ u, with P ⪰ 0. (Here l and u are bounds, not prices.) It runs ADMM on a split chosen so that the x-step's linear system has the same matrix at almost every iteration. OSQP factors that matrix once. After that, each iteration is cheap vector work plus clipping into the box [l, u]. It is the lasso's factor-once trick, applied to any convex QP. Stellato 2020

That suits model predictive control (Day 17). There a controller solves a QP of the same shape at every control step, and only the measured state changes. OSQP can warm-start from the last solution and reuse its factorization. In the paper's MPC benchmark with a 10-step horizon, this cut the solve time from 7 ms to 2 ms with 20 states, and from 67 ms to 22 ms with 80 states. The iteration counts barely changed, so the saving came from reusing the factorization. Stellato 2020 The OSQP paper page and the ADMM card have more.

The trade-offs are ADMM's own. On the paper's 1,400 test QPs at low accuracy (10⁻³), OSQP was the fastest solver overall and failed on none. It is still a first-order method, though. On hard, badly scaled problems GUROBI was faster, and the authors call the choice of ρ "still an open research question". The paper's defaults adapt ρ during the solve and over-relax with a factor of 1.6. The current default starting ρ is 0.1 (as of Oct 2026). Stellato 2020 OSQP docs