Linear Programming

Standard primal and dual problems

The standard linear programming problem is

\[ \text{(P)} \qquad \min cx \text{ s.t. } Ax=b;\ x\geq0. \]

Without loss of generality we can take $b\geq0$. We let $x^*$ be an optimal point.

The dual to the standard problem is

\[ \text{(D)} \qquad \max yb \text{ s.t. } yA\leq c \]

or alternatively

\[ \text{(D)} \qquad \max yb \text{ s.t. } yA+z=c;\ z\geq 0. \]

The variables $y$ are called dual variables and the variables $z$ are slack variables.

Theorem If $x$ is feasible for (P) and $y$ is feasible for (D), then

\[cx \geq yb.\]

Further, if $x^*$ is optimal for (P) and $(y^*,z^*)$ is optimal for (D), then

\[cx^* =y^*b \qquad \text{and} \qquad z^*\!\cdot\!x^*=0 .\]

The second condition is called complementary slackness.

Basic solutions

A basic solution is given by a set of column indices $B$ such that the square matrix $A_B$ formed by the $B$ columns of $A$ is nonsingular. Then

\[x_B=A_B^{-1}b, \ x_N=0; \quad y=c_BA_B^{-1}; \quad z_B=0, \ z_N=c_N-c_BA_B^{-1}A_N. \]


\[ cx = c_BA_B^{-1}b = c_B(A_B)^{-1}b = yb . \]

Suppose $x^*$ is optimal. Then $x_B\geq0$ and

\[cx = c_B x_B + c_N x_N = c_BA_B^{-1}(b-A_Nx_N)+c_Nx_N = c_BA_B^{-1}b + (c_N-c_BA_B^{-1}A_N)x_N = cx^* + (c_N-c_BA_B^{-1}A_N)x_N,\]

and hence $z_N=c_N-c_BA_B^{-1}A_N\geq0$.

Lower and upper bounds on variables

The constrained primal linear programming problem is

\[ \text{(CP)} \qquad \min cx \text{ s.t. } Ax=b;\ l\leq x\leq u . \]

Given lower bounded variables $x_L$ and upper bounded variables $x_U$, the problem becomes

\[ \min cx \text{ s.t. } Ax=b;\ x_L\geq l_L;\ x_U\leq u_U\]

and the dual problem is

\[ \text{(DCP)} \qquad \max\ y(b - A_Ll_L - A_Uu_U)+c_Ll_L + c_Uu_U \ \text{ s.t. }\ yA_L\leq c_L;\ yA_U\geq c_U . \]

Note that the objective function can be written $yb +(c_L - yA_L)l_L + (c_U - yA_U)u_U$

To define a basic solution, we need to split the variables $x$ into three sets, the basic variables $B$, the lower bounded variables $L$ and upper bounded variables $U$. A basic solution then satisfies

\[ x_L=l_L, \ x_U=u_U, \ x_B=A_B^{-1}(b-A_Ll_L-A_Uu_U); \quad y=c_B^{-1}A_B; \quad z_B=0,\ z_N=c_N-c_BA_B^{-1}A_N. \]

Lower and upper bounds on variables and constraints

The most general form of the linear programming problem is

\[ \text{(GCP)} \qquad \min cx \text{ s.t. } b_L \leq A_Cx \leq b_U;\ A_E = b_E;\ d_L\leq x\leq d_U . \]

