The Linear Programming Solver

The Interior Point Algorithm

The interior point LP algorithm in PROC OPTMODEL implements an infeasible primal-dual predictor-corrector interior point algorithm. To illustrate the algorithm and the concepts of duality and dual infeasibility, consider the following LP formulation (the primal):

\[ \begin{array}{rl} \mbox{min} & \mathbf{c}^\mr {T} \mathbf{x} \\ \mbox{subject to} & \mathbf{A} \mathbf{x} \ge \mathbf{b} \\ & \mathbf{x} \ge \mathbf{0} \end{array} \]

The corresponding dual is as follows:

\[ \begin{array}{rccccl} \mbox{max} & \mathbf{b}^\mr {T} \mathbf{y} & & & & \\ \mbox{subject to} & \mathbf{A}^\mr {T} \mathbf{y} & + & \mathbf{w} & = & \mathbf{c} \\ & & & \mathbf{y} & \ge & \mathbf{0} \\ & & & \mathbf{w} & \ge & \mathbf{0} \end{array} \]

where $\mathbf{y} \in \mathbb {R}^ m$ refers to the vector of dual variables and $\mathbf{w} \in \mathbb {R}^ n$ refers to the vector of dual slack variables.

The dual makes an important contribution to the certificate of optimality for the primal. The primal and dual constraints combined with complementarity conditions define the first-order optimality conditions, also known as KKT (Karush-Kuhn-Tucker) conditions, which can be stated as follows:

\[ \begin{array}{rclr} \mathbf{A}\mathbf{x} - \mathbf{s} & = & \mathbf{b} & \mr{(Primal~ Feasibility)} \\ \mathbf{A}^\mr {T} \mathbf{y} + \mathbf{w} & = & \mathbf{c} & \mr{(Dual~ Feasibility)}\\ \mathbf{W} \mathbf{X} \mathbf{e} & = & \mathbf{0} & \mr{(Complementarity)}\\ \mathbf{S} \mathbf{Y} \mathbf{e} & = & \mathbf{0} & \mr{(Complementarity)}\\ \mathbf{x},\, \mathbf{y},\, \mathbf{w}, \, \mathbf{s} & \ge & \mathbf{0} & \\ \end{array} \]

where $\mathbf{e} \equiv (1, \ldots , 1)^\mr {T}$ of appropriate dimension and $\mathbf{s} \in \mathbb {R}^ m$ is the vector of primal slack variables.

Note: Slack variables (the s vector) are automatically introduced by the algorithm when necessary; it is therefore recommended that you not introduce any slack variables explicitly. This enables the algorithm to handle slack variables much more efficiently.

The letters $\mathbf{X},$ $\mathbf{Y},$ $\mathbf{W},$ and $\mathbf{S}$ denote matrices with corresponding x, y, w, and s on the main diagonal and zero elsewhere, as in the following example:

\[ \mathbf{X} \equiv \left[ \begin{array}{cccc} x_1 & 0 & \cdots & 0 \\ 0 & x_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & x_ n \end{array} \right] \]

If $(\mathbf{x}^\ast , \mathbf{y}^\ast , \mathbf{w}^\ast , \mathbf{s}^\ast )$ is a solution of the previously defined system of equations representing the KKT conditions, then $\mathbf{x}^\ast $ is also an optimal solution to the original LP model.

At each iteration the interior point algorithm solves a large, sparse system of linear equations as follows:

\[ \left[ \begin{array}{cc} \mathbf{Y}^{-1}\mathbf{S} & \mathbf{A} \\ \mathbf{A}^\mr {T} & - \mathbf{X}^{-1}\mathbf{W} \end{array}\right] \left[ \begin{array}{c} \Delta \mathbf{y} \\ \Delta \mathbf{x} \end{array} \right] = \left[ \begin{array}{c} \Xi \\ \Theta \\ \end{array} \right] \]

where $\Delta \mathbf{x}$ and $\Delta \mathbf{y}$ denote the vector of search directions in the primal and dual spaces, respectively; $\Theta $ and $\Xi $ constitute the vector of the right-hand sides.

The preceding system is known as the reduced KKT system. The interior point algorithm uses a preconditioned quasi-minimum residual algorithm to solve this system of equations efficiently.

An important feature of the interior point algorithm is that it takes full advantage of the sparsity in the constraint matrix, thereby enabling it to efficiently solve large-scale linear programs.

The interior point algorithm works simultaneously in the primal and dual spaces. It attains optimality when both primal and dual feasibility are achieved and when complementarity conditions hold. Therefore it is of interest to observe the following four measures:

  • Relative primal infeasibility measure $\alpha $:

    \[ \alpha = \frac{\| \mathbf{Ax} - \mathbf{b} - \mathbf{s}\| _2}{\| \mathbf{b}\| _2 + 1} \]

  • Relative dual infeasibility measure $\beta $:

    \[ \beta = \frac{\| \mathbf{c} - \mathbf{A}^\textrm {T} \mathbf{y} - \mathbf{w} \| _2}{\| \mathbf{c}\| _2 + 1} \]

  • Relative duality gap $\delta $:

    \[ \delta = \frac{|\mathbf{c}^\textrm {T}\mathbf{x} - \mathbf{b}^\textrm {T}\mathbf{y}|}{|\mathbf{c}^\textrm {T}\mathbf{x}| + 1} \]

  • Absolute complementarity $\gamma $:

    \[ \gamma =\sum _{i=1}^ n x_ iw_ i + \sum _{i=1}^ m y_ i s_ i \]

where $\| v\| _2$ is the Euclidean norm of the vector v. These measures are displayed in the iteration log.