11

Linear Constrained Optimization

Linear programming involves solving optimization problems with linear objective functions and linear constraints. Many problems are naturally described bylinear programs, including problems from fields as diverse as transportation,communication networks, manufacturing, economics, and operations research.
Many problems that are not naturally linear can often be approximated by linear
programs. Several methods have been developed for exploiting the linear structure. Modern techniques and hardware can globally minimize problems withmillions of variables and millions of constraints.1

11.1 Problem Formulation
A linear programming problem, called a linear program,2 can be expressed inseveral forms. Each linear program consists of a linear objective function and aset of linear constraints:c⊤ x

minimizex

subject towLE x ≤ bi for i ∈ {1, 2, . . .}(i )⊤wGE x ≥ b j for j ∈ {1, 2, . . .}( j)⊤

(11.1)

wEQ x = bk for k ∈ {1, 2, . . .}(k)⊤where i, j, and k vary over finite sets of constraints. Such an optimization problem
is given in example 11.1. Transforming real problems into this mathematical form
is often nontrivial. This text focuses on the algorithms for obtaining solutions, butother texts discuss how to go about modeling real problems.3 Several interesting
conversions are given in example 11.2.

This chapter is a short introduction to linear programs andone variation of the simplex algorithm used to solve them. Several textbooks are dedicated entirely to linear programs, including R. J. Vanderbei, Linear Programming: Foundations and Extensions,
4th ed. Springer, 2014. There are avariety of packages for solving linear programs, such as Convex.jl
and JuMP.jl, both of which include interfaces to open-source
and commercial solvers.
1

A quadratic program is a generalization of a linear program,where the objective function isquadratic and the constraints arelinear. Common approaches forsolving such problems include
some of the algorithms discussed
in earlier chapters, including theinterior point method, augmented
Lagrange method, and conjugate
gradient. The simplex method, covered in this chapter, has also been
adapted for optimizing quadratic
programs. J. Nocedal and S. J.
Wright, Numerical Optimization,
2nd ed. Springer, 2006.
2

See, for example, H. P. Williams,
Model Building in Mathematical Programming, 5th ed. Wiley, 2013.
3

190 c ha p te r 11. l i n e a r c on stra i n e d op ti mi z a ti on

The following problem has a linear objective and linear constraints, making
it a linear program.

Example 11.1. An example linear
program.

2x1 − 3x2 + 7x3

minimize
x1 ,x2 ,x3

subject to 2x1 + 3x2 − 8x3 ≤5

4x1 + x2 + 3x3 ≤9x1 + x2 + 2x3 =1x1 − 5x2 − 3x3 ≥ −4

Many problems can be converted into linear programs that have the same
solution. Two examples are L1 and L∞ minimization problems:
minimizekAx − bk1

The first problem is equivalent to solving

Example 11.2. Common norm minimization problems that can be converted into linear programs.

minimizekAx − bk∞1⊤ sminimizex,ssubject to Ax − b ≤sAx − b ≥ −swith the additional variables s.
The second problem is equivalent to solving
minimizetsubject to Ax − b ≤

t1

x,t

Ax − b ≥ −t1

with the additional variable t.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

11. 1. p robl e m formu l a ti on

191

11.1.1 General Form
We can write linear programs more compactly using matrices and arrive at thegeneral form:4
minimizec⊤ x
x

ALE x ≤ bLE

subject to

AGE x ≥ bGE(11.2)

4
Here, each constraint is elementwise. For example, in writinga ≤ b,

we mean ai ≤ bi for all i.

AEQ x = bEQ

11.1.2 Standard Form
The general linear program given in equation (11.2) can be converted into standard
form where all constraints are less-than inequalities and the design variables arenonnegative
minimizec⊤ xx

subject to Ax ≤ b(11.3)

x≥0

Greater-than inequalities are inverted, and equality constraints are split in two
AGE x ≥ bGE→

AEQ x = bEQ→−AGE x ≤ −bGE
 A x≤ b
EQ
EQ
−AEQ x ≤ −bEQ

(11.4)

We must ensure that all x entries are nonnegative as well. Suppose we startwith a linear program where x is not constrained to be nonnegative:
minimizexc⊤ x

subject to Ax ≤ b(11.5)© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