Define auxiliary variables $w = w_C = A_C x$, and $w_E=A_E x = b_E$. The primal and auxiliary variables can be either basic , lower, upper or fixed. A lower or upper variable satisfies $x_N=d_{L/U}$ or $w_N=b_{L/U}$, and a fixed variable satisfies $w_E=b_E$ (alternatively, $w=b_L=b_U$.

After the introduction of auxiliary variables, the system can be written in the form (CP) above, with equations

\[ \begin{pmatrix} A_C & -I \\ A_E & 0 \end{pmatrix} \begin{pmatrix} x \\ w \end{pmatrix} = \begin{pmatrix} 0 \\ b_E \end{pmatrix} \]

If there are $m$ variables, $n_C$ inequality constraints and $n_E$ equality constraints, then the total number of basic variables $\#x_B+\#w_B=n_C+n_E$, since there are $n_C+n_E$ equations $A_Cx=w_C$ and $A_Ex=b_E$.

Given a basic solution, we can split the equations and variables into The basic variables $x_B$ are therefore given by

\[ \begin{pmatrix} A_{BB} & A_{BN} & -I & 0 \\ A_{NB} & A_{NN} & 0 & -I \\ A_{EB} & A_{EN} & 0 & 0 \end{pmatrix} \begin{pmatrix} x_B \\ x_N \\ w_B \\ w_N \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\b_E \end{pmatrix} \]

have the equations

\[ \begin{aligned} w_B &= A_{BB} x_B + A_{BN} x_N \\ A_{NB} x_B &= w_N - A_{NN} x_N \\ A_{EB} x_B &= b_E - A_{EN} x_N \end{aligned} \]

The basic variables $x_B$ are therefore given by

\[ x_B = {\begin{pmatrix}A_{NB}\\A_{EB}\end{pmatrix}}^{-1} \begin{pmatrix}w_N - A_{NN} x_N \\ b_E - A_{EN} x_N \end{pmatrix} \]


\[ B = {\begin{pmatrix}A_{NB}\\A_{EB}\end{pmatrix}}^{-1} = \begin{pmatrix} B_{BN} & B_{BE} \end{pmatrix} = {\begin{pmatrix}A^{-\dag}_{BN} & A^{-\dag}_{BE}\end{pmatrix}} \]

Since the objective function is given by $cx$, we have

\[ \begin{aligned} cx &= c_B x_B + c_N x_N = c_B A_{BN}^{-\dag} (w_N - A_{NN} x_N) + c_B A_{BE}^{-\dag}(b_E - A_{EN} x_N) \\ &= c_B A_{BE}^{\dag} b_E \ + \ c_B A_{BN}^{-\dag} w_N \ +\ \bigl(c_N - c_B (A_{BN}^{-\dag} A_{NN} + A_{BE}^{-\dag} A_{EN} ) \bigr) x_N \\ &= c_B A_{BE}^{\dag} b_E \ + \ c_B A_{BN}^{-\dag} w_N \ +\ (c_N - c_B A_{{}^N_E B}^{-1} A_{{}^N_E N} ) x_N \end{aligned} \]

Set dual variables

\[ \begin{aligned} y_N &= c_B A_{BN}^{-\dag} \\ y_E &= c_B A_{BE}^{-\dag} \end{aligned} \]


\[ \begin{aligned} cx &= c_B A_{BE}^{-\dag} b_E \ + \ y_N w_N \ +\ (c_N - (y_N A_{NN} + y_E A_{EN}) ) x_N \\ &= c_B A_{BE}^{-\dag} b_E \ + \ y_N w_N \ +\ (c_N - y_{{}^N_E} A_{{}^N_EN} ) x_N \end{aligned} \]

This shows how the objective changes with changes in $x_N$ and $w_N$. The main difficulty performing an update is that the number of basic $x$ variables may change during an algorithm.

Primal and dual feasibility problems

Certificates of infeasibility / Farka's Lemma

Basic feasibility testing

Given a partition of variables into upper, lower and basic, we can test for a solution of the constrained feasibility problem.

Robust feasibility problems

Due to numerical errors, the solution to a feasibility problem may be inaccurate. It is therefore important to be able to validate solutions to feasibility problems. One approach is to use rational arithmetic, but this is frequently too expensive. The alternative is to use interval arithmetic. For this, we need our problems to be robustly solvable.

Note that robustly solving an inequality is straightforward; we merely change $\leq$ to $<$. To robustly solve a system of linear equalities, we need to choose a basis and express the basic variables in terms of the non-basic variables.

Robust basic solutions

Given a basis $B$, and if necessary lower and upper variables $L$ and $U$, we can attempt to solve robust feasibility problems by finding a robust basis. We can also set the nonbasic variables to the exact values, as long as the basic variables robustly satisfy the constraints. In order to check that a basis gives a solution to the robust dual feasibility problem,

TODO: It is not clear whether a robustly solvable problem has a robust basic solution.

Certificates of robust (in)feasibility

A robust certificate for the primal feasibility problem $Ax=b;\ x\geq0$ is a base $B$ such that $A_B$ is nonsingular, and an $x_N>0$ such that $-A_B^{-1}A_Nx_N>0$; this is equivalent to $Ax=b;\ x>0$. A robust certificate of infeasibility is a point $y$ such that $A^Ty<0$ and $b^Ty>0$.

Theorem Suppose $Ax=b;\ x\geq0$ is robust. Then either $A$ has full row rank and there exists $x>0$ such that $Ax=b$, or there exists $y$ such that $A^T<0$ and $b^Ty>0$.

A robust certificate for the dual feasibility problem $A^Ty\leq c$ is a point $y$ such that $A^Ty<c$. A robust certificate of infeasibility is a base $B$ such that $A_B$ is nonsingular, and a vector $x>0$ such that $Ax=0$ and $c^Tx<0$. We prove $Ax=0$ by setting $x_B=-A_B^{-1}A_Nx_N$, so that $Ax=A_Bx_B+A_Nx_N=-A_BA_B^{-1}A_Bx_N+A_Nx_n=0$.

Theorem Suppose $A^Ty\leq c$ is robust. Then either there exists $y$ such that $A^Ty<c$, or $A$ has full row rank and there exists $x>0$ such that $Ax=0$ and $c^Tx<0$.

Proof If $A^Ty<c$, then this holds also for perturbations of $y$. If $Ax=0$ and $A_B$ is nonsingular, then $x_B=-x_NA_NA_B^{-1}$. Perturbing $A,b$ and keeping $x_N$ constant, we obtain a perturbation of $x_B$, and hence a certificate for the perturbed problem.
Conversely, suppose the problem is robustly solvable. Then the problem $A^Ty\leq c-\epsilon p$ is solvable for $p>0$ and $\epsilon$ sufficiently small. Hence there exists $y$ such that $A^Ty<c$.
Suppose the problem is robustly unsolvable Let $P$ be a matrix with all positive entries. Since $A^Ty\leq c$ is robustly unsolvable, $(I+\epsilon P)^T A^Ty\leq (I+\epsilon P)^T c$ is unsolvable for some $\epsilon>0$. Then there exists $x_\epsilon$ such that $A(I+\epsilon P)x_\epsilon=0$, $x_\epsilon\geq0$ and $(I+\epsilon P)x_\epsilon c<0$. Then if $x=(I+\epsilon P)x_\epsilon$, then $Ax=0$, $c^Tx<0$ and $x>0$ since $x_\epsilon\geq0;\ x_\epsilon\neq0$ and $(I+\epsilon P)>0$.

Converting robust feasibility problems to feasibility problems

To solve the robust primal feasibility problem,

\[ \text{(RPF)} \qquad Ax=b;\ x>0;\ A_B\ \text{nonsingular} \]

we choose a vector $p>0$ and consider the problem

\[ \min -s \text{ s.t. } Ax + Ap\,s = b;\ x,s\geq0 .\]

Let $\hat{x}^T=(x\;s)$, $\hat{A}=(A\;Ap)$ and $\hat{c}^T=(0\;\mbox{}-\!1)$. Then we obtain the standard primal optimisation problem $ \min \hat{c}^T\hat{x} \text{ s.t. } \hat{A}\hat{x} = b;\ \hat{x}\geq0 . $ If the optimal value is negative, then we have found $x^*,s^*$ such that $A(x^*+ps^*)=b; \ x^*\geq0,\ s^*>0$, so taking $\tilde{x}=x^*+ps^*$, we have $A\tilde{x}=b;\ \tilde{x} = x^*+s^*p \geq s^*p > 0$. If the optimal value is non-negative, then we can attempt to solve the dual robust optimisation problem $ \max yb \text{ s.t. } yA<0 $. Since if this problem is feasible, it has unbounded solutions, we can instead choose $q>0$ look for a positive optimal value of

\[ \max yb \text{ s.t. } Ay\leq-q . \]

To solve the robust dual feasibility problem,

\[ \text{(RDF)} \qquad A^T y + z = c; \ z>0 \quad \text{or} \quad A^Ty < c, \]

we choose $q>0$ and consider the problem

\[ \max t \text{ s.t. } A^T y + qt \leq c . \]

Let $ \hat{A}^T = (A^T \; q)$, $\hat{b}^T = (0^T \; 1)$ and $\hat{y}^T = ( y^T\;t )$. Then we obtain the standard dual optimisation problem $ \max \hat{b}^T \hat{y} \text{ s.t. } \hat{A}^T \hat{y} \leq c; $ If the optimal value is positive, then we have found $y^*,t^*$ such that $A^Ty^*\leq c-t^*q < c$. If the optimal value is zero or negative, then we can attempt to solve the primal robust feasibility problem $ \min b^T x \text{ s.t. } A x = 0, \ q^T x = 1, \ x\geq 0 . $ However, even if the dual feasibility problem is unsolvable (i.e. $t^* < 0$), the primal may become solvable by a perturbation of $A$. We therefore consider a robust version $ \min b^T x \text{ s.t. } A x = 0,\ x>0 . $ and introduce $p>0$ to make a problem $ \min b^T x \text{ s.t. } A x = 0,\ x-p\geq0 . $ Note that if this problem has negative value, then we have found $x^*$ such that $b^Tx^*<0$, $Ax^*=0$ and $x^*>0$, which implies that the original dual problem has no solution. Taking $\tilde{x} = x-p$, we obtain

\[ \min b^T \tilde{x} + b^T p \text{ s.t. } A\tilde{x} = -Ap,\ \tilde{x}\geq 0 . \]

The simplex algorithm

Suppose we wish to update a basis of the standard linear programming problem.

In general, if $A' = A + a \alpha$, then

\[ (A')^{-1} = A^{-1} - \frac{ A^{-1} a \; \alpha A^{-1} }{1+\alpha A^{-1} a} \]

The dual revised simplex algorithm

Suppose we wish to update a dual feasible basis of the standard linear programming problem; note that the reduced costs satisfy satisfy $z_N\geq0$.

The simplex algorithm with constraints

Suppose we wish to update a basis of the standard linear programming problem.

Algorithms for feasibility

Remark: Since we do not assume the existence of $\pm\infty$ in our number types, we use $l=0,\ u=-1$ for the constraint $x\geq0$; this is the only unbounded constraint we allow.

See Chvatal [Chapter 8, pp 129] for more details on constrained feasibility.

Efficiency of the simplex algorithm

For a linear programming problem of standard form, with $A$ an $m\times n$ matrix, the number of iterations of the simplex algorithm for practical problems grows approximately as $m\log n$.

The homogeneous self-dual algorithm

The homogeneous self-dual algorithm is a method for simultaneously computing optimality and feasibility. Consider the primal and dual problems

\[ \text{(P)}\ \min c^T x \mid Ax=b,\ x\geq 0; \qquad \text{(D)}\max b^Ty \mid A^Ty\leq c \]

The homogeneous self-dual problem is

\[ \text{(HSD)} \begin{gathered} \min (n+1)\theta \\ Ax - b\tau + (b-A\textbf{1})\theta = 0 \\ -A^Ty +c\tau - (c-\mathbf{1})\theta\geq 0 \\ b^Ty - c^T x + (c^T\mathbf{1}+1)\theta \geq 0 \\ -(b-A\mathbf{1})^Ty + (c-\mathbf{1})^Tx - (c^T\mathbf{1}+1)\tau = -(n+1) \\ x_j,\tau\geq 0; \ y_i,\theta\in\R . \end{gathered} \]

A feasible point is always given by $x=\mathbf{1},\ y=0,\ \tau=1,\ \theta=1$. When $\tau=1,\ \theta=0$, the first three equations reduce to $Ax=b$, $A^Ty\leq c$ and $c^Tx=b^Ty$. If $z$ and $\kappa$ denote the slack variables, the final constraint becomes

\[ \mathbf{1}^Tx + \mathbf{1}^Tz + \tau + \kappa - (n+1)\theta = (n+1) \]

We have the following result:

To relate (HSD) to (P) and (D), we have:

To related (HSD) to the central path, we have:

A potential function is given by

\[ \psi_{n+1+\rho}(x,\tau,z,\kappa) = (n+1+\rho)\log(x^Tz+\tau\kappa)-\sum_j \log(x_jz_j) -\log(\tau\kappa) \]

See: David G. Luenberger and Yinyu Ye "Linear and nonlinear programming".

Feasibility problems for geometric operations

In the Geometry module, we need to solve the following linear programming problems to test intersection.

\[ \begin{array}{|l||c|c|c|c|}\hline &\text{Polyhedron}&\text{Polytope}&\text{Zonotope}\\\hline\hline \text{Point} & Ap\leq b & p=Vs;\ 1\!\cdot\!s=1;\ s\geq0 & p=c+Ge;\ -1\leq e\leq1 \\\hline \text{Rectangle} & Ax\leq b;\ l\leq x\leq u & x=Vs;\ 1\!\cdot\!s=1;\ l\leq x\leq u;\ s\geq0 & x=c+Ge;\ l\leq x\leq u; \ -1\leq e\leq1 \\\hline \text{Zonotope} & A(c+Ge)\leq b;\ -1\leq e\leq 1 & Vs=c+Ge;\ 1\!\cdot s=1;\ -1\leq e\leq1;\ s\geq0 & c_1+G_1e_1=c_2+G_2e_2;\ -1\leq e_1,e_2\leq1 \\\cline{0-3} \text{Polytope} & AVs\leq b;\ 1\!\cdot\!s=1;\ s\geq0 & V_1s_1=V_2s_2;\ 1\!\cdot s_1=1;\ 1\cdot s_2=1;\ s_1,s_2\geq0 \\\cline{0-2} \text{Polyhedron} & A_1x\leq b_1;\ A_2x\leq b_2 \\\cline{0-1} \end{array} \]

We notice that by introducing slack variables, we can convert all problems into a standard linear programming problem with constraints.

We can convert the standard dual feasibility problem into a primal linear programming problem $\min b^Ty\text{ s.t. } A^Ty=0$, but it is not so straightforward to convert a constrained dual feasibility problem into its dual. Instead we add slack variables and solve $ Ax+z=b;\ l\leq x\leq u$. We can use the reduced simplex algorithm to take advantage of sparseness.

Feasibility using Interior Point Methods

Constrained Dual Feasibility

We shall be most concerned with the following problem:

Without the bounds on $y$, it is easy to show that no solution exists if there exists $x$ such that $Ax=0$, $e^Tx=1$, $x\geq0$ and $c^Tx<0$.

To solve the problem, introduce extended variables $x_L,x_U$, slack variables $s,s_L,s_U$ and a value $v$. The original (dual) feasibility problem becomes

\[ \max v \text{ s.t. } A^Ty + et + s = c; \ -y+t+s_L = -l; y+t+s_U = u; s,s_L,s_U \geq 0 \]

whereas the prime feasibility problem becomes

\[ \min c^Tx + l^Tx_L + u^Tx_U \text{ s.t. } Ax - x_L + x_U = 0; e^Tx+e^Tx_L+e^Tx_U = 1; x,x_L,x_U \geq 0\]

. The original problem is feasible if the optimum value $t$ is positive, and infeasible if it is negative.

Robust primal and dual basic feasible solutions are easily found. Take $x=1$, and $x_L,x_U>0$ such that $Ax=x_L-x_U$. Finally, scale so that the sum is equal to $1$. Take $y=0$, $t=\min(c_i,-l_i,u_i)-1$ and $s=c-et$, $s_L=-l-t$ and $s_U=u-t$

We can now apply the standard interior point update. Although the matrix $A$ has been augmented, we can still recover a problem such that we only need to invert a matrix with only one extra row. Hence the dimension of the matrix inversion problem is equal to the dimension of the underlying space plus one.

General feasibility problem

Consider feasibility of the linear programming problem

\[ Ax=b; \quad Cx\leq d . \]

Special cases of this occur in, for example testing whether the image of a constraint set intersects another constraint set.

Introduce slack variables $r$ to transform the problem to

\[ Ax=b; \quad Cx+r=d; \quad r\geq 0 . \]

The duality conditions are

\[ A^T y + C^T s = 0 \]

and the value function is

\[ b^T y + d^T s . \]

The complementary slackness conditions are

\[ r\cdot s = 0; \quad r,s\geq 0 . \]

Consider a step updating the variables in an attempt to solve the complementarity system. The Newton equations are

\[ \left(\begin{matrix} A & 0 & 0 & 0 \\ C & 0 & I & 0 \\ 0 & A^T & 0 & C^T \\ 0 & 0 & S & R \end{matrix} \right) \left(\begin{matrix} \Delta x \\ \Delta y \\ \Delta r \\ \Delta s \end{matrix}\right) = \left(\begin{matrix} 0 \\ \rho_r \\ \rho_ s \\ \tau \end{matrix}\right) \]

where $\rho_r$ and $\rho_s$ are the residuals $\rho_r=Cx+r-d$ and $\rho_s=A^Ty + C^T s$ which are typically zero.

Assuming the residuals are zero, we have

\[\Delta r=-C\Delta x, \]

\[\Delta s = R^{-1}(\tau-S\Delta r) = R^{-1}\tau + R^{-1}SC\Delta x\]

and hence $A^T\Delta y + C^T(R^{-1}\tau + R^{-1}SC\Delta x) = 0$ so $C^TR^{-1}SC\Delta x + A^T \Delta y = -C^TR^{-1}\tau $. Combined with the equation $A^T\Delta x=0$, we obtain the system

\[ \left(\begin{matrix} C^TR^{-1}SC & A^T \\ A & 0 \end{matrix}\right) \; \left(\begin{matrix} \Delta x \\ \Delta y \end{matrix}\right) = \left(\begin{matrix} -C^TR^{-1}\tau \\ 0 \end{matrix}\right) . \]

The left-hand matrix is symmetric, so can be (fairly) easily factorised. Further, the $A$ matrix is independent of the problem, so the initial part of the factorisation can be pre-computed.

Indeed, it may be better to completely eliminate $A$ at the start of the problem by explicitly solving for some variables of x. If we have a square submatrix $A_B$ of $A$ which is invertible, then

\[x_B = A_B^{-1} (b-A_Nx_N)\]

and the inequalities for $x$ become

\[(C_N-C_BA_B^{-1}A_N)x_N \leq d-A_B^{-1}b\]


 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines