26 - Constrained Optimisation: Lagrange Multipliers, Duality and the KKT Conditions

Machine Learning Engineer course · Module 02 · Mathematics and statistics · Lesson 11

Objective. Solving a minimisation problem when you are not free to move wherever you like. By the end you can set up and solve a problem with equality constraints using the Lagrangian, you can state and apply the four KKT conditions, you can prove weak duality and recognise when strong duality holds, and you can derive the dual of the support vector machine from its primal, reading the support vectors straight off the multipliers.

Myth to dispel. "A Lagrange multiplier is a housekeeping number: it helps you do the algebra and then you throw it away." No: λ is the derivative of the optimal value with respect to the constraint's right-hand side, that is, the marginal price of that constraint. It is often the most interesting number in the problem, and in an SVM it is literally the list of points that decide the boundary.

Introduction

In the lesson on the gradient and gradient descent we learnt how to go downhill: you compute ∇f, take a step in the opposite direction, repeat. That method, however, assumes something that almost no machine learning problem grants you: the freedom to move anywhere.

Look at the problems you have already met, or are about to. The PCA of the lesson on eigenvalues, SVD and PCA looks for a direction of norm 1: without that constraint the projected variance has no maximum, you just make the vector longer. The support vector machine looks for the widest margin subject to classifying everything correctly: the margin on its own is unbounded, it is the correct-classification constraints that turn it into a sensible problem. L2 and L1 regularisation limits how large the weights may get. A probability distribution has to sum to 1. A budget cannot be overspent.

In all these cases the minimum is not where the gradient vanishes. It is where the gradient can no longer make you better off while staying within the rules. This lesson builds the machinery that describes that situation, and it builds it all the way down: no "it can be shown that".

Final result of this lesson. You will know where the α values of an SVM's support vectors come from, deriving them yourself from the primal problem and checking them on six points you can draw by hand on a sheet of squared paper.

Lesson roadmap.

1. Equality constraints: why at the optimum the gradient must be parallel to the constraint's gradient.
2. The Lagrangian, and two problems solved all the way through, one of which is PCA.
3. Inequality constraints: active or inactive, and the four KKT conditions that follow.
4. The dual: the function that always sits below, weak duality proved, Slater stated.
5. The SVM: the dual derived from the primal, with the support vectors as a consequence.
6. Regularisation as a constraint in disguise, and why the L1 diamond produces exact zeros.
7. How they are solved in practice: penalty, projected gradient, SLSQP, with the three results compared.
8. The SVM experiment on six points, the typical mistakes, and when hard constraints are not worth it.

The eight steps should be read as a chain in which each link uses the previous one: tangency gives the multiplier, the multiplier becomes a price, the price becomes an inequality with the KKT conditions, the KKT conditions produce the dual, the dual produces the SVM. Skipping a step means ending up with formulas that seem to have fallen from the sky.

Equality constraints: the path on the hillside

The analogy: the path on the hillside

Picture a hill and, on its slope, a marked path. You cannot walk wherever you like: you can only follow the path. The question is: what is the highest point the path reaches?

It is not the summit of the hill, barring a lucky coincidence: the summit almost certainly lies off the path. And it is not a point where the ground is flat either. At the highest point of the path the terrain keeps rising — it just rises perpendicularly to the path, in a direction you are not allowed to take.

Try drawing the level curves of the hill, the lines joining points at the same height, as on a topographic map. Walk along the path. As long as the path crosses a level curve you are changing height: in one of the two directions you are going up. The highest point is the one where the path stops crossing level curves and grazes one: it is tangent to it.

From tangency to the formula

Let us translate. The height is a function f(x) with x ∈ ℝⁿ; the path is the set of points satisfying g(x) = 0. The problem is

max f(x) subject to g(x) = 0.

Two facts from calculus, both already seen: the gradient ∇f points in the direction of steepest ascent of f, and the gradient ∇g is perpendicular to the surface g(x) = 0 (along that surface g does not change, so its directional derivative is zero, so ∇g is orthogonal to every direction that stays on the surface).

Now the argument, which is a two-line proof by contradiction. Let x* be the optimum. Suppose ∇f(x*) did have a non-zero component along the constraint, that is, suppose there were a direction d tangent to the surface with ∇f(x*)·d > 0. Then by moving a small step ε along d you stay (to first order) on the surface and you gain f(x* + εd) ≈ f(x*) + ε·∇f(x*)·d > f(x*). But then x* was not the maximum: contradiction.

So ∇f(x*) has no component tangent to the surface: it is entirely perpendicular to it. And the other direction perpendicular to the surface that we know about is ∇g(x*). Two vectors perpendicular to the same (regular, n − 1 dimensional) surface are parallel. There therefore exists a number λ such that

∇f(x*) = λ ∇g(x*) ∎

This number is the Lagrange multiplier. Note carefully: we did not postulate it, we derived it. The condition that the gradients be parallel is not a rule to memorise, it is a rewriting of "there is no admissible direction left that improves the value".

The Lagrangian, and two problems solved all the way through

The Lagrangian: a trick for going back to a free problem

The condition ∇f = λ∇g plus the constraint g(x) = 0 form a system of n + 1 equations in n + 1 unknowns (the n components of x plus λ). There is an elegant way to write it all at once. Define the Lagrangian

L(x, λ) = f(x) − λ·g(x).

Setting the gradient of L with respect to all of its variables to zero gives exactly the conditions we want: ∇_x L = ∇f − λ∇g = 0 is the parallelism, and ∂L/∂λ = −g(x) = 0 is the constraint. The trick is this: a constrained problem in n variables becomes a free problem in n + 1 variables, and we already know how to work with free problems.

Beware one temptation: L must not be minimised with respect to λ. In λ the Lagrangian is linear, so it has neither a maximum nor a minimum; the point we are after is a saddle point. That is why a naive gradient descent on L in all variables does not converge, as we shall see in the section on practical methods.

Example (a): the rectangle of fixed perimeter, and the price of the constraint

Maximise the area f(x, y) = x·y subject to g(x, y) = x + y − 10 = 0. In words: among all rectangles of half-perimeter 10, which one has the largest area?

The Lagrangian is L = x·y − λ(x + y − 10). Let us set the three partial derivatives to zero:

∂L/∂x = y − λ = 0, ∂L/∂y = x − λ = 0, ∂L/∂λ = −(x + y − 10) = 0.

The first two give x = y = λ; the third gives 2λ = 10, that is λ = 5. The optimum is the square x = y = 5, of area f* = 25. Nothing surprising so far. The interesting part is λ = 5.

What is that 5? Let us write the optimal value as a function of the right-hand side c of the constraint x + y = c. Redoing the algebra gives x = y = c/2 and f*(c) = c²/4. Then df*/dc = c/2, which for c = 10 equals 5: exactly the multiplier. The multiplier is the derivative of the optimal value with respect to the right-hand side: it tells you what one more unit of constraint is worth, at the margin. In economics it is called the shadow price; in an allocation problem it tells you which constraint is worth buying first.

def ottimo_lagrange(c):
    """Maximum of x*y on the constraint x + y = c: from the Lagrangian, x = y = lambda = c/2."""
    lam = c / 2.0
    return lam, lam, lam

