Interior-point method for LP

From optimization
Jump to: navigation, search

Authors: John Plaxco, Alex Valdes, Wojciech Stojko. (ChE 345 Spring 2014)
Steward: Dajun Yue, Fengqi You
Date Presented: May 25, 2014


Introduction and Uses

Examples of logarithmic barrier functions.

Interior point methods are a type of algorithm that are used in solving both linear and nonlinear convex optimization problems that contain inequalities as constraints. The LP Interior-Point method relies on having a linear programming model with the objective function and all constraints being continuous and twice continuously differentiable. In general, a problem is assumed to be strictly feasible, and will have a dual optimal that will satisfy Karush-Kuhn-Tucker (KKT) constraints described below. The problem is solved (assuming there IS a solution) either by iteratively solving for KKT conditions or to the original problem with equality instead of inequality constraints, and then applying Newton's method to these conditions.

Interior point methods came about from a desire for algorithms with better theoretical bases than the simplex method. While the two strategies are similar in a few ways, the interior point methods involve relatively expensive (in terms of computing) iterations that quickly close in on a solution, while the simplex method involves usually requires many more inexpensive iterations. From a geometric standpoint, interior point methods approach a solution from the interior or exterior of the feasible region, but are never on the boundary.[1]

There are two important interior point algorithms: the barrier method and primal-dual IP method. The primal-dual method is usually preferred due to its efficiency and accuracy. Major differences between the two methods are as follows. There is only one loop/iteration in primal-dual because there is no distinction between outer and inner iterations as with the barrier method. In primal-dual, the primal and dual iterates do not have to be feasible.[3]

Barrier Method

For the barrier method algorithm, there a few approximations that must be made. Given a problem in the form of
Minimize f_0(x)
subject to  f_i(x) \leq 0
we must reformulate it to implicitly include the inequalities in the objective function. We can do this by creating a function that greatly increases the objective if a constraint is not met. Our conditions then changes to
minimize  f_0(x) + \sum_i^m I_-(f_i(x))
st  Ax=b
where  I_-(x)=\left|\begin{array}{cc}0& x \leq 0\\ \infty & x > 0\end{array}\right|

This problem, however, is not continuous. A modification can be made by approximating I_-(x) as a logarithm log(-x), which approaches infinity when x approaches 0 as we want, and makes all functions twice differentiable. We then put the logarithm over a variable that sets a level of accuracy for the approximation we make. Here we will call that variable t. We define
which blows up if any of our constraints are violated. Our LP problem now becomes
minimize  f_0(x) + \frac{1}{t}\phi(x)
st  Ax=b

This allows us to use Newton's method to follow what is called a Central Path, which is a series of points we iterate through that all satisfy the equality constraints Ax=b from the original problem, but give increasingly more optimized values for the objective function, with the inequality constraints f_i(x) not necessarily equal to 0.


Given strictly feasible x, t:=t^0>0,\mu >1, \epsilon<0
1. Compute x^*(t) by minimizing tf_0 + \phi subject to Ax = b, starting at x.
2. Update x:=x^*(t).
3. Quit if \frac{m}{t} \leq \epsilon, else
4. Increase t:=\mu t[3]


min  f(x) = x_1 - 2x_2
s.t.  1 + x_1 - (x_2)^2 >= 0

 x_2 >= 0

Use the Frisch barrier function to yield the unconstrained problem:

min  r(c,x) = x_1 - 2x_2 - clog(1+x_1-(x_2)^2) - clog(x_2)

For a specific parameter c, the first order necessary conditions for optimality are:

 1 - c/(1+x_1-(x_2)^2) = 0

 -2 + 2cx_2/(1+x_1-(x_2)^2) - c/x_2 = 0


 c = 1 + x_1 - (x_2)^2

substituting c in

 -2 + 2cx_2/(1+x_1-(x_2)^2) - c/x_2 = 0

we obtain

 (x_2)^2 - x_2 - (1/2)c = 0


 x_2 = (1 +/- sqrt(1+2c))/2 , of which the positive is the only feasible solution

Primal-Dual IP Algorithm

The primal-dual interior-point method can easily be understood by using the simplest NLP problem; one with only inequality constraints. Consider the following:

Number of iterations required for PD IP method to solve randomly generated standard LPs of different dimensions, with n = 2m. Error bars show standard deviation.

minimize f(x)

s.t. g_i(x)\ge 0 with i = 1,...,m.

We now introduce slack variables to turn all inequalities into non-negativities:

minimize f(x)

s.t. g(x)- s = 0 with s \ge 0.

The logrithmic barrier function is now introduced:

minimize  f(x) - \mu~  \sum_{i=1}^m\log(s_i)

s.t. h(x)-s=0

Now incorporate the equality constraint(s) into the objective function using Lagrange multipliers:

minimize f(x) - \mu~  \sum_{i=1}^m\log(s_i) - y^T(g(x) - s)

Next, set all of the derivatives equal to 0:

\nabla f(x) - \nabla g(x)^Ty = 0

- \mu W^{-1} e + y = 0

g(x) - s = 0


 \nabla f(x) - \nabla g(x)^Ty = 0

W Y e = \mu e

g(x) - s = 0

Utilize Newton's Method to determine the search directions,  \Delta x, \Delta s, \Delta y :

\begin{bmatrix} G(x,y) & 0 & -A(x)^T \\ 0 & Y & W \\ A(x) & -I & 0 \end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta s \\ \Delta y \end{bmatrix} = \begin{bmatrix} -\nabla f(x) + A(x)^T y \\ \mu e - W Y e \\ -g(x) + s \end{bmatrix}

where G(x,y) = \nabla^2 f(x) - \sum_{i=1}^m\ y_i \nabla^2 g_i(x)

and  A(x) = \nabla g(x)

Using the 2nd equation, we solve for \Delta s, the result of which is the reduced KKT system:

\begin{bmatrix} -G(x, y) & A^T(x) \\ A(x) & WY^{-1} \end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix} = \begin{bmatrix} \nabla f(x) - A^T(x) y \\ -g(x) + \mu Y^{-1} e \end{bmatrix}

From here, perform iterations:

 x^{k+1} = x^{k} + \alpha^k \Delta x^k

 s^{k+1} = s^{k} + \alpha^k \Delta s^k

 y^{k+1} = y^{k} + \alpha^k \Delta y^k


The Interior Point method approximates the constraints of a linear programming model as a set of boundaries surrounding a region. These approximations are used when the problem has constraints that are discontinuous or otherwise troublesome, but can me modified so that a linear solver can handle them. Once the problem is formulated in the correct way, Newton's method is used to iteratively approach more and more optimal solutions within the feasible region. Two practical algorithms exist in IP method: barrier and primal-dual. Primal-dual method is a more promising way to solve larger problems with more efficiency and accuracy. As shown in the figure above, the number of iterations needed for the primal-dual method to solve a problem increases logarithmically with the number of variables, and standard error only increases rapidly when a very large number of dimensions exist.


1. R.J. Vanderbei, Linear Programming: Foundations and Extensions (Chp 17-22). Springer, 2008.
2. J. Nocedal, S. J. Wright, Numerical optimization (Chp 14). Springer, 1999.
3. S. Boyd, L. Vandenberghe, Convex Optimization (Chp 11). Cambridge University Press, 2009