Dantzig G., Thapa M. Linear programming. Vol.1. Introduction
Подождите немного. Документ загружается.
This page intentionally left blank
CHAPTER 4
INTERIOR-POINT
METHODS
An interior-point algorithm is one that improves a feasible interior solution point of
the linear program by steps through the interior, rather than one that improves by
steps around the boundary of the feasible region, as the classical Simplex Algorithm
does.
Just as there are many variants of the Simplex Method (which we refer to as
pivot-step algorithms), so there are many variants of interior methods such as pro-
jective, affine, and path-following. Some interior-point methods inscribe an ellip-
soidal ball in the feasible region with its center at the current iterate, or first trans-
form the feasible space and then inscribe a hypersphere with the current iterate at
the center. Next an improving direction is found by joining the current iterate to
the point on the boundary of the ellipsoid or sphere that maximizes (or minimizes)
the linear objective function (obtained by solving a linear system of equations). A
point is then selected on the improving direction line as the next current iterate.
Sometimes this iterate is found along a line that is a linear combination of the
improving direction and some other direction.
During the period 1979–1996, there has been intensive interest in the develop-
ment of interior-point methods. These methods are related to classical least-square
methods used in numerical analysis for making fits to data or fitting simpler func-
tional forms to more complicated ones. Therefore interior research can tap a vast
literature of approximation theory. A theoretical breakthrough came in 1979; the
Russian mathematician L.G. Khachian discovered an ellipsoid algorithm whose run-
ning time in its worst case was significantly lower than that of the Simplex Algo-
rithm in its worst case—its iterates are not always required to be feasible. Other
theoretical results quickly followed, notably that of N. Karmarkar who discovered
an interior-point algorithm whose running time performance in its worst case was
113
114 INTERIOR-POINT METHODS
significantly lower than that of Kachiyan’s. This in turn was followed by more
theoretical results by others improving on the worst-case performance.
It soon became apparent that these worst-case running times were many, many
times greater than the actual running times of the Simplex Method on practical
problems and totally misleading as to what algorithm to choose for solving practical
problems. It appears that the problems encountered in practical applications belong
to a very specialized class (that is yet to be characterized in a simple way). Or to
put it in another way, the class of problems simply described by Ax = b, x ≥ 0,
c
T
x = min is much too broad, and the worst-case examples drawn from this very
broad class are far different from anything ever encountered in practical applications
or likely ever to be encountered in the future.
Attempts to characterize in a simple way the class (or classes) of practical prob-
lems from which one could derive a theoretical explanation of the excellent perfor-
mance times of various algorithms in practice have, in general, failed. In special
cases, such as the shortest path problem, the worst-case performance and perfor-
mance on actual problems are comparable, and the algorithms are very efficient.
There has been progress proving that average performance on classes of randomly
generated problems using a parametric variant of the Simplex Method resembles
that obtained on practical problems, but no one claims that these randomly gener-
ated problems belong to the class of practical problems.
Because the theoretical results are totally misleading as to what algorithm to
choose to solve a practical problem, a different approach is used. The linear pro-
gramming profession has accumulated a large collection of test problems drawn
from practical sources. These are used to compare running times of various pro-
posed algorithms. The general goal of these efforts is to develop algorithms that
surpass the performance of the Simplex Method on this collection of problems. For
example, Karmarkar claimed that on very large problems his technique is signifi-
cantly faster. As of this writing, as far as the authors of this book can ascertain,
there appears to be no algorithm that is a clear winner, i.e., solves all (or almost
all of) the test problems faster than all the other proposed methods. On problems
with many bounding hyperplanes in the neighborhood of the optimum point, an
interior-method will probably do better. On problems with relatively few boundary
planes (which is often the case in practice) an exterior method will be hard to beat.
For this reason, it is likely that the commercial software of the future will be some
sort of a hybrid, because one does not know which kind of problem is being solved or
because one wishes to obtain an extreme-point solution. Many specialized, efficient
codes have been proposed for solving structured linear programs such as network
problems, staircase problems, and block-angular systems.
Here we illustrate the primal affine method which is the same as Dikin’s method,
and which has the distinction of having been rediscovered by many.
4.1 BASIC CONCEPTS 115
q q
q
X
X
X
X
X
X
X
X
X
X
q
q
@
@
@q
s
:
¯x
t
s
:
ˆx
t
C
C
C
C
C
C
C
C
c
T
x = constant
Figure 4-1: Comparison of a Move from a Point ˆx
t
Near the Center Versus a Point
¯x
t
Near the Boundary.
4.1 BASIC CONCEPTS
In interior-point methods for solving linear programs, we are given, or must be
able to find, an initial starting interior feasible point from which we proceed to an
improved interior feasible point, and so on, eventually stopping by some criterion.
Our development here of an interior-point method will be for solving a linear
program in standard form:
Minimize c
T
x = z
subject to Ax = b, A : m × n,
x ≥ 0.
(4.1)
The rationale for the approach is based on the following observations. When
minimizing, one is first tempted to move from the current solution x
t
in the di-
rection of steepest descent of the objective function (i.e., in the negative gradient
of the objective function, which is the same as moving orthogonal to the objective
hyperplane c
T
x = constant). If the current iterate x
t
is an interior point, so that
x
t
> 0, such a move will in general violate the constraints Ax = b. To adjust for
this, one typically moves instead in the direction given by the projection of the
negative gradient of the objective onto the hyperplanes Ax = b. However, if x
t
is
close to the boundary hyperplanes, as x
t
=¯x
t
is in Figure 4-1, very little improve-
ment will occur. On the other hand, if the current iterate happens to be near the
“center,” such as x
t
=ˆx
t
in Figure 4-1, there could be a big improvement.
The key ideas behind interior-point methods are as follows, assuming an initial
feasible interior point is available and that all moves satisfy Ax = b:
1. Try to move through the interior in directions that show promise of moving
quickly to the optimal solution.
2. Recognize that if we move in a direction that sets the new point too “close” to
the boundary, this will be an obstacle that will impede our moving quickly to
an optimal solution. One way around this is to transform the feasible region
so that the current feasible interior point is at the center of the transformed
116 INTERIOR-POINT METHODS
feasible region. Once a movement has been made, the new interior point is
transformed back to the original space, and the whole process is repeated with
the new point as the center.
MAINTAINING FEASIBILITY
It is clear that the feasible point must satisfy Ax = b. Hence the only feasible
interior points x = x
t
are those that satisfy Ax
t
= b and x
t
> 0. As a result,
whenever we move from a point (in the original space or transformed space) x
t
to
x
t+1
= x
t
+ αp
t
, with α>0, we need to satisfy Ax
t+1
= b and x
t+1
> 0. This
implies that
Ax
t+1
= A(x
t
+ αp
t
)=b + αAp
t
. (4.2)
Since α = 0, the above equation implies that in order to stay feasible, p
t
must
satisfy
Ap
t
=0. (4.3)
In other words, p
t
= 0 must be orthogonal to the rows of A. A vector direction p
t
that satisfies Ap
t
= 0 is said to belong to the null space of A. It is easy to see that
the matrix
P = I − A
T
(AA
T
)
−1
A, (4.4)
called the projection matrix, will transform any vector v into Pv = p, and p will
be in the null space of A because AP = 0 is the zero matrix. Hence, in order to
maintain feasibility, every move from x
t
to x
t+1
must be in the null space of A, i.e.,
p
t
= Pv for some vector v.
STEEPEST DESCENT DIRECTION
The next thing to note is that the direction of maximum decrease or steepest descent
is the direction that is the negative of the gradient of the objective function z.
Hence, the direction of steepest descent is −c. However, the condition that we must
maintain feasibility implies that we must move in the direction of steepest descent
projected into the null space of A. That is,
p
t
= −Pc. (4.5)
Such a direction is called a projected gradient direction. To see that such a direction
is a direction of decrease, let x
t+1
= x
t
+ αp
t
, with α>0 and p
t
= −Pc. Then
c
T
x
t+1
= c
T
x
t
+ αp
t
= c
T
x
t
− αc
T
Pc
≤ c
T
x
t
, (4.6)
where the last inequality follows because P has the property that P = P
T
P , and
therefore
c
T
Pc = c
T
P
T
Pc =(Pc)
T
(Pc) ≥ 0.
Exercise 4.1 Prove P
T
P = P , where P is defined by (4.4).
4.2 PRIMAL AFFINE / DIKIN’S METHOD 117
STAYING IN THE INTERIOR
We plan to stay in the interior of x ≥ 0 by limiting the amount of movement
α>0 in the direction p
t
. The quantity α, called the steplength, is chosen to ensure
feasibility of the nonnegativity constraints. However, we do not want to move all the
way to the boundary because we want to avoid the need to move to the boundary
on subsequent iterations. Hence we move some portion of the way towards the
boundary. Based on the experience of others on practical problems, moving 50% to
90% of the way is a simple good rule for the primal affine/Dikin’s method that we
are about to describe.
CENTERING
What is left is the procedure for transforming the space so that an interior point x
t
is the “center” point in the transformed space. The basic idea is very simple—
we change the scaling of the solution so that each component is equidistant from
the constraint boundaries in the transformed feasible region; for example, after
rescaling, x
t
= e, where e =(1, 1,...,1)
T
. Let D = Diag (x
t
) be the diagonal
matrix having the components of the current iterate x
t
as the diagonal elements.
The simplest such scheme is to rescale x by setting x = D ˆx, so that ˆx
t
= e.
Substituting x = D ˆx in (4.1) we obtain
Minimize ˆc
T
ˆx = z
subject to
Aˆx = b,
A : m × n,
ˆx ≥ 0,
(4.7)
where ˆc = Dc and
A = AD. As a result of this transformation, clearly the projection
matrix for the system (4.7) is now defined by
P = I −
A
T
(
A
A
T
)
−1
A. (4.8)
Hence at each iteration t, we rescale x
t
to ˆx
t
= e and move to ˆx
t+1
in the
negative projected gradient direction −
P ˆc;thus
ˆx
t+1
= e − α
P ˆc, (4.9)
where e =(1, 1,...,1)
T
and
P and ˆc are defined as above. Finally, we rescale back
to the original space by setting x
t+1
= D ˆx
t+1
. Next we repeat the iterative process
replacing D by D = Diag (x
t+1
).
STOPPING RULE
The only item left out in our discussion is the stopping rule. The simple rule typi-
cally followed is to stop with an approximate optimal solution when the difference
between iterates x
t+1
and x
t
is “deemed” sufficiently small in the original space.
118 INTERIOR-POINT METHODS
4.2 PRIMAL AFFINE / DIKIN’S METHOD
Now we are prepared to write down the specific steps of the primal affine algorithm
and illustrate it with an example. The name affine method results from the following
definition:
Definition (Affine Transformation): Let F be an n ×n nonsingular matrix
and let d ∈
n
. The transformation y = Fx + d is called an affine transfor-
mation of x into y.
In the literature, an interior-point method that uses an affine transformation to
obtain the search direction is called an affine method.
The algorithm is described with
P computed explicitly. In practice
P is never
computed; instead, p = −
P ˆc is computed directly using a QR factorization or
Cholesky factorization.
Algorithm 4.1 (Primal Affine/Dikin’s Method) Assume that an initial feasible
interior point solution x = x
o
> 0 satisfying Ax
o
= b is known. The steps of the algorithm
are as follows:
1. Initialize Counter: Set t ← 0.
2. Create D: Set D = Diag (x
t
).
3. Compute Centering Transformation: Compute
A = AD,ˆc = Dc.
4. Determine the Projection Matrix: Compute
P = I −
A
T
(
A
A
T
)
−1
A.
5. Compute Steepest Descent Direction: Set p
t
= −
P ˆc.
6. Set θ = −min
j
p
t
j
.
7. Test for Unbounded Objective.Ifθ ≤ 0.0 report the objective as being unbounded
and stop.
8. Obtain ˆx
t+1
: Compute
ˆx
t+1
= e +
α
θ
p
t
,
where e =(1, 1,...,1)
T
and α is set to a number strictly between 0 and 1. Typically,
α is set to be between 0.5 and 0.95 for primal affine methods (in most implementa-
tions it is set between 0.9 and 0.95).
9. Transform Back to Original Space: Compute x
t+1
= D ˆx
t+1
.
10. Termination Check: If x
t+1
is “close” to x
t
, report x
t+1
as “close” to optimal and
stop.
11. Set t ← t + 1 and go to Step 2.
Example 4.1 (Illustration of the Primal Affine Method) Consider the two-variable
example
Minimize −2x
1
− x
2
= z
subject to x
1
+ x
2
≤ 5
2x
1
+3x
2
≤ 12
and x
1
≥ 0,x
2
≥ 0.
4.2 PRIMAL AFFINE / DIKIN’S METHOD 119
An initial starting feasible interior solution x
o
1
=1,x
o
2
= 2 is given.
We start by transforming this into the standard form by adding slack variables x
3
and
x
4
to obtain
Minimize −2x
1
− x
2
= z
subject to x
1
+ x
2
+ x
3
=5
2x
1
+3x
2
+ x
4
=12
and x
j
≥ 0 for j =1,...,4.
The corresponding feasible interior-point starting solution is
x
o
=(1224)
T
. (4.10)
The objective value at x
o
is z
o
= −4.
We start the first iteration by setting the scaling matrix D to be
D =
1
2
2
4
. (4.11)
The rescaled
A matrix and objective function coefficient ˆc are computed as
A = AD =
1110
2301
1
2
2
4
=
1220
2604
, (4.12)
ˆc = Dc =
1
2
2
4
−2
−1
0
0
=
−2
−2
0
0
. (4.13)
Next we compute the projection matrix
P :
P = I −
A
T
A
A
T
−1
A
=
1
1
1
1
−
12
26
20
04
1220
2604
12
26
20
04
−1
1220
2604
=
.8831 −.2597 −.1818 −.0519
−.2597 .3117 −.1818 −.3377
−.1818 −.1818 .2727 .3636
−.0519 −.3777 .3636 .5325
. (4.14)
Hence the projected gradient is
p
o
= −
P ˆc = −
.8831 −.2597 −.1818 −.0519
−.2597 .3117 −.1818 −.3377
−.1818 −.1818 .2727 .3636
−.0519 −.3777 .3636 .5325
−2
−2
0
0
=
1.2468
0.1039
−0.7273
−0.7792
. (4.15)
Next we pick θ by
θ = −min
j
p
o
j
= .7792. (4.16)
120 INTERIOR-POINT METHODS
Next we rescale the current iterate to ˆx
o
= D
−1
x
o
= e and move to ˆx
1
in the transformed
space:
ˆx
1
=
1
1
1
1
+
α
θ
p
o
=
1
1
1
1
+
0.5
.7792
1.2468
0.1039
−0.7273
−0.7792
=
1.8000
1.0667
0.5333
0.5000
. (4.17)
Transforming this point to the original space we obtain
x
1
= D ˆx
1
=
1
2
2
4
1.8000
1.0667
0.5333
0.5000
=
1.8000
2.1333
1.0667
2.0000
. (4.18)
The objective value at x
1
is z
1
= −5.7333. This completes the first iteration of the method.
Because x
1
is quite different from x
o
we proceed to the next iteration.
We update the scaling matrix D to Diag (x
t+1
):
D =
1.8000
2.1333
1.0667
2.0000
. (4.19)
The rescaled
A matrix and objective function coefficient ˆc are computed as
A = AD =
1110
2301
1.8000
2.1333
1.0667
2.0000
=
1.8000 2.1333 1.0667 0.0000
3.6000 6.4000 0.0000 2.0000
, (4.20)
ˆc = Dc =
1.8000
2.1333
1.0667
2.0000
−2
−1
0
0
=
−3.6000
−2.1333
0.0000
0.0000
. (4.21)
Next we need to compute the projection matrix
P . This is
P = I −
A
T
A
A
T
−1
A
=
.6203 −.3718 −.3031 .0733
−.3718 .2885 .0505 −.2539
−.3031 .0505 .4106 .3841
.0733 −.2539 .3841 .6806
. (4.22)
Hence the projected gradient is
p
1
=
P ˆc =
.6203 −.3718 −.3031 .0733
−.3718 .2885 .0505 −.2539
−.3031 .0505 .4106 .3841
.0733 −.2539 .3841 .6806
−3.6000
−2.1333
0.0000
0.0000
=
1.4399
−0.7231
−0.9836
−0.2779
.
(4.23)
4.3 INITIAL SOLUTION 121
Next we pick θ by
θ = −min
j
p
1
j
=0.9836 (4.24)
Next we rescale the current iterate to ˆx
1
= D
−1
x
1
= e and move to ˆx
2
in the transformed
space:
ˆx
2
=
1
1
1
1
+
α
θ
p =
1
1
1
1
+
0.5
0.9836
1.4399
−0.7231
−0.9836
−0.2779
=
1.7320
0.6324
0.5000
1.8588
. (4.25)
Transforming this point to the original space, we obtain
x
2
= D ˆx
2
=
1.8000
2.1333
1.0667
2.0000
1.7320
0.6324
0.5000
1.8588
=
3.1175
1.3492
0.5333
1.7175
. (4.26)
The objective value is z
2
= −7.5842. This completes the second iteration of the method.
Because x
2
is quite different from x
1
we proceed to the next iteration. We leave the rest
of the iterations as an exercise for the reader.
It is clear that the solution of the above problem using the primal affine method will
take many more iterations than the Simplex Method. However, the experience with very
large practical problems indicates that interior-point methods tend to solve faster than the
Simplex Method.
Exercise 4.2 Complete the steps of the algorithm on Example 4.1 using the Primal
Affine / Dikin software option.
4.3 INITIAL SOLUTION
So far we have assumed that an initial feasible interior-point solution is available.
If an initial feasible interior solution is not available, we can easily generate one by
picking an arbitrary x
o
> 0, say x
o
= e where e =(1, 1,...,1)
T
, and setting up the
following linear program with one artificial variable x
a
≥ 0 and M (called “Big-M ”
in the literature) a large scalar constant:
Minimize c
T
x + Mx
a
= z
subject to Ax +(b − Ax
0
)x
a
= bx≥ 0,x
a
≥ 0.
(4.27)
Then x = x
0
> 0 and x
a
= 1 is clearly a feasible solution to (4.27).
Example 4.2 (Initial Feasible Point) Consider the two-variable example described
in Example 4.2 with the slacks put in:
Minimize −2x
1
− x
2
= z
subject to x
1
+ x
2
+ x
3
=5
2x
1
+3x
2
+ x
4
=12
and x
j
≥ 0 for j =1,...,4.