192 c ha p te r 11. l i n e a r c on stra i n e d op ti mi z a ti on

We replace x with x+ − x− and constrain x+ ≥ 0 and x− ≥ 0:
" #
i x+
h
⊤
⊤
minimizec−cx−
x+ ,x−" #hi x+
subject to A −A
≤b
x−
" #
x+
≥0
x−

(11.6)The linear objective function c⊤ x forms a flat ramp. The function increases inthe direction of c, and, as a result, all contour lines are perpendicular to c andparallel to one another as shown in figure 11.1.
A single inequality constraint w⊤ x ≤ b forms a half-space, or a region on one
side of a hyperplane. The hyperplane is perpendicular to w and is defined byw⊤ x = b as shown in figure 11.2. The region w⊤ x > b is on the +w side of thehyperplane, whereas w⊤ x < b is on the −w side of the hyperplane.
w⊤
x=

c⊤
x=

co

ns
t

c

Figure 11.1. The contours of a linear objective function c⊤ x, which
increase in the direction of c.

Figure 11.2. A linear constraint.

bw

w⊤ x > bw⊤ x < bHalf-spaces are convex sets (see appendix C.3), and the intersection of convex
sets is convex, as shown in figure 11.3. Thus, the feasible set of a linear program
will always form a convex set. Convexity of the feasible set, along with convexityof the objective function, implies that if we find a local feasible minimum, it isalso a global feasible minimum.
The feasible set is a convex region enclosed by flat faces. Depending on theregion’s configuration, the solution can lie at a vertex, on an edge, or on an entire
face. If the problem is not properly constrained, the solution can be unbounded,X

Figure 11.3. The intersection of linear constraints is a convex set.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

11. 1. p robl e m formu l a ti on

193

and, if the system is over-constrained, there is no feasible solution. Several suchcases are shown in figure 11.4.
Unbounded Solution

Infinite Solutions

No Solution

−c

−c

−c

−c

x2

One Solution

x∗
x∗
x1

11.1.3x1

x1

x1
Figure 11.4. Several different linear problem forms with different
solutions.

Equality Form

Linear programs are often solved in a third form, the equality form
minimizexc⊤ x

subject to Ax = b(11.7)

x≥0
where x and c each have n components, A is an m × n matrix, and b has mcomponents. In other words, we have n nonnegative design variables and asystem of m equations defining equality constraints.
The equality form has constraints in two parts. The first, Ax = b, forces thesolution to lie in an affine subspace.5 Such a constraint is convenient because
search techniques can constrain themselves to the constrained affine subspace toremain feasible. The second part of the constraints requires x ≥ 0, which forces
the solution to lie in the positive quadrant. The feasible set is thus the nonnegative
portion of an affine subspace. Example 11.3 provides a visualization for a simple
linear program.
Any linear program in standard form can be transformed to equality form. The
constraints are converted to:Ax ≤ b→Ax + s = b, s ≥ 0

5Informally, an affine subspace is avector space that has been translated such that its origin in a higherdimensional space is not necessarily 0.

(11.8)© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

194c ha p te r 11. l i n e a r c on stra i n e d op ti mi z a ti on

Example 11.3. Feasible sets forthe equality form are hyperplanes.
sConsider the standard-form linear program:
minimizex

subject to

x≥1

x

2

When we convert this to equality form, we getx2minimizex

subject to

x−s = 1x,sx, s ≥ 0

The equality constraint requires that feasible points fall on the line x − s =1. That line is a one-dimensional affine subspace of the two-dimensional
Euclidean space.

by introducing slack variables s. These variables take up the extra slack to enforce
equality.
Starting with a linear program:
minimizexc⊤ x

subject to Ax ≤ b(11.9)x≥0

We introduce the slack variables:
minimizex,ssubject to

" #
i x
c⊤ 0⊤
s
" #
h
i x
=b
A I
s
" #
x
≥0
s
h

(11.10)© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

411. 2. si mp l e x a l g ori thm 195

Example 11.4 demonstrates converting from standard to equality form.
Example 11.4. Converting a linear
program to equality form.

Consider the linear programminimizex

5x1 + 4x2

subject to 2x1 + 3x2 ≤ 5
4x1 + x2 ≤ 11

To convert to equality form, we first introduce two slack variables:
minimizex,s

5x1 + 4x2

subject to 2x1 + 3x2 + s1 = 5
4x1 + x2 + s2 = 11
s1 , s2 ≥ 0
We then split x:5( x1+ − x1− ) + 4( x2+ − x2− )minimize
x+ ,x− ,ssubject to 2( x1+ − x1− ) + 3( x2+ − x2− ) + s1 = 5

4( x1+ − x1− ) + ( x2+ − x2− ) + s2 = 11
x1+ , x1− , x2+ , x2− , s1 , s2 ≥ 0

11.2 Simplex Algorithm
The simplex algorithm solves linear programs by moving from vertex to vertex of
the feasible set.6 The method is guaranteed to arrive at an optimal solution solong as the linear program is feasible and bounded.
The simplex algorithm operates on equality-form linear programs (Ax = b,x ≥ 0). We assume that the rows of A are linearly independent.7 We also assume
that the problem has no more equality constraints than it has design variables(m ≤ n), which ensures that the problem is not over constrained. A preprocessing
phase guarantees that A satisfies these conditions.

The simplex algorithm was originally developed in the 1940s by
George Dantzig. A history of thedevelopment can be found here:
G. B. Dantzig, “Origins of the Simplex Method,” in A History of Scientific Computing, S. G. Nash, ed.,ACM, 1990, pp. 141–151.
6

A matrix whose rows are linearly
independent is said to have fullrow rank. Linear independence isachieved by removing redundantequality constraints.
7© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

196

c ha p te r 11. l i n e a r c on stra i n e d op ti mi z a ti on

11.2.1 Vertices
Linear programs in equality form have feasible sets in the form of convex polytopes,
which are geometric objects with flat faces. These polytopes are formed by theintersection of the equality constraints with the positive quadrant. Associatedwith a polytope are vertices, which are points in the feasible set that do not lie
between any other points in the feasible set.
The feasible set consists of several different types of design points. Points on
the interior are never optimal because they can be improved by moving along−c. Points on faces can be optimal only if the face is perpendicular to c. Points
on faces not perpendicular to c can be improved by sliding along the face in the
direction of the projection of −c onto the face. Similarly, points on edges can be
optimal only if the edge is perpendicular to c, and can otherwise be improvedby sliding along the projection of −c onto the edge. Finally, vertices can also be
optimal.
The simplex algorithm produces an optimal vertex. If a linear program hasa bounded solution, then it also contains at least one vertex. Furthermore, atleast one solution must lie at a vertex. In the case where an entire edge or face is
optimal, a vertex-solution is just as good as any other.
Every vertex for a linear program in equality form can be uniquely defined byn − m components of x that equal zero. These components are actively constrainedby xi ≥ 0. Figure 11.5 visualizes the active constraints needed to identify a vertex.
The equality constraint Ax = b has a unique solution when A is square.
We have assumed that m ≤ n, so choosing m design variables and setting theremaining variables to zero effectively removes n − m columns of A, yielding an
m × m constraint matrix (see example 11.5).
For a problem with 5 design variables and 3 constraints, setting 2 variables
to zero uniquely defines a point.


a11

 a21
a31

a12
a22
a32

a13
a23
a33

a14
a24
a34

 
 x1


a15 
a11
0
  
a25   x3  =  a21 
a35  x4 
a31
0

a13
a23
a33

X
vertexFigure 11.5. In two dimensions, anyvertex will have at least two active
constraints.

Example 11.5. Setting n − m components of x to zero can uniquely
define a point.

   b1
x1
a14
   
a24   x3  = b2 
b3
x4
a34

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

11. 2. si mp l e x a l g ori thm 197

The indices into the components {1, . . . , n} of any vertex can be partitioned
into two sets, B and V , such that:• The design values associated with indices in V are zero:i ∈ V =⇒ xi = 0(11.11)• The design values associated with indices in B may or may not be zero:i ∈ B =⇒ xi ≥ 0(11.12)• B has exactly m elements and V has exactly n − m elements.

We use xB to refer to the vector consisting of the components of x that are in B
and xV to refer to the vector consisting of the components of x that are in V . Notethat xV = 0.
The vertex associated with a partition (B , V ) can be obtained using the m × m
matrix AB formed by the m columns of A selected by B :8
Ax = AB xB = b→

1xB = A−
B b(11.13)Knowing xB is sufficient to construct x; the remaining design variables are zero. Algorithm 11.1 implements this procedure, and example 11.6 demonstrates verifyingthat a given design point is a vertex.
mutable struct LinearProgram
A
b
c
end
function get_vertex(B, LP)A, b, c = LP.A, LP.b, LP.cb_inds = sort!(collect(B))
AB = A[:,b_inds]
xB = AB\bx = zeros(length(c))
x[b_inds] = xB
return x
end

If B and V identify a vertex, thenthe columns of AB must be linearly
independent because Ax = b musthave exactly one solution. Hence,
AB xB = b must have exactly one
solution. This linear independence
guarantees that AB is invertible.
8

Algorithm 11.1. A method for extracting the vertex associated witha partition B and a linear program
LP in equality form. We introduce
the special type LinearProgram forlinear programs.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

198

c ha p te r 11. l i n e a r c on stra i n e d op ti mi z a ti on

While every vertex has an associated partition (B , V ), not every partition corresponds to a vertex. A partition corresponds to a vertex only if AB is nonsingularand the design obtained by applying equation (11.13) is feasible.9 Identifying partitions that correspond to vertices is nontrivial, and we show in section 11.2.4 thatfinding such a partition involves solving a linear program! The simplex algorithm
operates in two phases—an initialization phase that identifies a vertex partition andan optimization phase that transitions between vertex partitions toward a partition
corresponding to an optimal vertex. We will discuss both of these phases later inthis section.
Consider the constraints:

1
1

0
−
1

2
1

1
2
2

 

2
1
 

3 x =  −1 ,3−1For example, B = {1, 2} for theconstraints
 

 x
 
1 2 0  1
1
x2 =
1 2 1
1
x3

9

corresponds to   
1 2 x11=
1 2 x2
1which does not produce an invertible AB and does not have a uniquesolution.
Example 11.6. Verifying that a design point is a vertex for constraints
in equality form.

x≥0

Consider the design point x = [1, 1, 0, 0]. We can verify that x is feasible
and that it has no more than three nonzero components. We can choose either
B = {1, 2, 3} or B = {1, 2, 4}. Both


1
1 1


A{1,2,3} = 0 −1 2
2
1 2

and



1

A{1,2,4} = 0
2

1
−1
1


1

3
−1

are invertible. Thus, x is a vertex of the feasible set polytope.

11.2.2 First-Order Necessary Conditions
The first-order necessary conditions (FONCs) for optimality are used to determine
when a vertex is optimal and to inform how to transition to a more favorable
vertex. We construct a Lagrangian for the equality form of the linear program:10

L(x, µ, λ) = c⊤ x − µ⊤ x − λ⊤ (Ax − b)(11.14)Note that in x ≥ 0 the polarityof the inequality must be inverted
by multiplying both sides by −1,
yielding the negative sign in frontof µ. The Lagrangian can be defined with either positive or negative λ.
10© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

11. 2. si mp l e x a l g ori thm

199with the following FONCs:
1. feasibility: Ax = b, x ≥ 0

2. dual feasibility: µ ≥ 0

3. complementary slackness: µ ⊙ x = 0

4. stationarity: A⊤ λ + µ = c

The FONCs are sufficient conditions for optimality for linear programs. Thus,
if µ and λ can be computed for a given vertex and all four FONC equations are
satisfied, then the vertex is optimal.
We can decompose the stationarity condition into B and V components: A⊤ λ + µ = c
B
BBA⊤ λ + µ = c →(11.15) A⊤ λ + µ = cVVV

We can choose µB = 0 to satisfy complementary slackness. The value of λ canbe computed from B :11
A⊤µB = cBB λ + |{z}(11.16)=0
We can use this to obtain

λ = A−⊤B cB

A⊤
V λ + µV = cV

µV = cV − A⊤
Vλ

⊤
1
µV = cV − A−cB
B AV

We use A−⊤ to refer to the transpose of the inverse of A: ⊤   −1
= A⊤
A−⊤ = A−1

11

(11.17)Knowing µV allows us to assess the optimality of the vertices. If µV contains
negative components, then dual feasibility is not satisfied and the vertex is suboptimal.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

200 c ha p te r 11. l i n e a r c on stra i n e d op ti mi z a ti on

