Birge J.R., Louveaux F. Introduction to Stochastic Programming
Подождите немного. Документ загружается.
224 5 Two-Stage Recourse Problems
L =
⎛
⎝
W
1l,J
1
00
0 ...W
kl,J
k
... 0
00W
Kl,J
K
⎞
⎠
.
Because D has rank at least m
1
, each W
kl,J
k
in L has rank at most m
2
.Thisgives
the result.
For Example 2 from Section 5.1, the solution, x
1
= 1andy
1
k
=
ξ
k
−1, k =
1,2, 3 , corresponds to the basis:
B
1
=
⎛
⎜
⎜
⎝
1100
1000
1010
1001
⎞
⎟
⎟
⎠
, (5.4)
which already has the form in Proposition 5 with D =
11
10
, C =
00
00
, F =
10
10
,and L =
10
01
. Note in this case that W
1
= W
1u
= W
1l
=[],anempty
matrix.
To show how the partition in Proposition 5 is used, consider the forward trans-
formation to find the basic values of (x
I
0
,x
I
1
,...,x
I
K
,y
J
1
,...,y
J
k
) , which we write
as (x
B
,y
B
) , that solve:
Dx
B
+Cy
B
= b
; Fx
B
+ Ly
B
= h
, (5.5)
where b
=
b
h
u
, h
= h
l
, h
u
corresponds to the components of the right-hand
side for rows of T in D ,and h
l
corresponds to the components with rows in F .
Note that L is invertible; so,
y
B
= L
−1
(h
−Fx
B
) . (5.6)
Substituting in the first system of equations yields
(D −CL
−1
F)x
B
= b
−CL
−1
h
. (5.7)
Hence, we use L to solve for the columns of L
−1
F and L
−1
h
, then form the work-
ing basis, (D−CL
−1
F) ,tosolvefor x
B
, and multiply x
B
again by L
−1
F and sub-
tract from L
−1
h
to obtain y
B
. Because most of the work involves just the square
block matrices in L and the working basis, substantial effort can be saved in the de-
composition procedure (see Exercise 1). The backward transformation can also be
performed by taking advantage of this structure (see Exercise 2). The other forward
transformation in the simplex method to find the leaving column is, of course, the
same as the operations used in (5.6)and(5.7).
5.5 Basis Factorization and Interior Point Methods 225
For basis B
1
in the Example 2, b
=[1,0, 1]
T
, D =
11
10
,and C = 0
2×2
,
yields a solution to 5.7 with x
B
=[1,9]
T
, where again the second component cor-
responds to the first-period slack variable. Now, with h
=[2,4] , F =
10
10
,and
L = I , the solution to (5.6)is y
B
=[1,3]
T
.
The entire simplex method then has the following form.
Basic Factorization Simplex Method
Step 0. Suppose that (x
0
B
0
,y
0
B
0
)=(x
0
I
0
0
,...,x
0
I
0
K
,y
0
J
0
0
,...,y
0
J
0
K
) is an initial basic feasible
solution for (1.1), with initial indices partitioned according to B
0
= {
β
0
1
,...,
β
0
l
0
}=
{I
0
i
,i = 0,...,K} and B
0
= {
β
0,
1,1
,...,
β
0,
1,l
1
,...,
β
0,
K,1
,...,
β
0,
K,l
K
}= J
0
j
, j = 1,...,K .
Let the initial permutation matrix be P
0
,andset
ν
= 0.
Step 1.Solve (
ρ
T
,
π
T
)
DC
FL
=(c
T
B
0
, ˆq
T
β
0
) ,where ˆq
k,i
= p
k
q
k,i
.
Step 2.Find ¯c
s
= min
j
{c
j
−(
ρ
T
|
π
T
)P
ν
(A
T
·j
| T
T
1,·j
|··· |T
T
K,·j
)
T
} and ¯q
k
,s
=
min
j,k
{p
k
q
k, j
−(
ρ
T
|
π
T
)P
ν
(0...W
k,·j
...0)}.If ¯c
s
≥ 0and¯q
k
,s
≥ 0 , then stop;
the current solution is optimal. Otherwise, if ¯c
s
< ¯q
k
,s
,gotoStep4.If ¯c
s
≥ ¯q
k
,s
,
go to Step 3.
Step 3. Solve for the entering column,
DC
FL
¯
W
k
,·s
= P
ν
(0...W
T
k
,·s
...0)
T
.Let
θ
= x
ν
B
ν
(r)
/
¯
W
k
,rs
= min
¯
W
k
,is
>0 , 1≤i≤l
ν
{x
ν
B
ν
(i)
/
¯
W
k
,is
} (5.8)
and
θ
= y
ν
B
ν
(r
)
/
¯
W
k
,r
s
= min
¯
W
k
,is
>0 , l
ν
+1≤i≤m
1
+Km
2
{y
ν
B
ν
(i)
/
¯
W
k
,is
} . (5.9)
If no minimum exists in either (5.8)or(5.9), then stop; the problem is unbounded.
Otherwise, if
θ
<
θ
,gotoStep5.If
θ
≥
θ
,gotoStep6.
Step 4. Solve for the entering column,
DC
FL
¯
A
·s
= P
ν
(A
T
·s
| T
T
1,·s
| ... | T
T
K,·s
)
T
.
Let
θ
= x
ν
B
ν
(r)
/
¯
A
rs
= min
¯
A
is
>0 , 1≤i≤l
ν
{x
ν
B
ν
(i)
/
¯
A
is
} (5.10)
and
θ
= y
ν
B
ν
(r
)
/
¯
A
r
s
= min
¯
A
is
>0 , l
ν
+1≤i≤m
1
+Km
2
{y
ν
B
ν
(i)
/
¯
A
is
} . (5.11)
If no minimum exists in either (5.10)or(5.11), then stop; the problem is unbounded.
Otherwise, if
θ
<
θ
,gotoStep5.If
θ
≥
θ
,gotoStep6.
226 5 Two-Stage Recourse Problems
Step 5.Let B
ν
+1
= B
ν
, B
ν
+1
= B
ν
, I
ν
+1
i
= I
ν
i
,and J
ν
+1
= J
ν
. Suppose B
ν
(r)=
I
ν
j,w
= t .If x
s
is entering, then let B
ν
+1
(r)=I
ν
+1
( j,w)=s .If y
k
s
is entering,
then let B
ν
+1
(i)=B
ν
(i+ 1) , i ≥ r , I
ν
+1
j,i
= I
ν
j,i+1
, i ≥w , J
ν
+1
k
,l
k
+1
= s
,and l
k
=
l
k
+ 1 . Update P
ν
to P
ν
+1
, the factorization correspondingly, let
ν
=
ν
+ 1,and
go to Step 1.
Step 6.Let B
ν
+1
= B
ν
, B
ν
+1
= B
ν
, I
ν
+1
i
= I
ν
i
,and J
ν
+1
= J
ν
. Suppose
B
ν
(r
)=J
ν
k,w
= t .If x
s
is entering, then let B
ν
+1
(
∑
k
j= 1
l
j
)=I
ν
+1
(k,l
k
+ 1)=s ,
B
ν
+1
(i)=B
ν
(i −1),i >
∑
k
j= 1
l
j
, l
k
= l
k+1
, J
ν
+1
k,i
= J
ν
k,i+1
, i ≥ w .If y
k
s
is en-
tering, then let B
ν
+1
(i)=B
ν
(i + 1) , i ≥ r , I
ν
+1
j,i
= I
ν
j,i+1
, i ≥ w , J
ν
+1
k
,l
k
+1
= s
,
J
ν
+1
k,i
= J
ν
+1
k,i+1
, i ≥ w , l
k
= l
k
−1, and l
k
= l
k
+ 1 . Update P
ν
to P
ν
+1
,the
factorization correspondingly, let
ν
=
ν
+ 1,andgotoStep1.
For updating a factorization of the basis as used in (5.6)and(5.7), several cases
need to be considered according to the possibilities in Steps 5 and 6 (see Exercise 3).
If the entering and leaving variables are both in x , then only D changes. Substantial
effort can again be saved. In other cases, only one block of L is altered by any
iteration so we can again achieve some savings by only updating the corresponding
parts of L
−1
F and L
−1
h .
As mentioned earlier, this procedure can also apply to the dual of (1.1)andthe
primal. In this case, the procedure can mimic decomposition procedures and entails
essentially the same work per iteration as the L -shaped method (see Birge [1988b])
or the inner linearization method applied to the dual. If choices of entering columns
are restricted in a special variant of a decomposition procedure, then factorization
and decomposition follow the same path.
In general, decomposition methods have been favored for this class of problems
because they offer other paths of solutions, require less overhead, and, by maintain-
ing separate subproblems, allow for parallel computation. The extensive form offers
little hope for efficient solution, so it is not surprising that even sophisticated fac-
torizations would not prove beneficial. Because most commercial methods already
have substantial capabilities for exploiting general matrix structure, it is difficult
to see how substantial gains could be obtained from basis factorization alone for
a direct extreme point approach. Combinations of decomposition and factorization
approaches may, however, be beneficial, as observed in Birge [1985b].
Factorization schemes also offer substantial promise for interior point methods,
where there is much speculation that the solution effort grows linearly in the size
of the problem. This observation is supported by the results we present here. For
this discussion, we assume that the interior point method follows a standard form
version of Karmarkar’s projective algorithm (Karmarkar [1984]).We also assume an
unknown optimal objective value and use Todd and Burrell’s [1986] method for up-
dating a lower bound on the optimal objective value. We use an initial lower bound,
as is often available in practice. An alternative is Anstreicher’s [1989] method to
obtain an initial lower bound.
5.5 Basis Factorization and Interior Point Methods 227
Many other interior point methods are available (see, e.g., Roos, Terlaky, and Vial
[2005] and Ye [1997]). In particular, many commercial solvers use the homogeneous
self-dual formulation of the standard linear program (see, e.g., Andersen [2009]).
Other interior point methods and interpretations include path-following, logarithmic
barrier, and affine scaling (see Roos, Terlaky, and Vial [2005] for descriptions of
alternatives). They all follow similar steps to the method given below.
We first describe the algorithm for a standard linear program:
min c
T
x
s. t. Ax = b ,
x ≥0 ,
(5.12)
where x ∈ ℜ
n
, c ∈ Z
n
(i.e., an n -vector of rationals), b ∈ Z
m
, A ∈ Z
m×n
with
optimal value c
T
x
∗
= z
∗
. In referring to the parameters in (5.12), we use ext as a
modifier, e.g., c
ext
, when necessary to distinguish the parameters in (5.12) from our
standard stochastic program form in (1.1).
Suppose we have a strictly interior feasible point x
0
of (5.12), i.e.,
Ax
0
= b , x
0
> 0 , (5.13)
a lower bound
β
0
on z
∗
, and the set of optimal solutions in (5.12) is bounded. Note
that if we do not have a feasible solution, we can solve a phase-one problem or use
the self-dual form of the problem. In that case, the goal becomes finding (x,t,
λ
) to
solve:
min 0
s. t. Ax −bt = 0 ,
−A
T
λ
+ct ≥ 0,
b
T
λ
−c
T
x ≥ 0,
x ≥0, t ≥ 0 ,
(5.14)
which can be solved by an interior point method initiated at any solution (x
0
,t
0
,
λ
0
)
with x
0
> 0andt
0
> 0 by iteratively choosing search directions to reduce the
infeasibility of the system (5.14) with solution (x
k
,t
k
,
λ
k
) at iteration k .
For exposition here, we follow the standard form variant of the projective scaling
algorithm, which creates a sequence of points x
0
, x
1
, ..., x
k
by the following
steps.
Standard Form Projective Scaling Method
Step 0.Set
ν
= 0 and lower bound
β
0
≤ z
∗
.
Step 1.If c
T
x
ν
−
β
ν
is small enough, i.e., less than a given positive number
ε
,then
stop. Otherwise, go to Step 2.
228 5 Two-Stage Recourse Problems
Step 2.Let D = diag{x
ν
1
,...,x
ν
n
},
ˆ
A :=[AD,−b] ,andlet
Π
ˆ
A
be the projection
onto the null space of
ˆ
A .Find
u =
Π
ˆ
A
Dc
0
, v =
Π
ˆ
A
0
1
, (5.15)
and let
μ
(
β
ν
)=min{u
i
−
β
ν
v
i
: i = 1,...,n + 1}.If
μ
(
β
ν
) ≤ 0,let
β
ν
+1
=
β
ν
.
Otherwise, let
β
ν
+1
= min{u
i
/v
i
: v
i
> 0,i = 1,...,n + 1}.GotoStep3.
Step 3.Let c
p
= u −
β
ν
+1
v −(c
T
x
ν
−
β
ν
+1
)e/(n + 1) ,where e =(1,...,1)
T
∈
ℜ
n+1
.Let
g
=
1
n + 1
e −
α
c
p
c
p
2
.
Let
g ∈ℜ
n
consist of the first n components of g
.Then x
ν
+1
= Dg/g
n+1
,
ν
=
ν
+ 1 , go to Step 1.
For the purpose of obtaining a worst-case bound, the step length
α
in the def-
inition of g
may be set equal to
1
3(n+1)
, (Gay [1987]), but better performance is
obtained by choosing
α
using a line search.
To show how the structure of a stochastic program can be exploited in these
methods, we consider the number of arithmetic operations in a complexity analysis.
The main computational effort in each iteration of the algorithm is to compute the
projections in (5.15), which requires, for n > m , O(n
3
) arithmetic operations (and,
on average, O(n
2.5
) , operations per iteration using a rank-one updating scheme).
In O(n/
ε
) iterations, or with some modifications in O(
√
n/
ε
) , the method finds a
solution with O(
ε
) precision.
In our case, if we consider the stochastic program (1.1) in the extensive form
(5.12), then n = n
2
+ Kn
2
, m = m
1
+ Km
2
,and x
ext
=
⎛
⎜
⎜
⎜
⎝
x
y
1
.
.
.
y
K
⎞
⎟
⎟
⎟
⎠
, c
ext
=
⎛
⎜
⎜
⎜
⎝
c
p
1
q
1
.
.
.
p
K
q
K
⎞
⎟
⎟
⎟
⎠
,
b
ext
=
⎛
⎜
⎜
⎜
⎝
b
h
1
.
.
.
h
K
⎞
⎟
⎟
⎟
⎠
,and
A
ext
=
⎛
⎜
⎜
⎜
⎝
A 0 ... 0
T
1
W ... 0
.
.
.0
.
.
.
0
T
K
0 ... W
⎞
⎟
⎟
⎟
⎠
. (5.16)
The main computational work at each step of the projective scaling algorithm is to
compute the projection in (5.15), which can be written as
Π
ˆ
A
=(I −
ˆ
A
T
(
ˆ
A
ˆ
A
T
)
−1
ˆ
A) , (5.17)
5.5 Basis Factorization and Interior Point Methods 229
where (
ˆ
A
ˆ
A
T
)=AD
2
A
T
+ bb
T
:= M + bb
T
. In this case, the work is dominated by
computing M
−1
(or solving systems with coefficient matrix, M = AD
2
A
T
). For
the general A in the formulation in (5.12), using the specific A
ext
in the stochastic
program extensive form as in (5.16) and letting D
0
= diag(x
ν
) , D
k
= diag(y
ν
k
) ,
k = 1,...,K , we would have
M =
⎛
⎜
⎜
⎜
⎝
AD
2
0
A
T
AD
2
0
T
T
1
... AD
2
0
T
T
K
T
1
D
2
0
A
T
T
1
D
2
0
T
T
1
+WD
2
1
W
T
... T
1
D
2
0
T
T
K
.
.
.
.
.
.
.
.
.
.
.
.
T
K
D
2
0
A
T
T
1
D
2
0
T
T
K
... T
K
D
2
0
T
T
K
+WD
2
K
W
T
⎞
⎟
⎟
⎟
⎠
, (5.18)
which is clearly much denser than the original constraint matrix in (1.1). In this case,
a straightforward implementation of an interior point method that solves systems
with M is quite inefficient.
To see the structure, we consider Example 2 from Section 5.1. Here,
A
ext
=
⎛
⎜
⎜
⎝
11000000
101−10 0 0 0
100 0 1−10 0
1000001−1
⎞
⎟
⎟
⎠
. (5.19)
Now, let x
0
ext
=(3, 7,1,3, 1,2,2, 1)
T
in (5.13)represent x
0
=(3, 7)
T
, y
0
1
=
(1,3)
T
, y
0
2
=(1,2)
T
,and y
0
2
=(2,1)
T
in Example 2. We then have D
0
=
diag(3,7) , D
1
= diag(1, 3) , D
2
= diag(1, 2) ,and D
3
= diag(2, 1) . The matrix
M in this case is:
M =
⎛
⎜
⎜
⎝
58999
9199 9
99149
99914
⎞
⎟
⎟
⎠
. (5.20)
While M is dense, it in fact has a great deal of structure that can be exploited in
any solution scheme. This is the object of the factorization scheme given by Birge
and Qi [1988] (see also Birge and Holmes [1992] for an implementation discussion).
The following proposition gives the essential characterization of that factorization.
Proposition 17. Let S
0
= I
2
∈ℜ
m
1
×m
1
,S
l
= W
l
D
2
l
W
T
l
,l= 1,...,K,
S = diag{S
0
,...,S
K
}.Then S
−1
= diag{S
0
,S
−1
1
,...,S
−1
N
}.Let I
1
and I
2
be iden-
tity matrices of dimensions n
1
and m
1
, respectively. Let
G
1
=(D
0
)
−2
+ A
T
S
−1
0
A +
K
∑
l=1
T
T
l
S
−1
l
T
l
, G
2
= −AG
−1
1
A
T
, (5.21)
U =
⎛
⎜
⎜
⎜
⎝
AI
2
T
1
0
.
.
.
.
.
.
T
K
0
⎞
⎟
⎟
⎟
⎠
, V =
⎛
⎜
⎜
⎜
⎝
A −I
2
T
1
0
.
.
.
.
.
.
T
K
0
⎞
⎟
⎟
⎟
⎠
.
230 5 Two-Stage Recourse Problems
If A , W
k
,k = 1,...,K have full row rank, then G
2
and M are invertible and
M
−1
= S
−1
−S
−1
U
I
1
−G
−1
1
A
T
0 I
2
I
1
0
0 −G
−1
2
I
1
0
AI
2
G
−1
1
0
0 I
2
V
T
S
−1
. (5.22)
Proof:Exercise6.
Following the assumptions, the number of arithmetic operations using this fac-
torization can be reduced from O((n
1
+ Kn
2
)
4
) as in the general projective scal-
ing method. Using the factorization, the effort is, in fact, dominated by O(K(n
3
2
+
n
2
2
n
1
+n
2
n
2
1
)) . It is also possible to reduce this bound further with a partial rank-one
updating scheme as mentioned earlier. In this case, for n = n
1
+ Kn
2
, the complex-
ity using the factorization in (5.22) becomes O((n
0.5
n
2
2
+ nmax{n
1
,n
2
}+ n
3
1
)n/
ε
)
for the entire algorithm, or, if K ∼ n
1
∼ n
2
, the full arithmetic complexity is
O(n
2.5
/
ε
) , compared to the general result of O(n
3.5
/
ε
) . Thus, the factorization in
(5.22) provides an order of magnitude improvement over a general solution scheme
if the number of realizations K approaches the number of variables in the first and
second stage.
In practice, we would not compute M
−1
explicitly, but solve a set of systems as
follows:
Mv = u (5.23)
using
v = p −r , (5.24)
where
Sp= u , Gq = V
T
p , Sr = Uq , (5.25)
where G is the inverse of the matrix between U and V
T
in (5.22):
G =
G
1
A
T
−A 0
=
G
1
0
0 I
2
I
1
0
AI
2
−1
I
1
0
0 −G
2
I
1
−G
−1
1
A
T
0 I
2
−1
. (5.26)
The systems in (5.25) require solving systems with S
l
, computation of G
1
and
G
2
, and solving systems with G
1
and G
2
. In practice, we find a Cholesky factor-
ization of each S
l
, use them to find G
1
and G
2
, and find Cholesky factorizations
of G
1
and G
2
.
For Example 2, with initial values (x
0
,y
0
1
,...,y
0
K
) given above, we have
S
0
=[1];S
1
=[10];S
2
=[5];S
3
=[5]; (5.27)
G
1
=
1
9
0
0
1
49
+
11
11
+
1
10
0
00
+
1
5
0
00
+
1
5
0
00
=
1.61 1
11.02
; (5.28)
5.5 Basis Factorization and Interior Point Methods 231
G2 = −[11]
1.61 1
11.02
−1
1
1
=[−0.98]; (5.29)
U =
⎛
⎜
⎜
⎝
111
100
100
100
⎞
⎟
⎟
⎠
;V =
⎛
⎜
⎜
⎝
11−1
10 0
10 0
10 0
⎞
⎟
⎟
⎠
. (5.30)
To solve for v in Mv = u ,wefirstsolve Sp= u as
⎛
⎜
⎜
⎝
1000
01000
0050
0005
⎞
⎟
⎟
⎠
⎛
⎜
⎜
⎝
p
1
p
2
p
3
p
4
⎞
⎟
⎟
⎠
=
⎛
⎜
⎜
⎝
u
1
u
2
u
3
u
4
⎞
⎟
⎟
⎠
(5.31)
to obtain:
p
1
= u
1
, p
2
= 0.1u
2
; p
3
= 0.2u
3
; p
4
= 0.2u
4
. (5.32)
Next, we find
V
T
p =
⎛
⎝
u
1
+ 0.1u
2
+ 0.2u
3
+ 0.2u
4
u
1
−u
1
⎞
⎠
. (5.33)
Next, we solve Gq = V
T
p as follows:
• find q
1
such that
G
1
0
0 I
2
q
1
= V
T
p as
⎛
⎝
1.61 1 0
11.02 0
001
⎞
⎠
q
1
=
⎛
⎝
u
1
+ 0.1u
2
+ 0.2u
3
+ 0.2u
4
u
1
−u
1
⎞
⎠
(5.34)
to obtain q
1
=
⎛
⎝
0.03u
1
+ 0.16u
2
+ 0.32u
3
+ 0.32u
4
0.95u
1
−0.16u
2
−0.31u
3
−0.31u
4
−u
1
⎞
⎠
;
• find q
2
such that q
2
=
I
1
0
AI
2
q
1
as
q
2
=
⎛
⎝
100
010
111
⎞
⎠
q
1
=
⎛
⎝
0.03u
1
+ 0.16u
2
+ 0.32u
3
+ 0.32u
4
0.95u
1
−0.16u
2
−0.31u
3
−0.31u
4
−0.02u
1
+ 0.003u
2
+ 0.006u
3
+ 0.001u
4
⎞
⎠
; (5.35)
• find q
3
such that
I
1
0
0 −G
2
q
3
= q
2
as
⎛
⎝
10 0
01 0
000.98
⎞
⎠
q
3
= q
2
(5.36)
232 5 Two-Stage Recourse Problems
to obtain q
3
=
⎛
⎝
0.03u
1
+ 0.16u
2
+ 0.32u
3
+ 0.32u
4
0.95u
1
−0.16u
2
−0.31u
3
−0.31u
4
−0.02u
1
+ 0.003u
2
+ 0.01u
3
+ 0.01u
4
⎞
⎠
;
• find q = q
4
such that q
4
=
I
1
−G
−1
1
A
T
0 I
2
q
3
as
q =
⎛
⎝
10−0.03
01−0.95
00 1
⎞
⎠
q
3
=
⎛
⎝
0.03u
1
+ 0.16u
2
+ 0.32u
3
+ 0.32u
4
0.97u
1
−0.16u
2
−0.32u
3
−0.32u
4
−0.02u
1
+ 0.003u
2
+ 0.01u
3
+ 0.01u
4
⎞
⎠
. (5.37)
The next step is to solve for Sr = Uq as
⎛
⎜
⎜
⎝
1000
01000
0050
0005
⎞
⎟
⎟
⎠
⎛
⎜
⎜
⎝
r
1
r
2
r
3
r
4
⎞
⎟
⎟
⎠
=
⎛
⎜
⎜
⎝
0.98u
1
+ 0.003u
2
+ 0.01u
3
+ 0.01u
4
0.03u
1
+ 0.16u
2
+ 0.32u
3
+ 0.32u
4
0.03u
1
+ 0.16u
2
+ 0.32u
3
+ 0.32u
4
0.03u
1
+ 0.16u
2
+ 0.32u
3
+ 0.32u
4
⎞
⎟
⎟
⎠
(5.38)
or r =
⎛
⎜
⎜
⎝
0.98u
1
+ 0.003u
2
+ 0.01u
3
+ 0.01u
4
0.003u
1
+ 0.02u
2
+ 0.03u
3
+ 0.03u
4
0.01u
1
+ 0.03u
2
+ 0.06u
3
+ 0.06u
4
0.01u
1
+ 0.03u
2
+ 0.06u
3
+ 0.06u
4
⎞
⎟
⎟
⎠
, which finally yields v = p −r as
v =
⎛
⎜
⎜
⎝
0.02u
1
−0.003u
2
−0.01u
3
−0.01u
4
−0.003u
1
+ 0.08u
2
−0.03u
3
−0.03u
4
−0.01u
1
−0.03u
2
+ 0.14u
3
−0.06u
4
−0.01u
1
−0.03u
2
−0.06u
3
−0.14u
4
⎞
⎟
⎟
⎠
, (5.39)
which can be seen as M
−1
u for M in (5.20).
Now, for the projection operation defined in (5.17), note that
(
ˆ
A
ˆ
A
T
)
−1
ˆ
A = M
−1
ˆ
A−M
−1
bb
T
M
−1
ˆ
A/(1 + b
T
M
−1
b), (5.40)
which requires finding V
1
and v
2
such that MV
1
=
ˆ
A and Mv
2
= b where
ˆ
A =
⎛
⎜
⎜
⎝
37000000−10
301−30000 −1
30001−20 0 −2
3000002−1 −4
⎞
⎟
⎟
⎠
and b =
⎛
⎜
⎜
⎝
10
1
2
4
⎞
⎟
⎟
⎠
, (5.41)
where note that v
2
is also the negative of the last column of V
1
.
Using (5.39) then yields
V
1
=[V
11
−v
2
]=
⎛
⎜
⎜
⎝
0.01 0.14 −0.003 0.01 −0.01 0.01 −0.01 0.01 −0.16
0.05 −0.02 0.08 −0.25 −0.03 0.06 −0.06 0.00.14
0.11 −0.05 −0.03 0.10 0.14 −0.27 −0.13 0.06 0.08
0.11 −0.05 −0.03 0.10 −0
.06 0.13 0.27 −0.14 −0.32
⎞
⎟
⎟
⎠
(5.42)
5.5 Basis Factorization and Interior Point Methods 233
From (5.40), (
ˆ
A
ˆ
A
T
)
−1
ˆ
A = V
1
−v
2
b
T
V
1
/(1 + b
T
v
2
) or
(
ˆ
A
ˆ
A
T
)
−1
ˆ
A =
⎛
⎜
⎜
⎝
−0.02 0.10 0.003 −0.01 −0.003 0.01 −0.04 0.02 −0.04
0.08 0.02 0.08 −0.24 −0.03 0.07 −0.04 0.02 0.04
0.12 −0.02 −0.03 0.10 0.14 −0.27 −0.11 0.06 0.02
0.03 −0.14 −0.02 0.06 −
0.06 0.11 0.21 −0.11 −0.09
⎞
⎟
⎟
⎠
.
(5.43)
Finally, the search direction components u and v in (5.15)arethen:
u=(I −
ˆ
A(
ˆ
A
ˆ
A
T
)
−1
ˆ
A)
Dc
ext
0
=
⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎝
0.31
0.17
0.53
0.42
0.43
0.47
0.24
0.55
0.22
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎠
;v =(I −
ˆ
A(
ˆ
A
ˆ
A
T
)
−1
ˆ
A)
⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎝
0
0
0
0
0
0
0
0
1
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎠
=
⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎝
0.22
0.32
−0.04
0.12
−0.02
0.04
0.18
−0.09
0.28
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎠
.
(5.44)
Note how these operations only required solutions with G
1
( n
1
×n
1
), G
2
( m
1
×
m
1
), and S ( K solutions using m
2
×m
2
matrices). After finding u and v ,the
other operations in the project scaling method only involve simple operations on
vectors of the same size. Exercise 7 asks for completion of these operations until
the objective value is within 0.01 of the bound.
Other factorizations or formulations can also yield advantages for interior point
methods. These include the following approaches:
1. Schur complement updates;
2. Column splitting;
3. Solution of the dual.
The Schur complement approach is used in many interior point method implemen-
tations. The basic idea is to write M as the sum of a matrix with sparse columns,
A
s
D
2
s
A
T
s
, and a matrix with dense columns, A
d
D
2
d
A
T
d
. Using a Cholesky factoriza-
tion of the sparse matrix, LL
T
= A
s
D
2
s
A
T
s
, the method involves solving Mu = v by:
LL
T
−A
d
D
d
D
d
A
T
d
I
v
w
=
u
0
, (5.45)
which requires solving [I+D
d
A
T
d
(LL
T
)
−1
A
d
D
d
]w = −D
d
A
T
d
(LL
T
)
−1
b and LL
T
v =
b + A
d
D
d
w ,where I + D
d
A
T
d
(LL
T
)
−1
A
d
D
d
is a Schur complement.
The Schur complement is thus quite similar to the factorization method given ear-
lier. If every column of x is considered a dense column, then the remaining matrix
is quite sparse but rank deficient. The factorization in (5.22) is a method for main-
taining an invertible matrix when A
s
D
2
s
A
T
s
is singular. It can thus be viewed as an
extension of the Schur complement to the stochastic linear program. Because of the
possible rank deficiency and the size of the Schur complement, the straightforward