Interior-Point Methods of
Linear Programming

Mathematica code is presented for
solving the linear programming problem
max(cx | Ax = b,  x ≥ 0) using an
interior-point projective scaling method.

by Harold Greenberg

The linear programming (LP) problem may be written in
the standard form

maximize    cx
                                      subject to    Ax = b, x ≥ 0

where c is an n-row vector, b is an m-column vector, and A is an, m × n  matrix, m <  n. The elements of c, b, and A are given constants. The problem is to find the optimal
nonnegative values of the components of
n-column vector x.

The simplex algorithm is the method of choice for
many LP applications. The Mathematica function
LinearProgramming uses the simplex algorithm to
solve the problem minimize cx subject to
Ax ≥ b, x ≥ 0. The standard form given
above can be easily converted to the necessary 
Mathematica format.

The constraints Ax=b, x≥0 define a polyhedron in
n-dimensional space. The solution point to
the LP problem occurs at a vertex of the polyhedron.
The simplex algorithm succeeds in determining the
optimal point by starting at a vertex and moving from
vertex to vertex, so that the value of cx is
generally increasing, until the optimal point is
reached. The starting vertex is usually found by a
simple variant of the algorithm. See any text on LP
for the details on the simplex algorithm.

In recent years, the LP problem has also been solved
by interior-point} methods. The methods succeed
in getting close to an optimal point of the polyhedron
by starting at a point in the interior of the
polyhedron and moving from interior point to interior
point, in direction of increasing objective, until a
point is reached near an optimal point. If a starting
interior point is not known, the standard form of the
problem is altered and an arbitrary initial point is
taken. The altered model will then eventually produce
a point close to an optimal one of the original model.
Both approaches will be given below.

Theoretically, an interior-point method is more
efficient than the simplex algorithm on problems where
any of m, n and the values of c, A, b are
large. It appears, however, that computational
comparisons have been mixed [Tomlin 1988]. Interior
methods complement the simplex algorithm on problems
that the simplex algorithm handles well. The interior
methods are most applicable to models that lead to
lengthy solution  times for the simplex algorithm.

This paper presents an interior-point method,
written in very few lines of Mathematica  code,
that follows the projective scaling approach
[Karmarkar 1984]. Other interior-point methods appear
in the literature [e.g., see Megiddo 1988].

Projective Scaling with a Given Feasible Solution

We assume that there is a known starting interior
feasible solution x, to the constraint polyhedron,
with all positive values. If such a feasible solution
is not known, the problem can be solved by taking an
arbitrary x point, as discussed in the next section.

As an example of the method, assume that the constants
of the problem are input as

c = {2, 3, 0, 0};
A = {{1, 1, 1, 0}, {1,5,0,1}};
b = {2, 5};

representing the problem: Determine
(Subscript[x, 1]) ≥ 0, (Subscript[x, 2]) ≥ 0, (Subscript[x, 3]) ≥ 0, (Subscript[x, 4]) ≥ 0 that
maximize 2x_1+3x_2 subject to (Subscript[x, 1]) + (Subscript[x, 2]) + (Subscript[x, 3])= 2,
(Subscript[x, 1 ])+ 5 (Subscript[x, 2]) + (Subscript[x, 4])= 5.

We continue with

d = Dimensions[A]; m = d[[1]]; n = d[[2]];

For the example, a starting interior solution
is easily determined as

x = {.5, .5, 1, 2};

The projective scaling method assumes there is an
upper bound, B, to the sum of the x variables.
Hence,  (Subscript[x, 1 ])+ (Subscript[x, 2]) + ... + (Subscript[x, n]) < B; a slack variable
(Subscript[x, n  + 1]) ≥ 0 is included to achieve an
additional constraint  (Subscript[x, 1 ])+ (Subscript[x, 2]) + ... + (Subscript[x, n])+(Subscript[x, n  + 1]) = B. The value of
B is selected to be large enough so that optimal
x is not cut off by the bound. (If (Subscript[x, n  + 1])
approaches zero during the calculations, then B is
not large enough. A larger B is then picked;
(Subscript[x, n  + 1])is increased accordingly. Matrix A, in
storage, is changed to the new A, as below,
and the method is continued from the current modified x.)