11.2.3 Optimization Phase
The simplex algorithm maintains a partition (B , V ), which corresponds to a vertexof the feasible set polytope. The partition can be updated by swapping indicesbetween B and V . Such a swap equates to moving from one vertex along an edgeof the feasible set polytope to another vertex. If the initial partition corresponds
to a vertex and the problem is bounded, the simplex algorithm is guaranteed toconverge to an optimum.
A transition x → x′ between vertices must satisfy Ax′ = b. Starting with apartition defined by B , we choose an entering index q ∈ V that is to enter B using
one of the heuristics described near the end of this section. The new vertex x′
must satisfy:
Ax′ = AB x′B + A{q} xq′ = AB xB = Ax = b(11.18)One leaving index p ∈ B in x′B becomes zero during the transition and is replacedby the column of A corresponding to index q. This action is referred to as pivoting.
We can solve for the new design point(11.19)1′x′B = xB − A−
B A{ q } x q

A particular leaving index p ∈ B becomes active when:
1x′B p = 0 = (xB ) p − A−
A
xq′
q
{
}
B
p

(11.20)and is thus obtained by increasing xq = 0 to xq′ with:xq′ = (xB ) p

1
A−
B A{ q } p(11.21)The leaving index is obtained using the minimum ratio test, which computes foreach potential leaving index and selects the one with minimum xq′ . We then swap
p and q between B and V . The edge transition is implemented in algorithm 11.2.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

11. 2. si mp l e x a l g ori thm 201

Algorithm 11.2. A method for computing the index p and the newcoordinate value xq′ obtained byincreasing index q of the vertex
defined by the partition B in theequality-form linear program LP.

function edge_transition(LP, B, q)
A, b, c = LP.A, LP.b, LP.cn = size(A, 2)
b_inds = sort(B)n_inds = sort!(setdiff(1:n, B))
AB = A[:,b_inds]
d, xB = AB\A[:,n_inds[q]], AB\bend

p, xq′ = 0, Inf
for i in 1 : length(d)if d[i] > 0v = xB[i] / d[i]if v < xq′
p, xq′ = i, v
end
end
end
return (p, xq′)

The effect that a transition has on the objective function can be computed using
xq′ . The objective function value at the new vertex is12′′c⊤ x′ = c⊤
B xB + c q x q


−1
′
′
= c⊤
B xB − AB A{ q } x q + c q x q
′
′
⊤ −1
= c⊤
B xB − cB AB A{ q } x q + c q x q

′
′
= c⊤
B xB − c q − µ q x q + c q x q= c⊤ x + µq xq′(11.22)

12Here, we used the fact that λ =⊤A−⊤
B cB and that A{q} λ = cq − µq .

(11.23)(11.24)(11.25)(11.26)

Choosing an entering index q decreases the objective function value byc⊤ x′ − c⊤ x = µq xq′(11.27)The objective function decreases only if µq is negative. In order to progress toward optimality, we must choose an index q in V such that µq is negative. If allcomponents of µV are non-negative, we have found a global optimum.
Since there can be multiple negative entries in µV , different heuristics can beused to select an entering index:13

Modern implementations use
more sophisticated rules. For example, see J. J. Forrest and D. Goldfarb, “Steepest-Edge Simplex Algorithms for Linear Programming,”
Mathematical Programming, vol. 57,
no. 1, pp. 341–374, 1992.

13© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

202 c ha p te r 11. l i n e a r c on stra i n e d op ti mi z a ti on• Greedy heuristic, which chooses a q that maximally reduces c⊤ x.
• Dantzig’s rule, which chooses the q with the most negative entry in µ. This ruleis easy to calculate, but it does not guarantee the maximum reduction in c⊤ x.
It is also sensitive to scaling of the constraints.14• Bland’s rule, which chooses the first q with a negative entry in µ. When used on
its own, Bland’s rule tends to result in poor performance in practice. However,this rule can help us prevent cycles, which occur when we return to a vertex
we have visited before without decreasing the objective function. This rule isusually used only after no improvements have been made for several iterationsof a different rule to break out of a cycle and ensure convergence.

For constraint A⊤ x = b →αA⊤ x = αb, α > 0, we do
not change the solution but theLagrange multipliers are scaled,
λ → α−1 λ.

14

One iteration of the simplex method’s optimization phase moves a vertex
partition to a neighboring vertex based on a heuristic for the entering index. Algorithm 11.3 implements such an iteration with the greedy heuristic. Example 11.7
demonstrates using the simplex algorithm starting from a known vertex partition
to solve a linear program.

11.2.4 Initialization Phase
The optimization phase of the simplex algorithm is implemented in algorithm 11.4.
Unfortunately, algorithm 11.4 requires an initial partition that corresponds to avertex. If we do not have such a partition, we must solve an auxiliary linear program
to obtain this partition as part of an initialization phase.
The auxiliary linear program to be solved in the initialization phase includes
extra variables z ∈ R m , which we seek to zero out:15
" #
i x
h
⊤
⊤
minimize
0
1
x,z
z
" #
h
i x
(11.28)=bsubject to A Z
z
" #
x
≥0
z

The values for z represent theamount by which Ax = b is violated. By zeroing out z, we finda feasible point. If, in solving theauxiliary problem, we do not finda vertex with a zeroed-out z, thenwe can conclude that the problem is infeasible. Furthermore, it
is not always necessary to add
all m extra variables, especiallywhen slack variables are included
when transforming between standard and equality form.
15© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

11. 2. si mp l e x a l g ori thm 203

Example 11.7. Solving a linear program with the simplex algorithm.

Consider the equality-form linear program with

A=

"

1
−4

1
2

1
0

#

" #
9
b=
,
2

0
,
1




3
 −1
 
c= 
 0 
0

and the initial vertex defined by B = {3, 4}. After verifying that B defines afeasible vertex, we can begin one iteration of the simplex algorithm.
We extract xB :1
xB = A−
B b=

"

1
0

0
1

# −1 " # " #
9
9
=
2
2

We then compute λ:λ = A−⊤
B cB =

"

1
0

0
1

#−⊤ " #
0=00and µV :

1µV = cV − A−
B AV

⊤

# "
1
3
−
cB =
0
−1
"

0
1

# −1 "

1−4#⊤ " # " #
3
1  0
=
−1
0
2

We find that µV contains negative elements, so our current B is suboptimal.
We will pivot on the index of the only negative element, q = 2. An edge
1transition is run from xB in the direction −A−
B A{q} = [1, 2].
Using equation (11.19), we increase xq′ until a new constraint becomes
active. In this case, xq′ = 1 causes x4 to become zero. We update our set ofbasic indices to B = {2, 3}.
In the second iteration, we find:
#
"
#
"
" #
1
0
1
.
, µV =
, λ=
xB =
1/2
−1/2
8
The vertex is optimal because µV has no negative entries. Our algorithm
thus terminates with B = {2, 3}, for which the design point is x∗ = [0, 1, 8, 0].
© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

204

c ha p te r 11. l i n e a r c on stra i n e d op ti mi z a ti on

function step_lp!(B, LP)A, b, c = LP.A, LP.b, LP.cn = size(A, 2)
b_inds = sort!(B)n_inds = sort!(setdiff(1:n, B))
AB, AV = A[:,b_inds], A[:,n_inds]
xB = AB\bcB = c[b_inds]
λ = AB' \ cB
cV = c[n_inds]
μV = cV - AV'*λAlgorithm 11.3. A single iterationof the simplex algorithm in which
the set B is moved from one vertex to a neighbor while maximallydecreasing the objective function.
Here, step_lp! takes a partition
defined by B and a linear programLP.

q, p, xq′, Δ = 0, 0, Inf, Inf
for i in 1 : length(μV)
if μV[i] < 0pi, xi′ = edge_transition(LP, B, i)
if μV[i]*xi′ < Δ
q, p, xq′, Δ = i, pi, xi′, μV[i]*xi′
end
end
end
if q == 0
return (B, true) # optimal point foundend
if isinf(xq′)
error("unbounded")
end

end

j = findfirst(isequal(b_inds[p]), B)
B[j] = n_inds[q] # swap indices
return (B, false) # new vertex but not optimalfunction minimize_lp!(B, LP)done = false
while !done
B, done = step_lp!(B, LP)
end
return B
end

Algorithm 11.4. Minimizing a linear program given a vertex partition defined by B and a linear program LP.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

