Optimization Methods - Barrier Interior Point Algorithms

1 Outline 1. 2. 3. 4. 5. 6. Slide 1 Barrier Methods The Central Path Approximating the Central Path The Primal Barrier Algorithm The Primal-Dual Barrier Algorithm Computational Aspects of IPMs 2 Barrier methods min f (x) s:t: gj (x) 0
hi (x) = 0
Slide 2 j = 1
: : :
p i = 1
: : :
m S = fxj gj (x) < 0
j = 1
: : :
p
hi(x) = 0
i = 1
: : :
mg 2.1 Strategy
A barrier function G(x) is a continous function with the property that is approaches 1 as one of gj (x) approaches 0 from negative values.
Examples: G(x) = ; p X j �1 log(;gj (x))
G(x) = ; Slide 3 p X 1 g (x) j �1 j Slide 4
Consider a sequence of k : 0 < k+1 < k and k ! 0.
Consider the problem xk = argminx2S f (x) + k G(x)

Theorem If Every limit point xk generated by a barrier method is a global minimum of the original constrained problem. 1 . . . . . . x* c x(0.01) central path x(0.1) x(1) x(10) analytic center 2.2 Primal path-following IPMs for LO (P ) min c x s:t: Ax = b (D) max b p s:t: A p + s = c 0 0 0 x0 Barrier problem: Slide 5 s0 min B (x) = c x ;  0 s:t: Ax = b n X j �1 log xj Minimizer: x() 3 Central Path
As  varies, minimizers x() form the central path
lim !0 x() exists and is an optimal solution x to the initial LP
For  = 1, x(1) is called the analytic center min ; n X j �1 Slide 6 log xj s:t: Ax = b Slide 7 2 x3 Q (1/2, 0, 1/2) the analytic center of Q . (1/3, 1/3, 1/3) the analytic center of P P . x2 the central path x1
Q = x j x = (x1
0
x3)
x1 + x3 = 1
x  0g
set of optimal solutions to original LP
The analytic center of Q is (1=2
0
1=2) 3.1 Example min x2 s:t: x1 + x2 + x3 = 1 x1
x2
x3  0 min x2 ;  log x1 ;  log x2 ;  log x3 s:t: x1 + x2 + x3 = 1 Slide 8 min x2 ;  log x1 ;  log x2 ;  log(1 ; x1 ; x2): x1() = 1 ; x22() p 1 + 92 + 2 1 + 3  ; x2() = 2 1 ; x () 2 x3() = 2 The analytic center: (1=3
1=3
1=3) Slide 9 3.2 Solution of Central Path Slide 10 � Barrier problem for dual: n X pb+ log s =1 s:t: p A + s = c max 0 0 3 j 0 0 j � Solution (KKT): Ax( ) = b x( ) � 0 A p( ) + s( ) = c s( ) � 0 X ( )S ( )e = e 0
Theorem: If x p and s satisfy optimality conditions, then they are ,  ,  Slide 11 optimal solutions to problems primal and dual barrier problems.
Goal: Solve barrier problem min B (x) = c x ;  0 s:t: Ax = b n X j �1 log xj 4 Approximating the central path @B (x) = c ;  i @xi xi 2 @ B (x) =  @x2i x2i @ 2 B (x) = 0
i = 6 j @xi @xj Given a vector x > 0: @B (x) d + @xi i i�1 n @ 2 B (x) d d 1X 2 i
j �1 @xi @xj i j B (x + d)  B (x) + X = diag(x1
: : :
xn) n X Slide 12 Slide 13 = B (x) + (c0 ; e0 X ;1 )d + 21 d0 X ;2 d Slide 14 Approximating problem: min (c0 ; e0X ;1)d + 12 d0X ;2d s:t: Ad = 0 Solution (from Lagrange): c ; X ;1e + X ;2d ; A0p = 0 Ad = 0 4 Slide 15 . . . . . . x* c x(0.01) central path x(0.1) x(1) x(10) analytic center
System of m + n linear equations, with m + n unknowns (dj
j = 1
: : :
n, and pi
i = 1
: : :
m).
Solution:   d() = I ; X 2A0(AX 2A0);1A xe ; 1 X 2c p() = (AX 2A0);1A(X 2c ; xe)  4.1 The Newton connection
d() is the Newton direction process of calculating this direction is called Slide 16 a Newton step
Starting with x, the new primal solution is x + d() ; 
The corresponding dual solution becomes (p
s) = p()
c ; A0p()
We then decrease  to  = , 0 <  < 1 4.2 Geometric Interpretation
Take one Newton step so that x would be close to x()
Measure of closeness      1 X Se ; e 
   0 <  < 1, X = diag(x1
: : :
xn) S = diag(s1
: : :
sn )
As  ! 0, the complementarity slackness condition will be satised 5 Slide 17 Slide 18 5 The Primal Barrier Algorithm Input (a) (A
b
c) A has full row rank (b) x0 > 0
s0 > 0
p0  (c) optimality tolerance > 0 (d) 0 , and , where 0 <  < 1. 1. (Initialization) Start with some primal and dual feasible x0 > 0
s0 > 0
p0, and set k = 0: 2. (Optimality test) If (sk )0xk < stop else go to Step 3. 3. Let X k = diag(xk1
: : :
xnk )
k+1 = k 4. (Computation of directions) Solve the linear system k+1X ;k 2 d ; A0p = k+1 X ;k 1 e ; c Ad = 0 5. (Update of solutions) Let Slide 19 Slide 20 Slide 21 xk+1 = xk + d
pk+1 = p
sk+1 = c ; A0p: 6. Let k := k + 1 and go to Step 2. 5.1 Correctness p Theorem Given  = 1 ; p  ;p < 1, (x0
s0
p0 ), (x0 > 0, s0 > 0):  + n
      1 X 0S 0 e ; e :   0 Then, after p p 00 0 K = p+; n log (s )(1x ;(1)+ ) iterations, (xK
sK
pK ) is found: (sK )0 xK : 6 Slide 22 5.2 Complexity
Work per iteration involves solving a linear system with m + n equations Slide 23 in m + n unknowns. Given that m n, the work per iteration is O(n3).
0 = (s0 )0 x0 : initial duality gap. Algorithm needs O p  n log 0 iterations to reduce the duality gap from 0 to , with O(n3) arithmetic operations per iteration. 6 The Primal-Dual Barrier Algorithm 6.1 Optimality Conditions Ax() x() A p() + s() s() 0 = b Slide 24  0 = c  0 sj ()xj () =  or X ()S()e = e ;  X () = diag x1()
: : :
xn() , S() = diag;s1()
: : :
sn() 6.2 Solving Equations z = (x
p
s), r = 2n + m Ax ; b 3 F (z) = 64 A0p + s ; c 57 X Se ; e Solve 2 Slide 25 F (z) = 0 6.2.1 Newton's method F (zk + d)  F (zk ) + J (zk)d Slide 26 Here J (z k ) is the r  r Jacobian matrix whose (i
j)th element is given by  @Fi (z )  @zj z �z k F (zk) + J (zk )d = 0 Set z k+1 = z k + d (d is the Newton direction) 7 Slide 27 (xk
pk
sk ) current primal and dual feasible solution Newton direction d = (dkx
dkp
dks ) 2 A 0 0 6 0 A0 I 4 Sk 0 X k 32 k 3 2 dx 76 7 6 5 4 dkp 5 = ; 4 dks Axk ; b 3 A0pk + sk ; c 57 X kSk e ; ke 6.2.2 Step lengths Slide 28 xk+1 = xk + Pk dxk pk+1 = pk + Dk dpk sk+1 = sk + Dk dks To preserve nonnegativity, take    xki = min 1
 fij(dmin ; (dkx )i
) <0g    ski ; Dk = min 1
 fij(min (dks )i
d ) <0g Pk k x i k s i 0<<1 6.3 The Algorithm Slide 29 1. (Initialization) Start with x0 > 0
s0 > 0
p0, and set k = 0 2. (Optimality test) If (sk )0xk < stop else go to Step 3. 3. (Computation of Newton directions) k 0 k k = (s n) x X k = diag(xk1
: : :
xnk ) Sk = diag(sk1
: : :
skn) Solve linear system 2 A 0 0 6 0 A0 I 4 Sk 0 X k 32 k 3 2 dx 76 7 6 5 4 dkp 5 = ; 4 dks 8 Axk ; b A0pk + sk ; c X kSk e ; k e 3 7 5 Slide 30 4. (Find step lengths)    xik Pk = min 1
 fij(dmin ; (dkx )i ) <0g   k  s i k D = min 1
 fij(dmin ; (dk ) ) <0g s i k x i k s i 5. (Solution update) xk+1 = xk + Pk dkx pk+1 = pk + Dk dkp sk+1 = sk + Dk dsk 6. Let k := k + 1 and go to Step 2 6.4 Insight on behavior � Slide 31 Ane Scaling   dane = ;X 2 I ; A (AX 2 A ) 1 AX 2 c � Primal barrier dprimal ; barrier = � � For = 1 Note that  0 �   I ; X 2 A (AX 2 A ) 1 A Xe ; 1 X 2 c 0 ; 0   dcentering = I ; X 2 A (AX 2 A ) 1 A X e dprimal 0 barrier = dcentering + ; � 0 ; 0 ; 1d ane When is large, then the centering direction dominates, i.e., in the beginning, the barrier algorithm takes steps towards the analytic center When is small, then the ane scaling direction dominates, i.e., towards the end, the barrier algorithm behaves like the ane scaling algorithm 7 Computational aspects of IPMs Simplex vs. Interior point methods (IPMs)
Simplex method tends to perform poorly on large, massively degenerate problems, whereas IP methods are much less aected.
Key step in IPMs ; AX 2k A0d = f 9 Slide 32
In implementations of IPMs AX 2k A0 is usually written as AX 2kA0 = LL0
where L is a square lower triangular matrix called the Cholesky factor
Solve system ; AX 2k A0d = f by solving the triangular systems Ly = f
L0d = y
The construction of L requires O(n3) operations but the actual compu- tational eort is highly dependent on the sparsity (number of nonzero entries) of L
Large scale implementations employ heuristics (reorder rows and columns of A) to improve sparsity of L. If L is sparse, IPMs are stronger. 8 Conclusions
IPMs represent the present and future of Optimization.
Very successful in solving very large problems.
Extend to general convex problems 10 Slide 33