Next, taking x as having n+1 components by
including (Subscript[x, n  + 1]), the original equations are turned
into a homogeneous form by Ax = (1/B)(ex)b,
where e is a row-vector} of n+1 ones.

Thus, to update the initial x and the coefficients
in the new n+1 dimensional problem, with B\/
having the value 10, form

B = 10;
AppendTo[x, B - Plus@@x];
AppendTo[c, 0];
A = Append[(B/b)*Apply[Join,
     Transpose[{A,Table[{0},{m}]}],
     {1}]-1, Table[1,{n+1}]];

Interior-point value of x is changed to
include the value of (Subscript[x, n  + 1]) > 0. The zero appended to
c is the coefficient of (Subscript[x, n  + 1]).  Matrix A is
changed to new matrix A by:

a)  adding column $n+1$ of m zeroes to A,
representing the coefficients of $x_{n+1}$;

(b) using the bound B to change to the
coefficients of the homogeneous equations;

(c) appending a row of $n+1$ ones to complete
the bottom row of A.

In the example, with B = 10, x, c, and
A become

x = {0.5, 0.5, 1, 2, 6}
                 c  = {2, 3, 0, 0, 0}
                 A = {{4, 4, 4, -1, -1}, {1, 9, -1, 1, -1},{1, 1, 1, 1, 1}}.

The main part of the method follows. Consider any
positive solution to Ax = b, where A is the updated
matrix, and b is now an (n+1)-vector} of n
zeroes with last component bound B.  This new b
is not used in the code continuation since the
constraints will be automatically satisfied. At each
step, the variables x are transformed by x = Dy,
where D is a diagonal matrix whose values along the
main diagonal are consecutively the current positive
feasible solution to Ax = b. The LP problem becomes

maximize    cDy
                    subject to    ADy=b, y ≥ 0

