(Fall EE236A 2007-08)
Lecture 14 Primal-dual interior-point methods
primal-dual path-following Mehrotras corrector step computing the search directions
141
Central path and complementary slackness
s + Ax b = 0 AT z + c = 0 zisi = 1/t, z 0, i = 1, . . . , m s0
continuous deformation of optimality conditions denes central path: solution is x = x(t), s = b Ax(t), zi = 1/tsi m + n linear and m nonlinear equations in the variables s Rm, x Rn, z Rm
Primal-dual interior-point methods
142
Interpretation of barrier method
apply Newtons method to s + Ax b = 0, AT z + c = 0, zi 1/(tsi) = 0, i = 1, . . . , m
i.e., linearize around current x, z , s: 0 AT X A I (Ax + s b) z x = (AT z + c) 0 0 s 0 X 1/t 1/t Xz
where X = diag(s) solution (for s + Ax b = 0, AT z + c = 0): determine x from AT X 2Ax = tc AT X 11 i.e., x is the Newton direction used in barrier method substitute to obtain s, z
Primal-dual interior-point methods 143
Primal-dual path-following methods
modications to the barrier method: dierent linearization of central path update both x and z after each Newton step allow infeasible iterates very aggressive step size selection (99% or 99.9% of step to the boundary) update t after each Newton step (hence distinction between outer & inner iteration disappears) linear or polynomial approximation to the central path limited theory, fewer convergence results work better in practice (faster and more reliable)
Primal-dual interior-point methods
144
Primal-dual linearization
apply Newtons method to s + Ax b = 0 AT z + c = 0 zisi 1/t = 0, i.e., linearize around s, 0A AT 0 X0 i = 1, . . . , m
x, z : I (Ax + s b) z 0 x = (AT z + c) s Z 1/t Xz
where X = diag(s), Z = diag(z ) iterates can be infeasible: b Ax = s, AT z + c = 0 we assume s > 0, z > 0
Primal-dual interior-point methods 145
computing x, z , s 1. compute x from AT X 1ZAx = AT z AT X 11/t rz AT X 1Zrx where rx = Ax + s b, rz = AT z + c 2. s = rx Ax 3. z = X 11/t z X 1Z s the most expensive step is step 1
Primal-dual interior-point methods
146
Ane scaling direction
z a AI (Ax + s b) 0 0 xa = (AT z + c) 0Z Xz sa where X = diag(s), Z = diag(z ) 0 AT X
limit of Newton direction for t Newton step for s + Ax b = 0 AT z + c = 0 zisi = 0, i = 1, . . . , m i.e., the primal-dual optimality conditions
Primal-dual interior-point methods
147
Centering direction
0AI z cent 0 AT 0 0 xcent = 0 scent 1 X0Z where X = diag(s), Z = diag(z ) limit of Newton direction for t 0 search direction is weighted sum of centering direction and ane scaling direction x = (1/t)xcent + xa z = (1/t)z cent + z a s = (1/t)scent + sa in practice: compute ane scaling direction rst choose t compute centering direction and add to ane scaling direction
Primal-dual interior-point methods 148
Heuristic for selecting t
compute ane scaling direction compute primal and dual step lengths to the boundary along the ane scaling direction x = max{ [0, 1] | s + sa 0} z compute = max{ [0, 1] | z + z a 0}
(s + xsa )T (z + z z a ) = sT z small means ane scaling directions are good search directions (signicant reduction in sT z ) use t = m/(sT z ) i.e., search direction will be the Newton direction towards the central point with duality gap sT z a heuristic, based on extensive experiments
Primal-dual interior-point methods 149
3
Mehrotras corrector step
0 AT X 0 AI z cor 0 0 0 xcor = cor a a s 0Z X z
higher-order correction to the ane scaling direction:
a cor (si + sa + scor)(zi + zi + zi ) 0 i i
computation can be combined with centering step, i.e., use x = xcc + xa , where 0 AT X 0 AI z cc 0 0 0 xcc = scc 0Z 1/t X a z a
1410
z = z cc + z a ,
s = scc + sa
Primal-dual interior-point methods
Step size selection
determine step to the boundary x = max{ 0 | s + s 0} z = max{ 0 | z + z 0}
update x, s, z x := x + min{1, 0.99x}x s := s + min{1, 0.99x}s z := z + min{1, 0.99z }z
Primal-dual interior-point methods
1411
Algorithm
choose starting points x, z, s with s > 0, z > 0 1. evaluate stopping criteria primal feasibility: Ax + s b 1(1 + b ) dual feasibility: AT z + c 2(1 + c ) maximum absolute error: cT x + bT z 3 maximum relative error: cT x + bT z 4|bT z | cT x + bT z 4|cT x| if bT z > 0 if cT x < 0
2. compute ane scaling direction (X = diag(s), Z = diag(z )) z a 0AI (Ax + s b) AT 0 0 xa = (AT z + c) X0Z Xz sa
Primal-dual interior-point methods 1412
3. compute steps to the boundary x = max{ [0, 1] | s + sa 0} z = max{ [0, 1] | z + z a 0}
4. compute centering-corrector steps 0 AT X AI z cc 0 0 xcc = sT z scc 0Z 1 X a z a m 0 0
where X a = diag(sa ), and (s + xsa )T (z + z z a ) sT z
3
=
Primal-dual interior-point methods
1413
5. compute search directions x = xa + xcc, s = sa + scc, z = z a + z cc
6. determine step sizes and update x = max{ 0 | s + s 0} z = max{ 0 | z + z 0}
x := x + min{1, 0.99x}x s := s + min{1, 0.99x}s z go to step 1 := z + min{1, 0.99z }z
Primal-dual interior-point methods
1414
Computing the search direction
most expensive part of one iteration: solve two sets of equations AT X 1ZAxa = r1, for some r1, r2 two methods sparse Cholesky factorization: used in all general-purpose solvers conjugate gradients: used for extremely large LPs, or LPs with special structure AT X 1ZAxcc = r2
Primal-dual interior-point methods
1415
Cholesky factorization
if B = B T Rnn is positive denite, then it can be written as B = LLT L lower triangular with lii > 0 L is called the Cholesky factor of B costs O(n3) if B is dense application: solve Bx = d with B = LLT solve Ly = d (forward substitution) solve LT x = y (backward substitution)
Primal-dual interior-point methods
1416
Sparse Cholesky factorization
solve Bx = d with B positive denite and sparse 1. reordering of rows and columns of B to increase sparsity of L 2. symbolic factorization: based on sparsity pattern of B , determine sparsity pattern of L 3. numerical factorization: determine L 4. forward and backward substitution: compute x only steps 3,4 depend on the numerical values of B ; only step 4 depends on the right hand side; most expensive steps: 2,3 in Mehrotras method with sparse LP: B = AT X 1ZA do steps 1,2 once, at the beginning of the algorithm (AT X 1ZA has same sparsity pattern as AT A) do step 3 once per iteration, step 4 twice
Primal-dual interior-point methods 1417
Conjugate gradients
solve Bx = d with B = B T Rnn positive denite iterative method requires n evaluations of Bx (in theory) faster if evaluation of Bx is cheap (e.g., B is sparse, Toeplitz, . . . ) much cheaper in memory than Cholesky factorization less accurate and robust (requires preconditioning) in Mehrotras method: B = AT X 1ZA evaluations Bx are cheap if evaluations Ax and AT y are cheap (e.g., A is sparse)
Primal-dual interior-point methods
1418