11. 2. si mp l e x a l g ori thm

where Z is a diagonal matrix whose diagonal entries are+1 if b ≥ 0iZii =−1 otherwise.

205(11.29)The auxiliary linear program is solved with a partition defined by B , which
selects only the z values. The corresponding vertex has x = 0, and each z-element
is the absolute value of the corresponding b-value: z j = |b j |. This initial vertex
can easily be shown to be feasible.
Example 11.8 demonstrates using an auxiliary linear program to obtain afeasible vertex.
Example 11.8. Using an auxiliary
linear program to obtain a feasible
vertex.

Consider the equality-form linear program:
minimize
x1 ,x2 ,x3

c1 x1 + c2 x2 + c3 x3

subject to 2x1 − 1x2 + 2x3 =

1x1 , x2 , x3 ≥0

5x1 + 1x2 − 3x3 = −2

We can identify a feasible vertex by solving:
minimize

x1 ,x2 ,x3 ,z1 ,z2

z1 + z2

subject to 2x1 − 1x2 + 2x3 + z1 =1x1 , x2 , x3 , z1 , z2 ≥0

5x1 + 1x2 − 3x3 − z2 = −2

with an initial vertex defined by B = {4, 5}.
The initial vertex has:# −1 " # " #
"
1
0
1
1
(1)
−1
xB = AB bB =
=
0 −1
−2
2
and is thus x(1) = [0, 0, 0, 1, 2]. Solving the auxiliary problem yields x∗ ≈[0.045, 1.713, 1.312, 0, 0]. Thus [0.045, 1.713, 1.312] is a feasible vertex in theoriginal problem.
© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

206

c ha p te r 11. l i n e a r c on stra i n e d op ti mi z a ti on

The partition obtained by solving the auxiliary linear program will produce a
feasible design point, Ax = b because z will have been zeroed out. If z is nonzero,then the original linear program is infeasible. If z is zero, then the resulting
partition can be used as the initial partition for the optimization phase of thesimplex algorithm. The original problem must be slightly modified to incorporate
the new z variables:
" #
i x
h
⊤
⊤
minimize
c
0
x,z
z
"
#" # " #A I xbsubject to(11.30)
=
0 I z
0
" #
x
≥0
z
The z values must be included. Despite their vector counterparts being zero, it
is possible that some indices in the components of z are included in the initial
partition B . One can inspect the initial partition and include only the specificcomponents that are needed.
The solution (x∗ , z∗ ) obtained by solving the second LP will have z∗ = 0. Thus,∗x will be a solution to the original linear problem.
An implementation for the complete simplex algorithm is given in algorithm 11.5.

11.3

Dual Certificates

Suppose we have a candidate solution and we want to verify that it is optimal.
Verifying optimality using dual certificates (algorithm 11.6) is useful in many cases,
such as debugging our linear program code.
We know from the FONCs for constrained optimization that the optimal valueof the dual problem d∗ is a lower bound of the optimal value of the primal problemp∗ . Linear programs are linear and convex, and one can show that the optimalvalue of the dual problem is also the optimal value of the primal problem, d∗ = p∗ .

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

11. 3. d u a l c e rti fi c a te sfunction minimize_lp(LP)A, b, c = LP.A, LP.b, LP.cm, n = size(A)z = ones(m)Z = Matrix(Diagonal([j ≥ 0 ? 1 : -1 for j in b]))207

Algorithm 11.5. The simplex algorithm for solving linear programsin equality form when an initial
partition is not known.

A′ = hcat(A, Z)
b′ = bc′ = vcat(zeros(n), z)
LP_init = LinearProgram(A′, b′, c′)
B = collect(1:m).+nminimize_lp!(B, LP_init)
if any(i-> i > n, B)
error("infeasible")
end

end

A′′ = [AMatrix(1.0I, m, m);zeros(m,n) Matrix(1.0I, m, m)]b′′ = vcat(b, zeros(m))
c′′ = c′
LP_opt = LinearProgram(A′′, b′′, c′′)minimize_lp!(B, LP_opt)
return get_vertex(B, LP_opt)[1:n]© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

208

c ha p te r 11. l i n e a r c on stra i n e d op ti mi z a ti on

