Linear programming is the problem of finding a vector *x* that
minimizes a linear function *f ^{T}x* subject
to linear constraints:

$$\underset{x}{\mathrm{min}}{f}^{T}x$$

such that one or more of the following hold:

*A·x* ≤ *b**Aeq·x* = *beq**l* ≤ *x* ≤ *u*.

`linprog`

AlgorithmThe `linprog`

`'interior-point'`

algorithm
is very similar to the `interior-point-convex`

`quadprog`

Algorithm.
It also shares many features with the `linprog`

`'interior-point-legacy'`

algorithm.
These algorithms have the same general outline:

Presolve, meaning simplification and conversion of the problem to a standard form.

Generate an initial point. The choice of an initial point is especially important for solving interior-point algorithms efficiently, and this step can be time-consuming.

Predictor-corrector iterations to solve the KKT equations. This step is generally the most time-consuming.

The algorithm begins by attempting to simplify the problem by removing redundancies and simplifying constraints. The tasks performed during the presolve step include:

Check if any variables have equal upper and lower bounds. If so, check for feasibility, and then fix and remove the variables.

Check if any linear inequality constraint involves just one variable. If so, check for feasibility, and change the linear constraint to a bound.

Check if any linear equality constraint involves just one variable. If so, check for feasibility, and then fix and remove the variable.

Check if any linear constraint matrix has zero rows. If so, check for feasibility, and delete the rows.

Check if the bounds and linear constraints are consistent.

Check if any variables appear only as linear terms in the objective function and do not appear in any linear constraint. If so, check for feasibility and boundedness, and fix the variables at their appropriate bounds.

Change any linear inequality constraints to linear equality constraints by adding slack variables.

If algorithm detects an infeasible or unbounded problem, it halts and issues an appropriate exit message.

The algorithm might arrive at a single feasible point, which represents the solution.

If the algorithm does not detect an infeasible or unbounded problem in the presolve step, it continues, if necessary, with the other steps. At the end, the algorithm reconstructs the original problem, undoing any presolve transformations. This final step is the postsolve step.

For simplicity, if the problem is not solved in the presolve step, the algorithm shifts all finite lower bounds to zero.

To set the initial point, `x0`

, the algorithm
does the following.

Initialize

`x0`

to`ones(n,1)`

, where`n`

is the number of elements of the objective function vector*f*.Convert all bounded components to have a lower bound of 0. If component

`i`

has a finite upper bound`u(i)`

, then`x0(i) = u/2`

.For components that have only one bound, modify the component if necessary to lie strictly inside the bound.

To put

`x0`

close to the central path, take one predictor-corrector step, and then modify the resulting position and slack variables to lie well within any bounds. For details of the central path, see Nocedal and Wright [7], page 397.

Similar to the `fmincon`

interior-point algorithm,
the `interior-point-convex`

algorithm tries to find a point
where the Karush-Kuhn-Tucker
(KKT) conditions hold. To describe these equations for the
linear programming problem, consider the standard form of the linear
programming problem after preprocessing:

$$\underset{x}{\mathrm{min}}{f}^{T}x\text{subjectto}\{\begin{array}{c}\overline{A}x=\overline{b}\\ x+t=u\\ x,t\ge 0.\end{array}$$

Assume for now that all variables have at least one finite bound. By shifting and negating components, if necessary, this assumption means that all

*x*components have a lower bound of 0.$$\overline{A}$$ is the extended linear matrix that includes both linear inequalities and linear equalities. $$\overline{b}$$ is the corresponding linear equality vector. $$\overline{A}$$ also includes terms for extending the vector

*x*with slack variables*s*that turn inequality constraints to equality constraints:$$\overline{A}x=\left(\begin{array}{cc}{A}_{eq}& 0\\ A& I\end{array}\right)\left(\begin{array}{c}{x}_{0}\\ s\end{array}\right),$$

where

*x*_{0}means the original*x*vector.*t*is the vector of slacks that convert upper bounds to equalities.

The Lagrangian for this system involves the following vectors:

*y*, Lagrange multipliers associated with the linear equalities*v*, Lagrange multipliers associated with the lower bound (positivity constraint)*w*, Lagrange multipliers associated with the upper bound

The Lagrangian is

$$L={f}^{T}x-{y}^{T}\left(\overline{A}x-\overline{b}\right)-{v}^{T}x-{w}^{T}\left(u-x-t\right).$$

Therefore, the KKT conditions for this system are

$$\begin{array}{c}f-{\overline{A}}^{T}y-v+w=0\\ \overline{A}x=\overline{b}\\ x+t=u\\ {v}_{i}{x}_{i}=0\\ {w}_{i}{t}_{i}=0\\ (x,v,w,t)\ge 0.\end{array}$$

The algorithm first predicts a step from the Newton-Raphson
formula, and then computes a corrector step. The corrector attempts
to reduce the residual in the nonlinear complementarity equations *s _{i}z_{i}* = 0. The Newton-Raphson
step is

$$\left(\begin{array}{ccccc}0& -{\overline{A}}^{T}& 0& -I& I\\ \overline{A}& 0& 0& 0& 0\\ -I& 0& -I& 0& 0\\ V& 0& 0& X& 0\\ 0& 0& W& 0& T\end{array}\right)\left(\begin{array}{c}\Delta x\\ \Delta y\\ \Delta t\\ \begin{array}{l}\Delta v\\ \Delta w\end{array}\end{array}\right)=-\left(\begin{array}{c}f-{\overline{A}}^{T}y-v+w\\ \overline{A}x-\overline{b}\\ u-x-t\\ \begin{array}{l}VX\\ WT\end{array}\end{array}\right)=-\left(\begin{array}{c}{r}_{d}\\ {r}_{p}\\ {r}_{ub}\\ \begin{array}{l}{r}_{vx}\\ {r}_{wt}\end{array}\end{array}\right),$$ | (8-1) |

where *X*, *V*, *W*,
and *T* are diagonal matrices corresponding to the
vectors *x*, *v*, *w*,
and *t* respectively. The residual vectors on the
far right side of the equation are:

*r*, the dual residual_{d}*r*, the primal residual_{p}*r*, the upper bound residual_{ub}*r*, the lower bound complementarity residual_{vx}*r*, the upper bound complementarity residual_{wt}

To solve Equation 8-1, first convert it to the symmetric matrix form

$$\left(\begin{array}{cc}-D& {\overline{A}}^{T}\\ \overline{A}& 0\end{array}\right)\left(\begin{array}{c}\Delta x\\ \Delta y\end{array}\right)=\left(\begin{array}{c}-R\\ {r}_{p}\end{array}\right),$$ | (8-2) |

where

$$\begin{array}{l}D={X}^{-1}V+{T}^{-1}W\\ R=-{r}_{d}-{X}^{-1}{r}_{vx}+{T}^{-1}{r}_{wt}+{T}^{-1}W{r}_{ub}.\end{array}$$

All the matrix inverses in the definitions of *D* and *R* are
simple to compute because the matrices are diagonal.

To derive Equation 8-2 from
Equation 8-1, notice
that the second row of Equation 8-2 is the same as the second matrix
row of Equation 8-1.
The first row of Equation 8-2 comes
from solving the last two rows of Equation 8-1 for Δ*v* and
Δ*w*, and then solving for Δ*t*.

Equation 8-2 is
symmetric, but it is not positive definite because of the –*D* term.
Therefore, you cannot solve it using a Cholesky factorization. A few
more steps lead to a different equation that is positive definite,
and hence can be solved efficiently by Cholesky factorization.

The second set of rows of Equation 8-2 is

$$\overline{A}\Delta x={r}_{p}$$

and the first set of rows is

$$-D\Delta x+{\overline{A}}^{T}\Delta y=-R.$$

Substituting

$$\Delta x={D}^{-1}{\overline{A}}^{T}\Delta y+{D}^{-1}R$$

gives

$$\overline{A}{D}^{-1}{\overline{A}}^{T}\Delta y=-\overline{A}{D}^{-1}R+{r}_{p}.$$ | (8-3) |

Usually, the most efficient way to find the Newton step is to
solve Equation 8-3 for
Δ*y* using Cholesky factorization. Cholesky
factorization is possible because the matrix multiplying Δ*y* is
obviously symmetric and, in the absence of degeneracies, is positive
definite. Afterwards, to find the Newton step, back substitute to
find Δ*x*, Δ*t*, Δ*v*,
and Δ*w*. However, when $$\overline{A}$$ has a dense column, it can be
more efficient to solve Equation 8-2 instead. The `linprog`

interior-point
algorithm chooses the solution algorithm based on the density of columns.

For more algorithm details, see Mehrotra [6].

After calculating the corrected Newton step, the algorithm performs more calculations to get both a longer current step, and to prepare for better subsequent steps. These multiple correction calculations can improve both performance and robustness. For details, see Gondzio [4].

The predictor-corrector algorithm iterates until it reaches a point that is feasible (satisfies the constraints to within tolerances) and where the relative step sizes are small. Specifically, define

$$\rho =\mathrm{max}\left(1,\Vert \overline{A}\Vert ,\Vert f\Vert ,\Vert \overline{b}\Vert \right).$$

The algorithm stops when all of these conditions are satisfied:

$$\begin{array}{c}{\Vert {r}_{p}\Vert}_{\infty}\le \rho \text{TolCon}\\ {\Vert {r}_{d}\Vert}_{\infty}\le \rho \text{TolFun}\\ {r}_{c}\le \text{TolFun,}\end{array}$$

where

$${r}_{c}=\underset{i}{\mathrm{max}}\left(\mathrm{min}\left(\left|{x}_{i}{v}_{i}\right|,\left|{x}_{i}\right|,\left|{v}_{i}\right|\right),\mathrm{min}\left(\left|{t}_{i}{w}_{i}\right|,\left|{t}_{i}\right|,\left|{w}_{i}\right|\right)\right).$$

*r _{c}* essentially measures
the size of the complementarity residuals

The default interior-point-legacy method is based on LIPSOL ([52]), which is a variant of Mehrotra's predictor-corrector algorithm ([47]), a primal-dual interior-point method.

The algorithm begins by applying a series of preprocessing steps (see Preprocessing). After preprocessing, the problem has the form

$$\underset{x}{\mathrm{min}}{f}^{T}x\text{suchthat}\{\begin{array}{c}A\cdot x=b\\ 0\le x\le u.\end{array}$$ | (8-4) |

The upper bounds constraints are implicitly included in the
constraint matrix *A*. With the addition of primal
slack variables *s*, Equation 8-4 becomes

$$\underset{x}{\mathrm{min}}{f}^{T}x\text{suchthat}\{\begin{array}{c}A\cdot x=b\\ x+s=u\\ x\ge 0,\text{}s\ge 0.\end{array}$$ | (8-5) |

which is referred to as the *primal* problem: *x* consists
of the primal variables and *s* consists of the primal
slack variables. The *dual* problem is

$$\mathrm{max}{b}^{T}y-{u}^{T}w\text{suchthat}\{\begin{array}{c}{A}^{T}\cdot y-w+z=f\\ z\ge 0,\text{}w\ge 0,\end{array}$$ | (8-6) |

where *y* and *w* consist
of the dual variables and* z* consists of the dual
slacks. The optimality conditions for this linear program, i.e., the
primal Equation 8-5 and the dual Equation 8-6, are

$$\begin{array}{l}F(x,y,z,s,w)=\left(\begin{array}{c}A\cdot x-b\\ x+s-u\\ {A}^{T}\cdot y-w+z-f\\ {x}_{i}{z}_{i}\\ {s}_{i}{w}_{i}\end{array}\right)=0,\\ \text{}x\ge 0,\text{}z\ge 0,\text{}s\ge 0,\text{}w\ge 0,\end{array}$$ | (8-7) |

where *x _{i}z_{i}* and

The quadratic equations *x _{i}z_{i}* = 0 and

*x ^{T}z* +

is the *duality gap*, which measures the residual of
the complementarity portion of *F* when (*x,z,s,w*) ≥ 0.

The algorithm is a *primal-dual algorithm*, meaning
that both the primal and the dual programs are solved simultaneously.
It can be considered a Newton-like method, applied to the linear-quadratic
system *F*(*x,y,z,s,w*) = 0 in Equation 8-7, while at the same time keeping the iterates *x*, *z*, *w*,
and *s* positive, thus the name interior-point method.
(The iterates are in the strictly interior region represented by the
inequality constraints in Equation 8-5.)

The algorithm is a variant of the predictor-corrector algorithm proposed
by Mehrotra. Consider an iterate *v* = [*x;y;z;s;w*],
where [*x;z;s;w*] > 0 First compute the so-called *prediction* direction

$$\Delta {v}_{p}=-{\left({F}^{T}(v)\right)}^{-1}F(v),$$

which is the Newton direction; then the so-called *corrector* direction

$$\Delta {v}_{c}=-{\left({F}^{T}(v)\right)}^{-1}F\left(v+\Delta {v}_{p}\right)-\mu \widehat{e},$$

where *μ* > 0 is called the *centering* parameter
and must be chosen carefully. $$\widehat{e}$$ is a zero-one vector with the ones corresponding to the
quadratic equations in *F*(*v*),
i.e., the perturbations are only applied to the complementarity conditions,
which are all quadratic, but not to the feasibility conditions, which
are all linear. The two directions are combined with a step length
parameter *α* > 0 and update *v* to
obtain the new iterate *v*^{+}:

$${v}^{+}=v+\alpha \left(\Delta {v}_{p}+\Delta {v}_{c}\right),$$

where the step length parameter *α* is
chosen so that

*v ^{+}* = [

satisfies

[*x ^{+};z^{+};s^{+};w^{+}*]
> 0.

In solving for the preceding predictor/corrector directions,
the algorithm computes a (sparse) direct factorization on a modification
of the Cholesky factors of *A·A ^{T}*.
If

`ldl`

function
reference page.) The algorithm then loops until the iterates converge. The main stopping criteria is a standard one:

$$\frac{\Vert {r}_{b}\Vert}{\mathrm{max}\left(1,\Vert b\Vert \right)}+\frac{\Vert {r}_{f}\Vert}{\mathrm{max}\left(1,\Vert f\Vert \right)}+\frac{\Vert {r}_{u}\Vert}{\mathrm{max}\left(1,\Vert u\Vert \right)}+\frac{\left|{f}^{T}x-{b}^{T}y+{u}^{T}w\right|}{\mathrm{max}\left(1,\left|{f}^{T}x\right|,\left|{b}^{T}y-{u}^{T}w\right|\right)}\le tol,$$

where

$$\begin{array}{c}{r}_{b}=Ax-b\\ {r}_{f}={A}^{T}y-w+z-f\\ {r}_{u}=x+s-u\end{array}$$

are the primal residual, dual residual, and upper-bound feasibility respectively, and

$${f}^{T}x-{b}^{T}y+{u}^{T}w$$

is the difference between the primal and dual objective values,
and *tol* is some tolerance. The sum in the stopping criteria measures the total relative errors in
the optimality conditions in Equation 8-7.

The algorithm begins by attempting to simplify the problem by removing redundancies and simplifying constraints. The tasks performed during the presolve step include:

Check if any variables have equal upper and lower bounds. If so, check for feasibility, and then fix and remove the variables.

Check if any linear inequality constraint involves just one variable. If so, check for feasibility, and change the linear constraint to a bound.

Check if any linear equality constraint involves just one variable. If so, check for feasibility, and then fix and remove the variable.

Check if any linear constraint matrix has zero rows. If so, check for feasibility, and delete the rows.

Check if the bounds and linear constraints are consistent.

Check if any variables appear only as linear terms in the objective function and do not appear in any linear constraint. If so, check for feasibility and boundedness, and fix the variables at their appropriate bounds.

Change any linear inequality constraints to linear equality constraints by adding slack variables.

If algorithm detects an infeasible or unbounded problem, it halts and issues an appropriate exit message.

The algorithm might arrive at a single feasible point, which represents the solution.

If the algorithm does not detect an infeasible or unbounded problem in the presolve step, it continues, if necessary, with the other steps. At the end, the algorithm reconstructs the original problem, undoing any presolve transformations. This final step is the postsolve step.

For simplicity, the algorithm shifts all lower bounds to zero.

While these preprocessing steps can do much to speed up the
iterative part of the algorithm, if the Lagrange multipliers are required, the preprocessing steps
must be undone since the multipliers calculated during the algorithm
are for the transformed problem, and not the original. Thus, if the
multipliers are *not* requested, this transformation
back is not computed, and might save some time computationally.

At a high level, the `linprog`

`'dual-simplex'`

algorithm
essentially performs a simplex algorithm on the *dual problem*.

The algorithm begins by preprocessing as described in Preprocessing. For details, see Andersen and Andersen [1] and Nocedal and Wright [7], Chapter 13. This preprocessing reduces the original linear programming problem to the form of Equation 8-4:

$$\underset{x}{\mathrm{min}}{f}^{T}x\text{suchthat}\{\begin{array}{c}A\cdot x=b\\ 0\le x\le u.\end{array}$$

*A* and *b* are transformed
versions of the original constraint matrices. This is the primal problem.

As explained in Equation 8-6, the dual problem is to find vectors *y* and *w*,
and a slack variable vector *z* that solve

$$\mathrm{max}{b}^{T}y-{u}^{T}w\text{suchthat}\{\begin{array}{c}{A}^{T}\cdot y-w+z=f\\ z\ge 0,\text{}w\ge 0.\end{array}$$

It is well known (for example, see [7]) that any finite solution of the dual problem gives a solution to the primal problem, and any finite solution of the primal problem gives a solution of the dual problem. Furthermore, if either the primal or dual problem is unbounded, then the other problem is infeasible. And if either the primal or dual problem is infeasible, then the other problem is either infeasible or unbounded. Therefore, the two problems are equivalent in terms of obtaining a finite solution, if one exists. Because the primal and dual problems are mathematically equivalent, but the computational steps differ, it can be better to solve the primal problem by solving the dual problem.

To help alleviate degeneracy (see Nocedal and Wright [7], page 366), the dual simplex algorithm begins by perturbing the objective function.

Phase 1 of the dual simplex algorithm is to find a dual feasible point. The algorithm does this by solving an auxiliary linear programming problem.

During Phase 2, the solver repeatedly chooses an entering variable and a leaving variable. The algorithm chooses a leaving variable according to a technique suggested by Forrest and Goldfarb [3] called dual steepest-edge pricing. The algorithm chooses an entering variable using the variation of Harris' ratio test suggested by Koberstein [5]. To help alleviate degeneracy, the algorithm can introduce additional perturbations during Phase 2.

The solver iterates, attempting to maintain dual feasibility while reducing primal infeasibility, until the solution to the reduced, perturbed problem is both primal feasible and dual feasible. The algorithm unwinds the perturbations that it introduced. If the solution (to the perturbed problem) is dual infeasible for the unperturbed (original) problem, then the solver restores dual feasibility using primal simplex or Phase 1 algorithms. Finally, the solver unwinds the preprocessing steps to return the solution to the original problem.

This section defines the terms *basis*, *nonbasis*,
and *basic feasible solutions* for a linear programming
problem. The definition assumes that the problem is given in the following
standard form:

$$\underset{x}{\mathrm{min}}{f}^{T}x\text{suchthat}\{\begin{array}{c}A\cdot x=b,\\ lb\le x\le ub.\end{array}$$

(Note that *A* and *b* are
not the matrix and vector defining the inequalities in the original
problem.) Assume that *A* is an *m*-by-*n* matrix,
of rank *m* < *n*, whose columns
are {*a*_{1}, *a*_{2}, ..., *a _{n}*}.
Suppose that $$\left\{{a}_{{i}_{1}},{a}_{{i}_{2}},\mathrm{...},{a}_{{i}_{m}}\right\}$$ is a basis for the column space
of

If *x* is a solution to *A·x* = *b* and
all the nonbasic variables in *x*_{N} are
equal to either their lower or upper bounds, *x* is
called a *basic solution*. If, in addition, the
basic variables in *x*_{B} satisfy
their lower and upper bounds, so that *x* is a feasible
point, *x* is called a *basic feasible solution*.

[1] Andersen, E. D., and K. D. Andersen. *Presolving
in linear programming*. Math. Programming 71, 1995, pp.
221–245.

[2] Applegate, D. L., R. E. Bixby, V. Chvátal
and W. J. Cook, *The Traveling Salesman Problem: A Computational
Study*, Princeton University Press, 2007.

[3] Forrest, J. J., and D. Goldfarb. *Steepest-edge
simplex algorithms for linear programming*. Math. Programming
57, 1992, pp. 341–374.

[4] Gondzio, J. "Multiple centrality
corrections in a primal dual method for linear programming." *Computational
Optimization and Applications*, Volume 6, Number 2, 1996,
pp. 137–156. Available at `http://www.maths.ed.ac.uk/~gondzio/software/correctors.ps`

.

[5] Koberstein, A. *Progress in the
dual simplex algorithm for solving large scale LP problems: techniques
for a fast and stable implementation*. Computational Optim.
and Application 41, 2008, pp. 185–204.

[6] Mehrotra, S. "On the Implementation
of a Primal-Dual Interior Point Method." *SIAM Journal
on Optimization*, Vol. 2, 1992, pp 575–601.

[7] Nocedal, J., and S. J. Wright. *Numerical
Optimization*, Second Edition. Springer Series in Operations
Research, Springer-Verlag, 2006.

Was this topic helpful?