x, y, lam = ottimo_lagrange(10.0)
print(f"optimum: x = {x:.4f}  y = {y:.4f}  f* = {x * y:.4f}  lambda = {lam:.4f}")

# the multiplier as a marginal price: we perturb the constraint by 0.01
passo = 0.01
xd, yd, _ = ottimo_lagrange(10.0 + passo)
salto = xd * yd - x * y
print(f"f*(10.01) - f*(10.00) = {salto:.8f}")
print(f"lambda * 0.01         = {lam * passo:.8f}")
print(f"gap                   = {abs(salto - lam * passo):.2e}")

Output obtained (Python 3.14.4, no random numbers):

optimum: x = 5.0000  y = 5.0000  f* = 25.0000  lambda = 5.0000
f*(10.01) - f*(10.00) = 0.05002500
lambda * 0.01         = 0.05000000
gap                   = 2.50e-05

What the code does. It computes the optimum in closed form for a generic c, then compares the true gain from relaxing the constraint by 0.01 with the prediction λ·0.01. The discrepancy is not numerical error: it is the second-order term. Since f*(c) = c²/4, the exact expansion is f*(c + h) − f*(c) = (c/2)·h + h²/4, and with h = 0.01 the second term equals 0.0001/4, exactly the number printed. The multiplier is a derivative, hence a linear approximation: reliable for small perturbations of the constraint, not for rewriting it from scratch.

Example (b): PCA is a constrained problem (and λ is the variance)

Second example, rather less of a toy. Given a covariance matrix Σ (symmetric, positive semi-definite), maximise the projected variance f(w) = wᵀΣw subject to g(w) = wᵀw − 1 = 0. This is, word for word, the first principal component of the lesson on eigenvalues, SVD and PCA.

Without the constraint the problem has no solution: doubling w quadruples the projected variance, and you can climb for ever. It is the unit-norm constraint that makes the question meaningful, and that is why PCA is a constrained optimisation problem, not a recipe from linear algebra.

The Lagrangian is L(w, λ) = wᵀΣw − λ(wᵀw − 1). The two gradients are ∇_w(wᵀΣw) = 2Σw (this holds because Σ is symmetric) and ∇_w(wᵀw) = 2w. Setting it to zero: 2Σw − 2λw = 0, that is

Σ w = λ w.

In words: the optimal w is an eigenvector of Σ, and the multiplier is the corresponding eigenvalue. And since on the constraint wᵀw = 1 the objective value is wᵀΣw = λ·wᵀw = λ, the multiplier is the explained variance. To maximise it you take the largest eigenvalue. PCA's eigenvalue equation is not a geometric intuition fallen from the sky: it is the first-order condition of a constrained problem.

import numpy as np

Sigma = np.array([[4.0, 1.0],
                  [1.0, 2.0]])                     # made-up covariance, symmetric

lam, V = np.linalg.eigh(Sigma)                     # eigh: eigenvalues in ascending order
w = V[:, -1] * np.sign(V[0, -1])                   # dominant eigenvector, sign fixed
print(f"eigenvalues: {lam[0]:.6f}  {lam[1]:.6f}")
print(f"w (norm 1): ({w[0]:.6f}, {w[1]:.6f})   w^T w = {w @ w:.6f}")
print(f"w^T Sigma w = {w @ Sigma @ w:.6f}   (must equal lambda_max)")
print("Sigma w == lambda w ?", np.allclose(Sigma @ w, lam[-1] * w))

# independent cross-check: brute force over all directions, a tenth of a degree
angoli = np.linspace(0.0, np.pi, 1801)
W = np.column_stack([np.cos(angoli), np.sin(angoli)])
valori = np.einsum("ij,jk,ik->i", W, Sigma, W)
print(f"brute force: angle {np.degrees(angoli[valori.argmax()]):.1f} degrees "
      f"-> w^T Sigma w = {valori.max():.6f}")

Output obtained (NumPy 2.3.5, no random numbers):

eigenvalues: 1.585786  4.414214
w (norm 1): (0.923880, 0.382683)   w^T w = 1.000000
w^T Sigma w = 4.414214   (must equal lambda_max)
Sigma w == lambda w ? True
brute force: angle 22.5 degrees -> w^T Sigma w = 4.414214

What the code does. The first half applies the result of the derivation: eigh gives the eigenvalues of Σ in ascending order and the eigenvectors as columns; we take the last one, fix its sign (an eigenvector and its negative are both valid) and check that Σw = λw really holds. The second half uses no theoretical result at all: it sweeps 1801 unit directions with np.einsum and keeps the best. The two methods agree, and the cross-check is the point: a derivation with multipliers should not be believed, it should be checked with a method that does not use it. The eigenvalues are 3 ± √2: you can get them by hand from the characteristic polynomial λ² − 6λ + 7 = 0 and compare them with the printout.

Inequality constraints and the KKT conditions

Active or inactive: there are two cases, not one

Let us move from paths to fences. An inequality constraint does not force you onto a line: it shuts you inside a region. Let us write every constraint in standard form

min f(x) subject to gᵢ(x) ≤ 0, i = 1, …, m

(an inequality h(x) ≥ c is rewritten as c − h(x) ≤ 0; this is Boyd and Vandenberghe's convention, and it is worth adopting always because it pins down the sign of the multipliers).

At the optimum x* each constraint is in one of two states, and only two.

• Inactive (gᵢ(x*) < 0): the optimum sits inside the fence, away from that wall. If you remove the constraint, nothing changes. The wall is not pushing: its multiplier is zero.
• Active (gᵢ(x*) = 0): the optimum is leaning against the wall. Were the wall not there, you would go further down. The wall pushes, and how hard it pushes is the value of the multiplier.

Everything follows from this dichotomy. If all constraints are inactive, we are at the free minimum and ∇f(x*) = 0. If some are active, those constraints behave like equality constraints and we are back to the previous case — with one crucial difference about the sign.

The four conditions, one at a time

Let us define the Lagrangian with the multipliers μᵢ of the inequality constraints:

L(x, μ) = f(x) + Σᵢ μᵢ gᵢ(x).

The Karush-Kuhn-Tucker conditions are four. We state each with the sentence that makes it memorable.

1. Stationarity: ∇f(x*) + Σᵢ μᵢ ∇gᵢ(x*) = 0. The objective's gradient is balanced by the push of the walls that touch. It is the same parallelism as before, written for several constraints at once.
2. Primal feasibility: gᵢ(x*) ≤ 0 for every i. The solution has to stay inside the fence. Trivial to say, forgotten constantly: a point that violates the constraints is not a solution, however small the value of f.
3. Dual feasibility: μᵢ ≥ 0. A wall can only push inwards, never pull outwards.
4. Complementary slackness: μᵢ·gᵢ(x*) = 0 for every i. Either the constraint touches, or its price is zero. It is the exact formalisation of the active/inactive dichotomy: a product of two numbers is zero if at least one of them is.

Why μ ≥ 0, geometrically. Take a single active constraint. Stationarity says ∇f(x*) = −μ ∇g(x*). The gradient ∇g points in the direction in which g increases, that is, towards the outside of the feasible region (where g becomes positive). If μ > 0, then ∇f points inwards: the way to make f go down would be to leave, and the wall prevents it. If instead μ < 0, ∇f would point outwards: you could make f go down by moving back into the region, staying feasible. But then x* would not be the minimum. The sign is not a convention: it is the statement that the constraint is genuinely blocking the descent.

The KKT conditions are necessary under a constraint regularity condition, and they become sufficient as well when f and the gᵢ are convex. In the convex case, then, finding a point that satisfies the four conditions means you are done: it is the global optimum, not one candidate among many. That is why convex optimisation is a discipline in its own right.

A problem solved by hand. Minimise f(x, y) = x² + y² subject to x + 2y ≥ 4. In standard form: g(x, y) = 4 − x − 2y ≤ 0, and L = x² + y² + μ(4 − x − 2y). There are two cases to examine. Case μ = 0: stationarity gives 2x = 0 and 2y = 0, that is the origin; but g(0, 0) = 4 > 0, the origin violates the constraint, case discarded. Case g = 0: stationarity gives 2x − μ = 0 and 2y − 2μ = 0, that is x = μ/2 and y = μ; substituting into the active constraint, μ/2 + 2μ = 4 gives μ = 8/5 = 1.6 ≥ 0, hence x* = 0.8, y* = 1.6, f* = 3.2. Both checks pass: this is the solution.

import numpy as np

def risolvi_kkt(a, c):
    """Minimum of ||v||^2 under a . v >= c, that is g(v) = c - a . v <= 0.

    It examines the two cases and returns (v, mu) of the case that survives."""
    a = np.asarray(a, dtype=float)
    if c <= 0.0:                       # case mu = 0: the free minimum is feasible
        return np.zeros_like(a), 0.0
    mu = 2.0 * c / (a @ a)             # case g = 0: stationarity 2v = mu*a
    return mu * a / 2.0, mu

a = np.array([1.0, 2.0])
for c in (4.0, -1.0):
    v, mu = risolvi_kkt(a, c)
    g = c - a @ v
    print(f"c = {c:5.1f} -> v* = ({v[0]:.4f}, {v[1]:.4f})  f* = {v @ v:.4f}  "
          f"mu = {mu:.4f}  g(v*) = {g:+.4f}  |mu*g| = {abs(mu * g):.2e}")

Output obtained:

c =   4.0 -> v* = (0.8000, 1.6000)  f* = 3.2000  mu = 1.6000  g(v*) = +0.0000  |mu*g| = 0.00e+00
c =  -1.0 -> v* = (0.0000, 0.0000)  f* = 0.0000  mu = 0.0000  g(v*) = -1.0000  |mu*g| = 0.00e+00

What the code does. It implements the discussion of the two cases, not a numerical search: the first if is the inactive constraint, the two lines after it the active case solved in closed form, and the choice between them is a feasibility condition, not a parameter. Complementary slackness is zero in both cases, but for different reasons: with c = 4 because g = 0, with c = −1 because μ = 0. That is the meaning of the condition: it does not say the two factors are small, it says at least one of them is exactly zero. The boundary between the two regimes is c = 0.

The dual: the function that always sits below

The dual function

So far we have used the Lagrangian as a system of equations. Now we turn it into something different, and more powerful. Fix the multipliers μ ≥ 0 and minimise the Lagrangian over all of x, with no constraints at all. The result is a number, which depends on μ:

d(μ) = infₓ L(x, μ) = infₓ [ f(x) + Σᵢ μᵢ gᵢ(x) ].

This is the dual function. The problem of maximising d(μ) for μ ≥ 0 is called the dual problem; the original problem is called the primal, and its optimal value is written p*.

One remarkable fact, and it comes free: d is concave no matter what f and the gᵢ look like, even if the primal is hideous. The reason is that d is the infimum of a family of functions that are affine in μ (for each fixed x, L(x, μ) is a straight line in μ), and the infimum of affine functions is concave. So the dual problem is always a concave maximisation, that is, convex: global maximum, no local minima to get stuck in.

Weak duality, proved in three lines

Let x̃ be any feasible point of the primal, that is, with gᵢ(x̃) ≤ 0 for every i, and let μ ≥ 0. Then:

Σᵢ μᵢ gᵢ(x̃) ≤ 0, because it is a sum of products of non-negative numbers (μᵢ) and non-positive numbers (gᵢ(x̃)).

L(x̃, μ) = f(x̃) + Σᵢ μᵢ gᵢ(x̃) ≤ f(x̃), by the previous point.

d(μ) = infₓ L(x, μ) ≤ L(x̃, μ) ≤ f(x̃), because the infimum over all x is less than or equal to the value at a specific x̃.

The chain holds for every feasible x̃, so it holds for the best of them too: d(μ) ≤ p*, for every μ ≥ 0. ∎

This is weak duality. In words: any value of the dual is a guaranteed lower bound on the primal. It is a practical result, not merely an aesthetic one: if you are solving a hard problem with a method that produces feasible solutions, every μ ≥ 0 you happen to hold gives you a certificate of the form "I shall never be able to get below d(μ)". The difference p* − d* between the primal and dual optima is called the duality gap, and it is always ≥ 0.

def duale(mu):
    """d(mu) for min x^2 + y^2 with 4 - x - 2y <= 0.

    The Lagrangian L = x^2 + y^2 + mu*(4 - x - 2y) is a paraboloid in (x, y):
    its free minimum is at x = mu/2, y = mu."""
    x, y = mu / 2.0, mu
    return x * x + y * y + mu * (4.0 - x - 2.0 * y)

p_stella = 3.2                                   # worked out by hand with the KKT conditions
for mu in (0.0, 0.5, 1.0, 1.6, 2.0, 3.0, 4.0):
    d = duale(mu)
    print(f"mu = {mu:4.1f}   d(mu) = {d:7.4f}   d(mu) <= p* ? {d <= p_stella + 1e-12}")

Output obtained:

mu =  0.0   d(mu) =  0.0000   d(mu) <= p* ? True
mu =  0.5   d(mu) =  1.6875   d(mu) <= p* ? True
mu =  1.0   d(mu) =  2.7500   d(mu) <= p* ? True
mu =  1.6   d(mu) =  3.2000   d(mu) <= p* ? True
mu =  2.0   d(mu) =  3.0000   d(mu) <= p* ? True
mu =  3.0   d(mu) =  0.7500   d(mu) <= p* ? True
mu =  4.0   d(mu) = -4.0000   d(mu) <= p* ? True

What the code does. The free minimum of the Lagrangian is worked out once and for all by hand (x = μ/2, y = μ) and substituted in: the dual function, written out in full, is d(μ) = −1.25μ² + 4μ, a concave parabola with its vertex at μ = 1.6. No value exceeds 3.2, exactly as weak duality guarantees, and the maximum touches it: the duality gap is zero, and the multiplier that maximises the dual is the same one found with the KKT conditions. Note the behaviour for large μ: the dual collapses, until it becomes negative. A lower bound stays valid even when it is useless.

Strong duality and Slater's condition: what we state and where it is proved

When the duality gap is zero we say strong duality holds: d* = p*, and solving the dual is equivalent to solving the primal. It does not always hold. The most widely used sufficient condition is Slater's condition.

Statement. If the primal problem is convex (f convex, the gᵢ convex, any equality constraints affine) and there exists a strictly feasible point, that is, an x̄ in the relative interior of the domain with gᵢ(x̄) < 0 for all the non-affine inequality constraints, then strong duality holds and the dual optimum is attained. Moreover, if p* is finite, an optimal μ* exists and the pair (x*, μ*) satisfies the KKT conditions.

We do not prove it here, and I say so rather than writing "it can be shown that": the proof goes through the separating hyperplane theorem applied to the set of achievable values {(g(x), f(x))}, and it is given in full in Boyd and Vandenberghe, Convex Optimization, §5.3.2, with the geometric reading in §5.3.1. The condition is named after Morton Slater, who introduced it in 1950. What we shall do here is check it numerically on the lesson's problems: we have already seen the gap close to zero in the output above, and we shall see it again on the SVM.

Why solving the dual can cost less. Three concrete reasons, all three of which hold for the SVM. First, the number of variables changes: the SVM primal has d + 1 unknowns (the weights plus the intercept), the dual has n, one per example — and when d ≫ n, that is, with few observations and a great many features, the dual is drastically smaller. Second, the constraints get simpler: n complicated inequality constraints become n simple positivity constraints plus a single equality, and specialised algorithms exist for those. Third, and decisively, in the dual the data appear only through inner products, and that opens the door to the kernel trick. The flip side: when n ≫ d, that is, with many examples and few features, the dual has more variables than the primal and the primal is preferable — which is exactly the choice scikit-learn's LinearSVC makes with the parameter dual="auto".

The SVM: the dual derived from the primal

The hard-margin primal

Given n examples (xᵢ, yᵢ) with yᵢ ∈ {−1, +1} and linearly separable, the support vector machine looks for the hyperplane wᵀx + b = 0 that separates them with the widest margin. Normalising the scale so that the closest points satisfy |wᵀx + b| = 1, the margin equals 1/‖w‖ and maximising it is equivalent to minimising ‖w‖². The primal problem is

min_(w, b) ½‖w‖² subject to yᵢ(wᵀxᵢ + b) ≥ 1, i = 1, …, n.

In standard form the constraints are gᵢ(w, b) = 1 − yᵢ(wᵀxᵢ + b) ≤ 0. The Lagrangian, with a multiplier αᵢ ≥ 0 for every example, is

L(w, b, α) = ½‖w‖² + Σᵢ αᵢ [1 − yᵢ(wᵀxᵢ + b)] = ½‖w‖² − Σᵢ αᵢ yᵢ(wᵀxᵢ + b) + Σᵢ αᵢ.

Setting the derivatives with respect to w and b to zero

To build the dual function we have to minimise L with respect to the primal variables. L is a convex parabola in w and linear in b, so setting the derivatives to zero is enough.

∇_w L = w − Σᵢ αᵢ yᵢ xᵢ = 0, whence the first relation:

w = Σᵢ αᵢ yᵢ xᵢ.

∂L/∂b = −Σᵢ αᵢ yᵢ = 0, whence the second:

Σᵢ αᵢ yᵢ = 0.

Read the first one: the weight vector is a linear combination of the examples, weighted by their multipliers and by the sign of the class. It is not just any vector in the space: it lives in the span of the data. The second says that the weights of the two classes balance out.

Now let us substitute. The term −Σᵢ αᵢ yᵢ b vanishes by the second relation. What remains is ½‖w‖² − wᵀ(Σᵢ αᵢ yᵢ xᵢ) + Σᵢ αᵢ, and since the sum in brackets is w, this becomes ½‖w‖² − ‖w‖² + Σᵢ αᵢ = Σᵢ αᵢ − ½‖w‖². Expanding ‖w‖² with the first relation gives the dual problem:

max_α Σᵢ αᵢ − ½ Σᵢ Σⱼ αᵢ αⱼ yᵢ yⱼ (xᵢᵀxⱼ) subject to αᵢ ≥ 0, Σᵢ αᵢ yᵢ = 0.

Look at what is left. Neither w nor b appears any more. There are n variables, one per example, two extremely simple constraints, and the data appear only inside the matrix of inner products Kᵢⱼ = xᵢᵀxⱼ. This is not a cosmetic rewriting: it is the structural result of the lesson.

Two consequences you can read straight off the result

First consequence: the support vectors. Complementary slackness says αᵢ·[1 − yᵢ(wᵀxᵢ + b)] = 0 for every i. So, for every example, either αᵢ = 0, or yᵢ(wᵀxᵢ + b) = 1, that is, the point lies exactly on the margin. Turning the sentence round: the points that lie beyond the margin have a zero multiplier and, by the relation w = Σ αᵢ yᵢ xᵢ, contribute nothing to the boundary. You could delete them from the dataset, retrain, and get the very same hyperplane. The others — those with αᵢ > 0 — are the support vectors. They are not a design choice nor a heuristic: they are an arithmetic consequence of complementary slackness.

Second consequence: the kernel trick is legitimate. In the dual the data appear only as xᵢᵀxⱼ. If you replace that inner product with a function k(xᵢ, xⱼ) that is the inner product in some higher-dimensional space (Mercer's theorem guarantees this), the rest of the problem does not change by a single line: same constraints, same structure, same solver. You get a non-linear classifier without ever computing the coordinates in the new space, which could have infinite dimension. The kernel trick of the SVM lesson works because this dual form exists: in the primal you would have nowhere to slot it in.

Soft margin, in two lines. With non-separable data you allow violations ξᵢ ≥ 0 and pay for them in the objective: min ½‖w‖² + C·Σᵢ ξᵢ subject to yᵢ(wᵀxᵢ + b) ≥ 1 − ξᵢ. Redoing the same algebra (there are now two families of multipliers, those of the margin constraints and those of ξᵢ ≥ 0, and their sum is tied to C) the dual stays identical apart from one detail: 0 ≤ αᵢ ≤ C. The regularisation parameter C becomes a ceiling on the multiplier. Direct interpretation: no single point may push the boundary beyond a certain amount, and a small C means "do not let one outlying example decide on its own".

Regularisation is a constraint in disguise

Let us go back to the lesson on overfitting and regularisation, where L2 and L1 appear as terms added to the cost function. Where do they come from? From a constraint.

The problem "fit the data but do not let the weights get large" is naturally written as a constrained problem: min_β L(β) subject to ‖β‖² ≤ t. In standard form, g(β) = ‖β‖² − t ≤ 0, and the Lagrangian is L(β) + λ(‖β‖² − t). For fixed λ the term −λt is a constant and does not move the argmin: minimising the Lagrangian is therefore equivalent to minimising

L(β) + λ‖β‖²,

which is exactly the penalised form you know. The correspondence between t and λ is complementary slackness. If the constraint is active (‖β*‖² = t) then λ > 0 and the two problems have the same argmin; if the constraint is inactive, that is, if the unregularised solution is already small enough, then λ = 0 and the penalty is pointless. They are not two different techniques: they are the same problem seen from the two sides of duality. The value of λ is not "chosen" in the abstract — it corresponds to a budget t on the size of the weights.

import numpy as np
from scipy.optimize import minimize

Xd = np.array([[2.0, 0.0],
               [0.0, 1.0]])                      # made-up design matrix
yd = np.array([4.0, 1.0])
lam = 1.0

beta_ridge = np.linalg.solve(Xd.T @ Xd + lam * np.eye(2), Xd.T @ yd)
t = float(beta_ridge @ beta_ridge)               # the budget that corresponds to this lambda
print(f"penalized,   lambda = {lam}: beta = ({beta_ridge[0]:.6f}, {beta_ridge[1]:.6f})")
print(f"                             ||beta||^2 = {t:.6f}")

perdita = lambda b: float(((Xd @ b - yd) ** 2).sum())
ris = minimize(perdita, np.zeros(2), method="SLSQP",
               constraints=[{"type": "ineq", "fun": lambda b: t - b @ b}],
               options={"ftol": 1e-12})
print(f"constrained, t = {t:.6f}:  beta = ({ris.x[0]:.6f}, {ris.x[1]:.6f})")
print(f"same argmin up to 1e-6 ? {np.abs(ris.x - beta_ridge).max() < 1e-6}")

Output obtained (NumPy 2.3.5, SciPy 1.16.3):

penalized,   lambda = 1.0: beta = (1.600000, 0.500000)
                             ||beta||^2 = 2.810000
constrained, t = 2.810000:  beta = (1.600000, 0.500000)
same argmin up to 1e-6 ? True

What the code does. It solves ridge in closed form with λ = 1, measures the squared norm of the resulting solution and uses it as the budget t of the constrained problem, which it then solves numerically with SLSQP. The two argmins coincide. The reverse direction is equally true but less convenient: given t, the corresponding λ is the multiplier of the constraint.

Ball versus diamond: why L1 produces exact zeros. The same reading explains the geometric difference between L2 and L1, which the regularisation lesson shows but does not derive. The L2 constraint ‖β‖₂² ≤ t draws a ball, smooth everywhere. The L1 constraint ‖β‖₁ ≤ t draws a diamond (an octahedron in three dimensions, a polytope in d), with 2d vertices sitting on the coordinate axes. The constrained solution is the point where the level curves of the loss, which are concentric ellipses around the unregularised solution, first touch the feasible region. On a ball the point of contact is generic: the tangent is defined everywhere and contact almost always happens at a point with all coordinates non-zero. On a diamond, by contrast, the vertices stick out, and a vertex has a whole cone of normal directions: an expanding ellipse has a non-negligible probability of touching there first. And a vertex of the L1 diamond has, by construction, all coordinates zero but one. That is why L1 zeroes coefficients out and L2 merely shrinks them: it is not a magic property of the norm, it is the shape of the corners.

Solving them in practice

Penalty: simple, and the coefficient is not a multiplier

The most immediate method is to turn the constraint into a cost: instead of min f(x) subject to g(x) ≤ 0, you minimise f(x) + ρ·max(0, g(x))² with ρ large. It works, and it is one line of code. But it has two drawbacks that need to be known.

The first: the solution is always slightly infeasible. The minimum of the penalised problem violates the constraint by an amount that goes to zero like 1/ρ, never exactly to zero for finite ρ. The second: ρ is not the multiplier, and this is a very widespread conceptual error. The implicit multiplier is μ̂ = 2ρ·s, where s is the residual violation; since s falls like 1/ρ, the product converges to μ* while ρ grows without bound. Let us check it.

def penalizzazione(rho):
    """Minimum of x^2 + y^2 + rho*max(0, 4 - x - 2y)^2, in closed form.

    Setting the gradient to zero: x = rho*s, y = 2*rho*s, with s = 4/(1 + 5*rho),
    the residual violation of the constraint."""
    s = 4.0 / (1.0 + 5.0 * rho)
    x, y = rho * s, 2.0 * rho * s
    return x, y, x * x + y * y, 2.0 * rho * s        # the last value is the estimated mu

for rho in (1e1, 1e2, 1e3, 1e4, 1e5):
    x, y, f, mu = penalizzazione(rho)
    print(f"rho={rho:>7.0f} | x={x:.6f} | y={y:.6f} | f={f:.6f} | "
          f"mu={mu:.6f} | err={abs(f - 3.2):.2e}")

Output obtained:

rho=     10 | x=0.784314 | y=1.568627 | f=3.075740 | mu=1.568627 | err=1.24e-01
rho=    100 | x=0.798403 | y=1.596806 | f=3.187238 | mu=1.596806 | err=1.28e-02
rho=   1000 | x=0.799840 | y=1.599680 | f=3.198720 | mu=1.599680 | err=1.28e-03
rho=  10000 | x=0.799984 | y=1.599968 | f=3.199872 | mu=1.599968 | err=1.28e-04
rho= 100000 | x=0.799998 | y=1.599997 | f=3.199987 | mu=1.599997 | err=1.28e-05

What the code does. It solves the penalised problem in closed form, for five values of ρ one decade apart; the last column measures the distance from the exact value 3.2 obtained with the KKT conditions. Three readings. The error is divided by ten at every decade: convergence is of order 1/ρ, which is slow. The estimated multiplier 2ρs converges to 1.6, the exact μ*, while ρ reaches 100000: if anyone tells you the penalty coefficient "is" the multiplier, the two numbers in this table differ by five orders of magnitude. And there is a hidden cost: the Hessian of the penalised problem has a condition number that grows like ρ, so a generic solver slows down and loses digits precisely when ρ becomes large enough to be useful. Trade-off. Penalty is worth it when the constraint is soft and a violation of 1e-4 is irrelevant; it is not worth it when the constraint is hard or when you need the multiplier.

Projected gradient

If the feasible region has an easily computed projection — a half-space, a ball, a box, the probability simplex — you can take an ordinary descent step and then put the point back inside: x ← Proj_C(x − η∇f(x)). The solution stays feasible at every iteration, something the penalty method never guarantees.

def gradiente_proiettato(v0, passi=200, eta=0.1):
    """Projected descent for min ||v||^2 with a . v >= c."""
    a, c = np.array([1.0, 2.0]), 4.0
    v = np.array(v0, dtype=float)
    for _ in range(passi):
        v = v - eta * 2.0 * v                     # descent step along the gradient of f
        violazione = c - a @ v
        if violazione > 0.0:                      # projection onto the half-space
            v = v + (violazione / (a @ a)) * a
    return v

v = gradiente_proiettato([3.0, 3.0])
print(f"projected gradient: v* = ({v[0]:.9f}, {v[1]:.9f})  f* = {v @ v:.9f}")
print(f"distance from the KKT optimum (0.8, 1.6) below 1e-12 ? "
      f"{np.hypot(v[0] - 0.8, v[1] - 1.6) < 1e-12}")

Output obtained:

projected gradient: v* = (0.800000000, 1.600000000)  f* = 3.200000000
distance from the KKT optimum (0.8, 1.6) below 1e-12 ? True

What the code does. It starts from a distant feasible point, (3, 3), and alternates a descent step and a projection. The projection onto a half-space aᵀv ≥ c is explicit: if the point violates it, you add ((c − aᵀv)/‖a‖²)·a, the shortest hop that brings it back onto the boundary. For the first few steps the point stays feasible and the projection does nothing; after that, at each iteration it slides along the boundary towards the optimum with a constant contraction factor. Trade-off. The method is simple, keeps feasibility and scales well to large dimensions, but it requires the projection to be computable: on a set defined by generic non-linear constraints the projection is itself a constrained problem, and you go round in circles. Moreover, like the penalty method, it does not return the multipliers.

SLSQP, and where the multipliers are

For small and medium problems the practical choice is scipy.optimize.minimize with method="SLSQP" (sequential least squares programming): it approximates the problem with a sequence of constrained quadratic subproblems and solves them. Constraints are passed as a list of dictionaries with "type" equal to "eq" or "ineq", remembering that SciPy uses the convention fun(x) ≥ 0 for inequalities: this is the opposite of the standard form g(x) ≤ 0 used in this lesson, and getting that sign wrong is mistake number one with this API.

from scipy.optimize import minimize

obiettivo = lambda v: float(v[0] ** 2 + v[1] ** 2)
vincoli = [{"type": "ineq", "fun": lambda v: v[0] + 2.0 * v[1] - 4.0}]   # SciPy: fun >= 0

ris = minimize(obiettivo, x0=[0.0, 0.0], method="SLSQP",
               constraints=vincoli, options={"ftol": 1e-12})
print("success:", ris.success, "| evaluations of f:", ris.nfev > 0)
print(f"SLSQP:  v* = ({ris.x[0]:.6f}, {ris.x[1]:.6f})   f* = {ris.fun:.6f}")
print(f"matches the KKT up to 1e-9 ? {abs(ris.fun - 3.2) < 1e-9}")
print(f"mu from stationarity 2x - mu = 0: {2.0 * ris.x[0]:.6f}")

Output obtained (SciPy 1.16.3):

success: True | evaluations of f: True
SLSQP:  v* = (0.800000, 1.600000)   f* = 3.200000
matches the KKT up to 1e-9 ? True
mu from stationarity 2x - mu = 0: 1.600000

What the code does. It sets the problem up with SciPy's convention and solves it. The last line deserves attention: minimize with SLSQP does not return the multipliers in its OptimizeResult. If you need them explicitly, the option is method="trust-constr", which exposes them in the field v. Alternatively — and this is what we do here — you recover them from the stationarity condition, which always holds at the optimum and does not depend on the solver: from 2x* − μ = 0 it follows that μ = 1.6, the same number we found by hand.

What it tells us. The section's three methods — KKT by hand, projected gradient and SLSQP — agree on the same point (0.8, 1.6) and the same value 3.2 to within 1e-9. The penalty method does not: at the last row of its table it is still far off, and to get below 1e-9 it would need ρ of the order of 10⁹, with the conditioning that comes with it. This is the practical conclusion of the section: if the constraint is hard, use a method that treats it as such.

The experiment: six points on a sheet of paper

Let us close the circle on the SVM. The data are six two-dimensional points with integer coordinates, chosen by hand: three of class −1 and three of class +1, linearly separable, with two support vectors you can see at a glance. You can draw them on squared paper in ten seconds, and that is the point: no data generator, no seed, no surprises.

Let us solve the dual exactly as it was derived: maximise Σᵢ αᵢ − ½ αᵀGα with G the matrix Gᵢⱼ = yᵢyⱼ xᵢᵀxⱼ, subject to αᵢ ≥ 0 and Σᵢ αᵢyᵢ = 0. Since solvers minimise, we pass the negative.

import numpy as np
from scipy.optimize import minimize

X = np.array([[1.0, 1.0], [0.0, 1.0], [1.0, 0.0],        # class -1
              [3.0, 3.0], [4.0, 4.0], [5.0, 3.0]])       # class +1
y = np.array([-1.0, -1.0, -1.0, 1.0, 1.0, 1.0])

Z = y[:, None] * X                                       # rows y_i x_i
G = Z @ Z.T                                              # G_ij = y_i y_j (x_i . x_j)

duale_negato = lambda al: float(0.5 * al @ G @ al - al.sum())

ris = minimize(duale_negato, np.zeros(6), method="SLSQP",
               bounds=[(0.0, None)] * 6,                 # hard margin: alpha_i >= 0
               constraints=[{"type": "eq", "fun": lambda al: float(al @ y)}],
               options={"ftol": 1e-12, "maxiter": 500})

alpha = np.where(ris.x > 1e-8, ris.x, 0.0)               # below the threshold it is zero
print("alpha:", " ".join(f"{a:.4f}" for a in alpha))
print(f"dual value d(alpha*) = {-ris.fun:.6f}")
print(f"constraint sum alpha_i y_i satisfied ? {abs(alpha @ y) < 1e-9}")

Output obtained (SciPy 1.16.3, no random numbers):

alpha: 0.2500 0.0000 0.0000 0.2500 0.0000 0.0000
dual value d(alpha*) = 0.250000
constraint sum alpha_i y_i satisfied ? True

What the code does. It builds the matrix G in one line, exploiting the fact that Gᵢⱼ = (yᵢxᵢ)·(yⱼxⱼ), then minimises the negative of the dual with the two constraints derived by hand. The 1e-8 threshold is needed because a numerical solver does not return exact zeros: in the raw ris.x vector of this run there is a 2.2e-16 alongside two genuine zeros, and the two non-zero multipliers come out as 0.249999982, not 0.25. The distinction between "zero multiplier" and "small multiplier" has to be imposed by us, and it is a choice to be declared. Two multipliers out of six come out non-zero, and they are equal to each other because the constraint Σαᵢyᵢ = 0 forces it when there is one support vector per class. Now let us reconstruct the classifier's parameters and compare them with scikit-learn.

from sklearn.svm import SVC

w = (alpha * y) @ X                                      # w = sum of alpha_i y_i x_i
supporti = np.flatnonzero(alpha > 1e-8)
b = float(np.mean(y[supporti] - X[supporti] @ w))        # from y_k (w.x_k + b) = 1
print(f"w = ({w[0]:.6f}, {w[1]:.6f})   b = {b:.6f}")
print(f"margin 1/||w|| = {1.0 / np.linalg.norm(w):.6f}   (square root of 2 = {np.sqrt(2):.6f})")
print("margins y_i (w.x_i + b):", " ".join(f"{m:.4f}" for m in y * (X @ w + b)))
print("support vector indices:", supporti.tolist())

mdl = SVC(kernel="linear", C=1e6).fit(X, y)              # huge C = hard margin
print(f"w matches SVC up to 1e-6 ? {np.abs(w - mdl.coef_[0]).max() < 1e-6}")
print(f"b matches SVC up to 1e-6 ? {abs(b - mdl.intercept_[0]) < 1e-6}")
print("scikit-learn support_:", mdl.support_.tolist())

Output obtained (scikit-learn 1.8.0):

w = (0.500000, 0.500000)   b = -2.000000
margin 1/||w|| = 1.414214   (square root of 2 = 1.414214)
margins y_i (w.x_i + b): 1.0000 1.5000 1.5000 1.0000 2.0000 2.0000
support vector indices: [0, 3]
w matches SVC up to 1e-6 ? True
b matches SVC up to 1e-6 ? True
scikit-learn support_: [0, 3]

What the code does. It applies the two relations derived by setting the Lagrangian's derivatives to zero: w as a weighted combination of the examples, and b from the margin condition on the support vectors alone (we average over all of them for numerical robustness). Then it trains a linear SVC with C = 10⁶, which is the standard way of approximating the hard margin with a soft-margin implementation.

What it tells us. Everything checks out, to 1e-6. The boundary is x₁ + x₂ = 4, the margin is √2, that is, half the distance between (1, 1) and (3, 3), which equals 2√2. The column of margins is the most instructive check: it equals exactly 1 at the two points with α > 0 and more than 1 at the other four — complementary slackness, measured. And the indices that scikit-learn calls support_ are the same ones we found by solving the dual: a library attribute and a KKT condition are the same thing. One last note of honesty about the tolerance: it is 1e-6 and not 1e-9 because the real discrepancy is of the order of 1e-7, and that discrepancy comes from our numerically solved dual, not from scikit-learn, which on these six points returns an exact 0.5 and −2.

A control experiment worth running. Delete rows 1, 2, 4 and 5 from the dataset — the four points with α = 0 — and retrain on the two support vectors alone. You get the same w and the same b. Move (3, 3) to (3.5, 3.5), by contrast, and everything changes: the margin widens and the multipliers shrink. The model depends on two points out of six, and those two points are the hardest ones.

When hard constraints are not worth it

All this machinery has a cost, and it should not always be paid. The question to ask is: is the constraint genuinely hard, or is it the approximation of a wish?

• It is the approximation of a wish — "the weights must not get too large", "the predictions ought to be smooth", "I would rather the model did not use too many features". Here no natural threshold exists: you would pick the value of t by cross-validation exactly as you would pick λ. In such cases the penalty is simpler: one more term in the cost function, no constrained solver, no multiplier to interpret, and the gradient stays the one you know.
• It is genuinely hard — a budget that cannot be overspent, a probability that must sum to exactly 1, a unit norm which, if violated, makes the object meaningless, a physical capacity limit, a legal constraint. Here the penalty method produces an answer that almost obeys the rule, and that is useless: a distribution summing to 1.0003 is not a distribution. Explicit constraints are needed.

There is also a diagnostic cost. A constrained problem can be infeasible — no point satisfies all the constraints at once — and in that case the solver fails without telling you which constraint is the odd one out. With the penalty method this does not happen: you always get an answer, which may however violate everything merrily. Neither behaviour is better in absolute terms: the first is noisy and honest, the second is silent and dangerous.

The rule to take away. Use hard constraints when a violation makes the answer meaningless, not when it merely makes it less elegant.

Typical mistakes (symptom, cause, remedy)

1. The sign of λ flipped. Symptom: the multiplier comes out negative in a problem where it should be positive, or the price interpretation points the wrong way ("relaxing the constraint makes the value worse"). Cause: having written L = f − λg with constraints in the form g ≤ 0, or L = f + μg with constraints in the form g ≥ 0. Remedy: always put everything into the standard form gᵢ(x) ≤ 0 before writing the Lagrangian, and use L = f + Σ μᵢgᵢ. For equality constraints the sign is free, for inequalities it is not.
2. Slater not checked, and the dual stops below the primal. Symptom: you solve the dual, get a value, use it as though it were the solution, and the reconstructed primal point is not optimal (or is not even feasible). Cause: strong duality assumed rather than verified — a feasible region with no strictly interior points, or a non-convex problem. Remedy: check convexity and the existence of a strictly feasible point beforehand; if you cannot, compute the gap p* − d* explicitly and declare it.
3. Penalty too weak. Symptom: the solution violates the constraint by a visible amount and nobody notices, because the solver reports success. Cause: ρ chosen once and far too small; the violation falls like 1/ρ. Remedy: increase ρ in stages (10, 100, 1000…) reusing the previous solution as the starting point, and always print the residual violation next to the objective value.
4. A multiplier read as zero when it equals 1e-9. Symptom: the number of support vectors changes from run to run, or is absurdly high. Cause: the comparison alpha != 0 on floating-point numbers. Remedy: an explicit, declared threshold (alpha > 1e-8), chosen by looking at the histogram of the values, not at random. It is the same problem as the numerical rank of a matrix.
5. A constraint written g(x) ≥ 0 in a solver that expects g(x) ≤ 0. Symptom: the solver converges quickly to a solution that plainly violates the constraint, or declares the problem infeasible when it obviously has solutions. Cause: differing conventions between libraries — SciPy wants fun(x) ≥ 0 for "ineq" constraints, the literature uses g(x) ≤ 0. Remedy: after solving, evaluate the constraints at the point you found and print them. Three lines that expose a wrong sign in a second.
6. Gradient descent on the Lagrangian in all variables. Symptom: the values oscillate and never converge, with λ diverging. Cause: the point you are looking for is a saddle, a minimum in x and a maximum in λ; descending in both directions makes no sense. Remedy: descent in x and ascent in λ (Uzawa's method, or the augmented Lagrangian), or a solver that handles the constraints explicitly.
7. Forgetting to check primal feasibility. Symptom: the solution has a surprisingly good objective value. Cause: only stationarity and complementary slackness were solved, forgetting that the point must also stay inside the fence. Remedy: there are four KKT conditions, and you check all four, in order.

The seven mistakes fall into three families. Four are about sign or convention (sign of the multiplier, constraint written the wrong way round in the solver, primal feasibility forgotten, threshold on zero), two concern assumptions that were not checked (Slater, penalty too weak) and one is a method applied where it does not converge. For the first four it is enough to print the constraints and the four KKT conditions at the point found.

Best practice

• Put everything into standard form before writing a single line. Constraints as gᵢ(x) ≤ 0, equalities as hⱼ(x) = 0: half the sign errors disappear right here.
• Check the four KKT conditions on the solution, always. Four lines of code printing stationarity, primal feasibility, dual feasibility and complementary slackness are worth more than any success message from the solver.
• Declare the threshold you use to tell a zero multiplier from a small one, and pick that number by looking at the distribution of the values.
• Use the dual as a certificate. Even when you solve the primal, any μ ≥ 0 gives you a free lower bound: printing it next to the current value tells you how much room for improvement is left.
• Test the problem on a case whose solution you know by hand before launching it on the real one: six points with integer coordinates are worth more than ten thousand randomly generated ones.
• Declare the library versions when you report numbers: the code in this lesson was run with Python 3.14.4, NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.8.0, with no random numbers.

The six rules cost very little and can all be checked with a few lines of code: the standard form and the four KKT checks eliminate most of the mistakes in the previous section, the declared threshold and the declared versions make the result reproducible by whoever rereads it. If you apply only one, apply the check of the four conditions on the solution.

Minimal glossary

• Active constraint (vincolo attivo): an inequality constraint that holds with equality at the optimum, gᵢ(x*) = 0. If it is not active it is inactive, and removing it does not change the solution.
• Lagrangian (lagrangiano): the function L(x, μ) = f(x) + Σ μᵢgᵢ(x) that absorbs the constraints into the objective, turning a constrained problem into a free one with more variables.
• Multiplier (moltiplicatore): the variable associated with a constraint; it equals the derivative of the optimal value with respect to that constraint's right-hand side, that is, its marginal price.
• Dual problem (problema duale): max_μ d(μ) with μ ≥ 0, where d(μ) = infₓ L(x, μ). It is always a concave problem, even when the primal is not.
• Duality gap: the difference p* − d* between the primal and dual optima. It is always ≥ 0 (weak duality); when it is zero, strong duality holds.
• Complementary slackness (scarto complementare): the condition μᵢ·gᵢ(x*) = 0. Either the constraint touches, or its price is zero.
• Slater's condition: a sufficient hypothesis for strong duality in a convex problem — there exists a point that strictly satisfies all non-affine inequality constraints.
• Support vector (vettore di supporto): an example with multiplier αᵢ > 0 in the SVM dual; it lies exactly on the margin and is one of the few points that determine the boundary.

The glossary has a hidden order: starting from the active constraint you arrive at the Lagrangian, the multiplier, the dual, the gap and complementary slackness, that is, the same reasoning as the lesson compressed into eight definitions. If an entry sounds opaque, go back to the one that precedes it and reread the corresponding section.

Mini-project

Build a script that handles the same problem in both languages, constrained and penalised, and compares them honestly.

1. Write a function that solves min ‖v − z‖² with v constrained to the probability simplex (vᵢ ≥ 0, Σ vᵢ = 1), using SLSQP. Try it with z = (0.6, 0.1, −0.2, 0.9) and check the four KKT conditions on the solution, printing them one by one;
2. solve the same problem with a quadratic penalty on the sum constraint plus a truncation of the negative components to zero, for ρ ∈ {10, 10³, 10⁵}, and tabulate, for each ρ, the residual violation of Σ vᵢ = 1 and the distance from the constrained solution;
3. document the limiting case: what happens when z has all negative components? Which constraint becomes active, what are the multipliers, and does the penalty method notice?

The three requests should be carried out in order: first the constrained solution with its four checks, then the comparison with the penalty as ρ grows, finally the limit case, which is where you see whether you have understood which constraints become active. The comparison between the two methods must stay honest, with the residual violation printed next to each row.

Deliverable: a script that prints the constrained solution, the four KKT checks, the table of penalties and one line of comment on the limiting case.

Challenge

• Easy. Pen and paper: solve min x² + y² subject to x + y = 2. Find x, y and λ, then write the optimal value as a function of the right-hand side c and check that λ is exactly its derivative at c = 2.
• Medium. Pen and paper: write the four KKT conditions for min x² subject to x ≥ 1 and discuss both cases (active and inactive constraint), saying for each which condition discards or confirms it. Then repeat with the constraint x ≥ −1 and compare the two multipliers.
• Hard. Simulate 300 rolls of a loaded die (choose the six true probabilities yourself and declare them) and estimate the probabilities with projected gradient on the simplex, maximising the log-likelihood. The projection onto the simplex is not obvious: implement the sort-and-threshold algorithm. Compare the estimate with the empirical frequencies, which are the closed-form solution, and measure how far apart they are.
• Expert. Derive the soft-margin dual by hand, starting from the primal with the hinge loss and the slack variables, showing where the constraint 0 ≤ αᵢ ≤ C comes from. Then add a seventh, misclassified point to the lesson's six and check numerically that as C grows from 0.01 to 100 the number of support vectors decreases. Explain why, in terms of complementary slackness on the two groups of multipliers.
• Enterprise. Take the real allocation of your working time across four or five activities. Write the problem in standard form: an objective function (the value you produce, even roughly estimated, as long as you declare it), the total-hours constraint, the contractual minima on some activities, the sustainability maxima on others. Solve it with SLSQP and produce a report listing which constraints are active at the optimum and with what multiplier: those are the constraints that are costing you something, in that order.

The five levels climb from a calculation by hand to a problem with your own data: the first two are solved with pen and paper and consolidate multipliers and KKT, the third and the fourth require code and a comparison with the closed-form solution, the last brings the method to your work. It is best to do them in order, because each provides the oracle for the next.

The coding exercises for this lesson are in a single file, esercizi.py (about 13 KB): the multiplier as a marginal price, the discussion of the two KKT cases and the SVM dual solved on the six points. It requires only Python 3.12 or later and NumPy 2.3, with no SciPy and no random numbers. Run it with python3 esercizi.py: it works even with the parts to complete left empty and tells you which ones are missing.

Recap

Mental map. Equality constraint: at the optimum ∇f = λ∇g, because otherwise an admissible improving direction would be left. Lagrangian L = f − λg (equalities) or L = f + Σμᵢgᵢ (inequalities in the form gᵢ ≤ 0). Four KKT conditions: stationarity, primal feasibility, μ ≥ 0, μᵢgᵢ = 0. Dual d(μ) = infₓ L, always concave, always ≤ p*; gap zero under Slater. SVM: w = Σαᵢyᵢxᵢ, Σαᵢyᵢ = 0, dual in the inner products alone, αᵢ > 0 only on the margin points. Regularisation: ‖β‖² ≤ t and λ‖β‖² are the same problem, and t ↔ λ is complementary slackness.

Cheat sheet.

3 concepts: the multiplier is the constraint's marginal price; complementary slackness reduces the cases to two; the dual is always a lower bound, and under Slater it is exact. 3 mistakes: getting the sign wrong by not using standard form; confusing the penalty coefficient with the multiplier; reading a multiplier of 1e-9 as zero. 3 best practices: standard form first of all; check the four conditions on the solution; test on a case whose solution you know by hand. 1 exercise: the mini-project on projection onto the simplex.

Next lesson: entropy, and the discovery that almost every cost function you will use is the same formula written in three different ways.

Questions to take further

• If the dual is always concave even when the primal is not, why not always solve the dual? (Because the gap may be large, and reconstructing a feasible primal x from an optimal μ is not automatic: Lagrangian relaxation, combinatorial optimisation.)
• What happens to the multipliers when two active constraints have parallel gradients? (Constraint qualification conditions — LICQ, MFCQ — and the cases in which the KKT conditions stop being necessary: Nocedal and Wright, ch. 12.)
• How do you write a constraint "at most k non-zero coefficients"? (It is not convex: it is the L0 norm, and this lesson's L1 diamond is its most widely used convex relaxation. Lesson 44, module 03.)

The three questions point beyond the lesson: the first on the price you pay to solve the dual, the second on the assumptions that make the KKT conditions valid, the third on the link between the L0 norm and the L1 diamond. Each has a short answer in brackets and a reference to start from for going deeper.

The same series also includes: 25 - Matrix Calculus and Automatic Differentiation: From the Jacobian to an Autograd Written From Scratch (the previous lesson), 27 - Entropy, Cross-Entropy and KL Divergence: Why Loss Functions Are All the Same Thing (the next lesson).