The primal linear program can be converted to its dual form:
max min L(x, µ, λ) = max min c⊤ x − µ⊤ x − λ⊤ (Ax − b)(11.31)= max min(c − µ − A⊤ λ)⊤ x + λ⊤ b

(11.32)

x

µ≥0,λ

µ≥0,λ

x

µ≥0,λ

x

From the FONC, we know c − µ − A⊤ λ = 0, which allows us to drop the firstterm in the objective above. In addition, we know µ = c − A⊤ λ ≥ 0, which
implies A⊤ λ ≤ c. In summary, we have:Primal Form (equality)
minimizex

Dual Formc⊤ xsubject to Ax = b
x≥0

maximize
λb⊤ λsubject to A⊤ λ ≤ c

If the primal problem has n variables and m equality constraints, then the dual
problem has m variables and n constraints.16 Furthermore, the dual of the dual is
the primal problem.
Optimality can be assessed by verifying three properties. If someone claims∗(x , λ∗ ) is optimal, we can quickly verify the claim by checking whether all threeof the following conditions are satisfied:
1. x∗ is feasible in the primal problem.
2. λ∗ is feasible in the dual problem.

An alternative to the simplex algorithm, the self-dual simplex algorithm, tends to be faster in practice. It does not require that the ma1
trix AB satisfy xB = A−B b ≥ 0.
The self-dual simplex algorithm isa modification of the simplex algorithm for the dual of the linear
programming problem in standard
form.
16

3. p∗ = c⊤ x∗ = b⊤ λ∗ = d∗ .
Dual certificates are used in example 11.9 to verify the solution to a linear
program.
function dual_certificate(LP, x, λ, ϵ=1e-6)
A, b, c = LP.A, LP.b, LP.c
primal_feasible = all(x .≥ 0) && A*x ≈ b
dual_feasible = all(A'*λ .≤ c)return primal_feasible && dual_feasible &&isapprox(c⋅x, b⋅λ, atol=ϵ)end

Algorithm 11.6. A method forchecking whether a candidate solution given by design point x anddual point λ for the linear program
LP in equality form is optimal. The
parameter ϵ controls the tolerance
for the equality constraint.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

11. 3. d u a l c e rti fi c a te s 209

Consider the standard-form linear program with

 

1
1
1 −1

 

A =  −120 , b =  −2 ,5
1
2
3

Example 11.9. Verifying a solution
using dual certificates.


1 
c =  1−1


We would like to determine whether x∗ = [2, 0, 1] and λ∗ = [1, 0, 0] are an
optimal solution pair. We first verify that x∗ is feasible:
Ax∗ = [1, −2, 5] = b,x∗ ≥ 0We then verify that λ∗ is dual-feasible:A⊤ λ∗ = [1, 1, −1] ≤ cFinally, we verify that p∗ and d∗ are the same:p ∗ = c⊤ x∗ = 1 = b⊤ λ∗ = d ∗We conclude that (x∗ , λ∗ ) are optimal.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

210

11.4

c ha p te r 11. l i n e a r c on stra i n e d op ti mi z a ti on

Summary• Linear programs are problems consisting of a linear objective function andlinear constraints.
• The simplex algorithm can optimize linear programs globally in an efficientmanner.
• Dual certificates allow us to verify that a candidate primal-dual solution pairis optimal.

11.5 Exercises
Exercise 11.1. Suppose you do not know any optimization algorithm for solving
a linear program. You decide to evaluate all the vertices and determine, by inspection, which one minimizes the objective function. Give a loose upper boundon the number of possible minimizers you will examine. Furthermore, does thismethod properly handle all linear constrained optimization problems?
Exercise 11.2. If the program in example 11.1 is bounded below, argue that thesimplex method must converge.
Exercise 11.3. Suppose we want to minimize 6x1 + 5x2 subject to the constraint
3x1 − 2x2 ≥ 5. How would we translate this problem into a linear program inequality form with the same minimizer?
Exercise 11.4. Suppose your optimization algorithm has found a search directiond and you want to conduct a line search. However, you know that there is a linear
constraint w⊤ x ≥ 0. How would you modify the line search to take this constraint
into account? You can assume that your current design point is feasible.
Exercise 11.5. Reformulate the linear programminimizexc⊤ x

subject to Ax ≥ 0(11.33)

into an unconstrained optimization problem with a log barrier penalty.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com