where c, A and b are the new values developed
in the code. A feasible y, for the positive feasible
x along the main diagonal of D, is then
y = (x/D). Thus, y is composed of a vector of
n+1 ones. The steepest ascent direction from
feasible y is along the gradient cD. A move
from current y to next value \(TraditionalForm\`y\&_\), in the
gradient direction, must be made so that \(TraditionalForm\`AD\ y\&_ = b\). To insure the constraints are satisfied, the
gradient is modified by the projection matrix

\(TraditionalForm
    \`P = i - AD\^\(\(T(AD\ AD\^T)\)\^{\(-1\)}\)\ AD\)

where I is the (n+1)×(n+1)identity matrix.
The next feasible solution is \(TraditionalForm\`y\&_\) as obtained as
\(TraditionalForm\`y\&_ = y + g\ α\) for constant  α and modified
gradient vector g=PcD. The value of α
is chosen so that the components of \(TraditionalForm\`y\&_\) remain
positive and far enough away from zero. Expanding
ADy with the above projection P gives \(TraditionalForm\`AD\ y\&_ = b\),
so that \(TraditionalForm\`\(y\ \)\&_\)remains feasible. The next feasible
x is given by \(TraditionalForm\`x\&_ = D\ y\&_\). This new feasible
\(TraditionalForm\`x\&_\) remains positive and replaces the current
one; the main part of the code is repeated until the
feasible x values are virtually unchanged. Hence, in
the code the modified gradient vector, g, is given
by

a = N[Transpose[A]];
g := 
(P = IdentityMatrix[n+1]-(a*x).PseudoInverse[a*x]; 
      P.(c*x));

Matrix A is changed to its transpose for
P; the N function insures that
real numbers are used to speed the calculations of
P. The function PseudoInverse gives the rest
of the P calculation. The modified gradient is then
given by P.(c x). The products a x and c x are
the same as (AD^T) and cD, respectively.
Finally,

While[g.g > 0.0001, x = x*(1 - 0.9*g/Min[g])];
Print["x=", Take[x,n], "z=", Take[c,n].Take[x,n]];

The While statement is the iterating step; the
modified gradient, g, and feasible x are
calculated repeatedly. The iterations stop when there
is little change in the feasible x values, which
occurs when g is small with g.g ≤ ϵ. Here
the choice was  ϵ=0.0001, but other values
of ϵ will do depending on the accuracy
required. If g.g > ϵ, x(1+ α g) is the
next value for x. The α value taken is
α = -0.9/Min[g]. The coefficient 0.9
prevents any component of x from becoming zero. The
final x and objective value z=cx are printed out
to end the method. In the example problem, four new
feasible x values are calculated; the output is
x={1.24892, 0.749284, 0.00179152, 0.00465532},
z=4.7457.

Figure 1 shows the feasible ((Subscript[x, 1]),(Subscript[x, 2])) region and two
successive ((Subscript[x, 1]),(Subscript[x, 2])) values calculated starting
from (0.5, 0.5); these are:  (0.918056, 0.776389) and
(1.25559, 0.713855). The values are seen to approach
(Subscript[x, 1])=1.25,  (Subscript[x, 2])=0.75 with objective value z=4.75.

The steps of the above coding of the method are
written so that the description of the method and the
code coincide. An overall Mathematica
code is given by Program Listing 1 that is somewhat
faster and preserves the original values of c, A, x.

Program Listing 1

InteriorLPMax1[c_, A_, b_, B_, x_] :=
  (d = Dimensions[A]; 
   m = d[[1]]; 
   n = d[[2]];
  x1 =  Flatten[{x, B-Plus@@x}];
  c1 = Flatten[{c, 0}];
  A1 = Append[B/b*Apply[Join, 
        Transpose[{A, Table[{0}, {m}]}], 
        {1}]-1, Table[1,{n+1}]];
   a = N[Transpose[A1]];
  g := (IdentityMatrix[n+1]-
         (a*x1).PseudoInverse[a*x1]).(c1*x1);
  While[g.g > .0001, x1 *= 1-.9*g/Min[g]];
  Print[''x='', Take[x1, n],
        ''z='', c.Take[x1,n]])

Projective Scaling with
No Apparent Feasible Solution

When no feasible \mbox{interior-point} solution is
apparent, the following \mbox{well-known} change
may be made; define a \mbox{big-M} type of problem as

maximize      cx - M (Subscript[x, n + 1])
subject to A x+(Subscript[x, n + 1])(b-A \(TraditionalForm\`\(x\^'\)\))=b,
x ≥ 0, (Subscript[x, n + 1]) ≥ 0,

where M is a sufficiently large real number,
(Subscript[x, n + 1]) is a scalar, \(TraditionalForm\`\(x\^'\)\) is
any n-vector with all positive components.
Thus, x is an optimal solution to the original
LP problem if, and only if, x, (Subscript[x, n + 1])= 0 is an
optimal solution to this big-M problem. The
idea is to solve by the interior-point
method. Obviously, x=\(TraditionalForm\`\(x\^'\)\), (Subscript[x, n + 1])=1 is a
feasible interior-point solution to the
constraints. The large value of M forces (Subscript[x, n + 1])
to approach zero. (If the method ends with (Subscript[x, n + 1])
not approaching zero, then the code must be restarted
from the current x by increasing M and changing c
as below.)

Taking \(TraditionalForm\`\(x\^'\)\) equal to a vector of ones works well.
That is the approach taken here. Hence, to solve the
big-M problem, we follow the code given
above. We repeat with the example problem

c = {2, 3, 0, 0};
A = {{1, 1, 1, 0}, {1, 5, 0, 1}};
b = {2, 5};
d = Dimensions[A]; 
m = d[[1]]; 
n = d[[2]];

The feasible n-dimensional x is given by

x = Table[1, {n}];

After allowing for (Subscript[x, n + 1]), the sum of the
x variables is again bounded by B. With
(Subscript[x, n + 2]) ≥ 0 as the slack variable, the resulting
equation allows for the change to the homogeneous form
for the big-M problem. For the example problem,
with B and M as having the values 10 and 20,
respectively, the initial x and the coefficients
are updated. Thus,

B = 10; 
M = 20;
A = Append[(B/b)*Apply[Join,
     Transpose[{A,Partition[b-A.x,1],
      Table[{0},{m}]}],
      {1}]-1, Table[1,{n+2}]];
x = Flatten[{x, 1, B-Plus@@x-1}];
c = Flatten[{c, -M, 0}];

The new matrix A has column
n+1 as the coefficients of (Subscript[x, n + 1]) as produced by
the Partition function and column n+2 from the
zero coefficients of the slack variable (Subscript[x, n + 2]), all
before using the bound B to change to the
coefficients of the homogeneous equations. A row of
n+2 ones is then appended as a bottom row for A.
The value of x is changed to include initial
(Subscript[x, n + 1]) and the value of (Subscript[x, n + 2]) > 0. The -M and
0 joined to c are the coefficients of (Subscript[x, n + 1]) and
(Subscript[x, n + 2]), respectively.
In the example, A,  x, and c become

A = {{4, 4, 4, -1, -6, -1}, {1, 9, -1, 1, -5, -1},{1, 1, 1, 1, 1, 1}}
x = {1, 1, 1, 1, 1, 5}
c = { 2, 3, 0, 0, -20, 0}.

The main part of the method continues as before with

a = N[Transpose[A]];
g := (IdentityMatrix[n+2]-
      (a*x).PseudoInverse[a*x]).(c*x);
While[g.g>.0001, x*=1-.9*g/Min[g]];
Print[''x='', Take[x, n],
      ''z='', Take[c, n].Take[x, n]];

The example problem produces
x={1.24934, 0.749407, 0.00134252, 0.00380264}
with z = 4.7469.

Figure 2 shows the feasible ((Subscript[x, 1]), (Subscript[x, 2])) region and
four successive  ((Subscript[x, 1]), (Subscript[x, 2])) values calculated starting
from infeasible (1.0, 1.0); these are: (0.763251,
0.680977), (0.921308, 0.655407), (1.32207, 0.641855),
(1.24838, 0.732185). The values are seen to approach
(Subscript[x, 1])=1.25, (Subscript[x, 2])=0.75 with objective value z=4.75.

Figure 2. ((Subscript[x, 1]), (Subscript[x, 2])) values starting at arbitrary (1, 1) and then approaching
(1.25, .75).

An overall Mathematica  code is given by Program Listing 2.

Program Listing 2:

InteriorLPMax2[c_, A_, b_, B_, M_] :=
 (d=Dimensions[A]; m=d[[1]]; n=d[[2]];
  x = Table[1, {n}];
 A1 = Append[(B/b)*Apply[Join,
       Transpose[{A,
        Partition[b-A.x,1],
        Table[{0},{m}]}],
       {1}]-1, Table[1,{n+2}]];
 x1 = Flatten[{x, 1, B-Plus@@x-1}];
 c1 = Flatten[{c, -M, 0}];
  a = N[Transpose[A1]];
  g:= (IdentityMatrix[n+2]-
       (a*x1).PseudoInverse[a*x1]).(c1*x1);
 While[g.g > .0001, x1 *= 1-.9*g/Min[g]];
 Print[''x='', Take[x1, n],
       ''z='', c.Take[x1,n]])

REFERENCES

Karmarkar N. A New Polynomial Time Algorithm for
Linear Programming. Combinatorica 4 (1984), 373-395.

Megiddo N.  (ed.). Progress in
Mathematical Programming.  Interior-Point  and
Related Methods.Springer-Verlag, New York (1988).

Tomlin J. A.  A Note on Comparing Simplex
andInterior-Point Methods for Linear
Programming. In N. Megiddo, loc. cit., 91-103.

About the Author

Harold Greenberg uses Mathematica in his work on networks, cryptography, and 
mathematical programming. He is the author of the book ``Integer Programming.'' 
In his spare time he writes humorous short  stories.

Harold Greenberg
Dept. of Statistics and Computer Information
Systems
Baruch College/CUNY
17 Lexington Avenue
New York, NY 10010
dlobb@cunyvm.cuny.edu