High-Order Methods for Overcoming Degeneracy in Interior-Point Methods for Quadratic Programming Nicholas Gould 1

Dominique Orban 2 1 2

3

Daniel Robinson 3

Rutherford Appleton Laboratory

École Polytechnique de Montréal

Oxford University Mathematical Institute

Santiago, Chile July 29, 2010

Outline

Introduction and Motivation

Problem of interest QP

minimize cTx + 12 xTHx x

subject to Ax = b, x ≥ 0

Interested in general H, but require H  0 for some results Assume that A has full row rank

First order optimality conditions Hx∗ + c = AT y∗ + z∗ Ax∗ = b min(x∗ , z∗ ) = 0 max(x∗i , z∗i ) > 0 for all i =⇒ nondegenerate max(x∗i , z∗i ) = 0 for some i =⇒ degenerate

Problem of interest QP

minimize cTx + 12 xTHx x

subject to Ax = b, x ≥ 0

Interested in general H, but require H  0 for some results Assume that A has full row rank

First order optimality conditions Hx∗ + c = AT y∗ + z∗ Ax∗ = b min(x∗ , z∗ ) = 0 max(x∗i , z∗i ) > 0 for all i =⇒ nondegenerate max(x∗i , z∗i ) = 0 for some i =⇒ degenerate

Problem of interest QP

minimize cTx + 12 xTHx x

subject to Ax = b, x ≥ 0

Interested in general H, but require H  0 for some results Assume that A has full row rank

First order optimality conditions Hx∗ + c = AT y∗ + z∗ Ax∗ = b min(x∗ , z∗ ) = 0 max(x∗i , z∗i ) > 0 for all i =⇒ nondegenerate max(x∗i , z∗i ) = 0 for some i =⇒ degenerate

Interior point subproblem T

minimize c x + x

1 T x Hx−µ 2

n X

ln(xi )

subject to Ax = b

i=1

µ > 0 is the barrier parameter
if x(µ), y(µ) is the minimizer, then generally
lim x(µ), y(µ), z(µ) → (x∗ , y∗ , z∗ ) µ→0

def

def

where z(µ) = µX(µ)−1 e and X = diag(x1 , x2 , . . . , xn )
Notation: v = (x, y, z), v(µ) = x(µ), y(µ), z(µ) , . . . etc. Key result (Stoer/Wechs/Mizuno 1998, Wright/Orban 2002) ( Θ(µ) nondegenerate ∗ v(µ) − v = √ Θ( µ) degenerate

Interior point subproblem T

minimize c x + x

1 T x Hx−µ 2

n X

ln(xi )

subject to Ax = b

i=1

µ > 0 is the barrier parameter
if x(µ), y(µ) is the minimizer, then generally
lim x(µ), y(µ), z(µ) → (x∗ , y∗ , z∗ ) µ→0

def

def

where z(µ) = µX(µ)−1 e and X = diag(x1 , x2 , . . . , xn )
Notation: v = (x, y, z), v(µ) = x(µ), y(µ), z(µ) , . . . etc. Key result (Stoer/Wechs/Mizuno 1998, Wright/Orban 2002) ( Θ(µ) nondegenerate ∗ v(µ) − v = √ Θ( µ) degenerate

Our motivation ( Θ(µ) v(µ) − v∗ = √ Θ( µ) 1

nondegenerate degenerate

“one off” solve of problem QP - µ = 1.0−6 =⇒ kv(µ) − v∗ k ≈ 1.0−3 in the degenerate case

2

improve cross-over methods for solving QP - use x(µ) to predict optimal active set for an active set phase

3

improve large-scale SQP methods (S2QP) - solve a convex QP to predict an optimal active set, which is used to compute an additional accelerator step

4

efficiency - higher-order methods =⇒ fewer factorizations

Nondegenerate QP

Trajectory of minimizers (Fiacco and McCormick) Let v∗ be a nondegenerate second-order sufficient local minimizer of QP and assume that A has full row rank, then there exists a continuously differentiable path v(µ) such that lim v(µ) = v∗

µ→0

and v(µ) minimizes the barrier problem for all µ sufficiently small 0.6

0.5

0.4

0.3

0.2

0.1

0

−0.1 −0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

Question: Given v(µ) and µ ¯ ∈ (0, µ), how do we find v(µ) ¯ efficiently? Answer: Apply Newton’s Method to perturbed optimality conditions

Barrier optimality conditions 1

Hx + c =

2

Ax = b µ ¯ Z−1 e

3

x=

4

min(x, z) > 0

H −AT  A 0 Z(µ) 0 

ATy

+z

Equivalent conditions 1

Hx + c = ATy + z

2

Ax = b

3

Xz = µe ¯

4

min(x, z) > 0

  N   −I ∆x Hx(µ) + c − ATy(µ) − z(µ)  0  ∆yN  = −  Ax(µ) − b N X(µ) ∆z X(µ)z(µ) − µe ¯

Question: Given v(µ) and µ ¯ ∈ (0, µ), how do we find v(µ) ¯ efficiently? Answer: Apply Newton’s Method to perturbed optimality conditions

Barrier optimality conditions 1

Hx + c =

2

Ax = b µ ¯ Z−1 e

3

x=

4

min(x, z) > 0

H −AT  A 0 Z(µ) 0 

ATy

+z

Equivalent conditions 1

Hx + c = ATy + z

2

Ax = b

3

Xz = µe ¯

4

min(x, z) > 0

  N   −I ∆x Hx(µ) + c − ATy(µ) − z(µ)  0  ∆yN  = −  Ax(µ) − b N X(µ) ∆z X(µ)z(µ) − µe ¯

Question: Given v(µ) and µ ¯ ∈ (0, µ), how do we find v(µ) ¯ efficiently? Answer: Apply Newton’s Method to perturbed optimality conditions

Barrier optimality conditions 1

Hx + c =

2

Ax = b µ ¯ Z−1 e

3

x=

4

min(x, z) > 0

H −AT  A 0 Z(µ) 0 

ATy

+z

Equivalent conditions 1

Hx + c = ATy + z

2

Ax = b

3

Xz = µe ¯

4

min(x, z) > 0

  N   −I ∆x Hx(µ) + c − ATy(µ) − z(µ)  0  ∆yN  = −  Ax(µ) − b N X(µ) ∆z X(µ)z(µ) − µe ¯

Question: Given v(µ), is there a method to compute accurate estimates of v(µ) ¯ for 0 < µ ¯  µ? Answer: Yes, use a higher-order Taylor approximation for v(µ), ¯ i.e., v(µ) ¯ ≈ v(µ) +

p X v[j] (µ)(µ ¯ − µ)j j=1

j!

def

= vp (µ) ¯

v1 (µ) ¯ ≡ ∆vN existence of derivatives follows from Fiacco and McCormick derivatives are obtained by repeated differentiation of Hx(µ) + c = ATy(µ) + z(µ) Ax(µ) = b X(µ)z(µ) = µe these results assume that v(µ) is computed exactly

Question: Given an approximation to v(µ), is there a method to compute accurate approximations to v(µ) ¯ for 0 < µ ¯  µ? Answer: Yes, use approximations to multiple points on the trajectory compute (µj , xj ) such that xj ≈ x(µj ) for j = 1, . . . p solve the least-squares problem minimize c

1 2

p X

pr (µj ) − xj


2

j=1 def

for the coefficient vector c, where pr (µ) =

Pr

l=0 cl µ

l

and r ≤ p

use the approximation pr (µ) ¯ ≈ x(µ) ¯ For most of the talk we focus on higher-order approximations, but key concepts may be adapted to a least-squares approach

Question: Given an approximation to v(µ), is there a method to compute accurate approximations to v(µ) ¯ for 0 < µ ¯  µ? Answer: Yes, use approximations to multiple points on the trajectory compute (µj , xj ) such that xj ≈ x(µj ) for j = 1, . . . p solve the least-squares problem minimize c

1 2

p X

pr (µj ) − xj


2

j=1 def

for the coefficient vector c, where pr (µ) =

Pr

l=0 cl µ

l

and r ≤ p

use the approximation pr (µ) ¯ ≈ x(µ) ¯ For most of the talk we focus on higher-order approximations, but key concepts may be adapted to a least-squares approach

Nondegenerate example minimize x

1 2 x 2

subject to x ≥ 1

Errors for various approximation schemes µ = 1.0e−5 µ 10 −10

x(0)

x µ2

Interp 2

9.99e−10

9.99e−10

Interp 3

−12

1.99e

−12

1.99e

1.77e−12

Interp 4

4.86e−14

4.84e−14

4.30e−14

Hermite 2

0.00e+00

0.00e+00

0.00e+00

Spline 2

9.99e−10

9.99e−10

8.90e−10

Spline 3

1.99e−12

1.99e−12

1.77e−12

Spline 4

−14

4.86e

−14

4.86e

4.30e−14

Taylor 2

1.99e−15

1.99e−15

1.33e−15




x




8.90e

Degenerate QP

Degenerate example minimize x

1 2 x 2

subject to 0 ≤ x ≤ 1

Errors for various approximation schemes µ = 1.0e−5

Interp 4

x(0)

x µ2

2.34e−03

2.33e−03

−03

−03




x

µ 10 −03




1.42e

Hermite 3

1.40e

1.39e

5.99e−04

Hermite 4

1.40e−03

1.39e−03

5.98e−04

Spline 4

2.34e−03

2.33e−03

1.42e−03

Taylor 2

1.18e−03

1.17e−03

4.19e−04

Taylor 3

9.88e−04

9.78e−04

2.75e−04

Taylor 4

−04

8.64e

−04

8.54e

1.93e−04

Taylor 5

7.78e−04

7.68e−04

1.42e−04

The main idea – go Puiseux!

√ Based on the result that v(µ) − v∗ = Θ( µ), we postulate v(µ) = v∗ +

∞ X

cj µ j/2 ,

j=1

and by defining ρ =

√

µ, we have def

v¯(ρ) = v∗ +

∞ X

cj ρ j = v(µ)

j=1

Hope that v¯(ρ) is analytic Use higher-order or least-squares approximations to v¯(ρ)

The main idea – go Puiseux!

√ Based on the result that v(µ) − v∗ = Θ( µ), we postulate v(µ) = v∗ +

∞ X

cj µ j/2 ,

j=1

and by defining ρ =

√

µ, we have def

v¯(ρ) = v∗ +

∞ X

cj ρ j = v(µ)

j=1

Hope that v¯(ρ) is analytic Use higher-order or least-squares approximations to v¯(ρ)

Degenerate example revisited minimize x

1 2 x 2

subject to 0 ≤ x ≤ 1

Errors for various approximation schemes µ = 1.0e−5

Interp 4

x(0)

x µ2

6.99e−09

6.96e−09

−13

−13




x

µ 10 −09




4.15e

Hermite 3

4.87e

4.82e

1.73e−13

Hermite 4

1.27e−14

1.26e−14

4.44e−15

Spline 4

6.99e−09

6.96e−09

4.15e−09

Taylor 2

1.19e−08

1.17e−08

3.79e−09

Taylor 3

5.49e−13

5.42e−13

1.29e−13

Taylor 4

−13

1.38e

−13

1.36e

2.07e−14

Taylor 5

4.97e−16

4.88e−16

5.05e−17

Comparison of “µ” and “ρ” parametrization. µ = 1.0e−5

Interp 4

x(0)

x(µ)
x µ2

x

2.34e−03

2.33e−03

1.42e

−03

−03

µ 10 −03




−04

x(0)

x¯(ρ)
x µ2

x

6.99e−09

6.96e−09

4.15e

−13

−13

µ 10 −09




Hermite 3

1.40e

1.39e

5.99e

4.87e

4.82e

1.73e−13

Hermite 4

1.40e−03

1.39e−03

5.98e−04

1.27e−14

1.26e−14

4.44e−15

Spline 4

2.34e−03

2.33e−03

1.42e−03

6.99e−09

6.96e−09

4.15e−09

Taylor 2

−03

1.18e

−03

1.17e

−04

4.19e

−08

1.19e

−08

1.17e

3.79e−09

Taylor 3

9.88e−04

9.78e−04

2.75e−04

5.49e−13

5.42e−13

1.29e−13

Taylor 4

8.64e−04

8.54e−04

1.93e−04

1.38e−13

1.36e−13

2.07e−14

Taylor 5

7.78e−04

7.68e−04

1.42e−04

4.97e−16

4.88e−16

5.05e−17

Give credit where credit is due We are not the first, but few seem aware of this result Stoer/Wechs/Mizuno (1998), proved that x¯(ρ) is analytic for the sufficient horizontal complementarity problem Stoer/Wechs/Mizuno (1998) suggest different strategies for using x¯(ρ) and indicator functions as the basis of an algorithm Zhao/Sun (1999), suggested a “predictor-corrector” method and a single step method based on perturbed trajectories We do not know of any numerical results We do not know of any software that incorporates this idea We have our own ideas for implementation

Implementation

Strategy 1: the “lunge” 1

Compute v(µ) for µ ≈ 10−6 using an interior point solver

2

Approximate v∗ by a sixth degree Puiseux expansion

3

Obtain an estimate of an optimal active set

4

Verify optimality

5

Use “cross-over” to resolve multipliers for general constraints May use Puiseux expansion to “skip down the path” Strategy has been added to QPB Highly accurate solutions are typically returned Active set identification appears to be improved Requires staying “close” to the path QPB is a two stage method for solving indefinite QP

Strategy 1: the “lunge” 1

Compute v(µ) for µ ≈ 10−6 using an interior point solver

2

Approximate v∗ by a sixth degree Puiseux expansion

3

Obtain an estimate of an optimal active set

4

Verify optimality

5

Use “cross-over” to resolve multipliers for general constraints May use Puiseux expansion to “skip down the path” Strategy has been added to QPB Highly accurate solutions are typically returned Active set identification appears to be improved Requires staying “close” to the path QPB is a two stage method for solving indefinite QP

Strategy 2: perturbed trajectory + Zhang + Zhao/Sun 1

Given v ≈ v(µ) where µ = (xT z)/n, compute σ ∈ (0, 1)

2

Approximate v(σµ) by a Puiseux expansion

3

Perform linesearch to ensure + +

X z

+

≥ γµ

+

and µ

 + 

Hx + c − AT y+ − z+

≥ν

Ax+ − b

v

Single phase v(µ)

Intended for convex QP Higher order friendly

v

Dynamic updating of µ

+

v(σµ)

CQP in GALAHAD

May use final “lunge”

v∗

Strategy 2: perturbed trajectory + Zhang + Zhao/Sun 1

Given v ≈ v(µ) where µ = (xT z)/n, compute σ ∈ (0, 1)

2

Approximate v(σµ) by a Puiseux expansion

3

Perform linesearch to ensure + +

X z

+

≥ γµ

+

and µ

 + 

Hx + c − AT y+ − z+

≥ν

Ax+ − b

v

Single phase v(µ)

Intended for convex QP Higher order friendly

v

Dynamic updating of µ

+

v(σµ)

CQP in GALAHAD

May use final “lunge”

v∗

Algorithm CQP on degenerate problem using Taylor expansion It 0 1 2 3 4 5 . 20 21 22 23 24 25 26 27 28 29

p-feas 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 . 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00

d-feas 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 . 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00 3.2E-27 8.1E-28 4.0E-28 2.0E-28

com-slk 1.0E+00 1.5E-01 2.1E-02 3.1E-03 4.6E-04 6.7E-05 . 1.2E-17 1.7E-18 2.3E-19 3.3E-20 4.6E-21 6.5E-22 9.1E-23 1.3E-23 1.8E-24 2.5E-25

obj 5.0E-01 7.3E-02 1.1E-02 1.6E-03 2.3E-04 3.4E-05 . 5.9E-18 8.3E-19 1.2E-19 1.6E-20 2.3E-21 3.2E-22 4.5E-23 6.4E-24 9.0E-25 1.3E-25

Optimal solution value = 5.0259E-13

step 1.0E+00 1.0E+00 1.0E+00 1.0E+00 1.0E+00 . 1.0E+00 1.0E+00 1.0E+00 1.0E+00 1.0E+00 1.0E+00 1.0E+00 1.0E+00 1.0E+00 1.0E+00

mu 1.0E-02 1.5E-03 2.1E-04 3.1E-05 4.6E-06 5.5E-07 . 4.0E-26 2.1E-27 1.1E-28 5.9E-30 3.1E-31 1.6E-32 8.7E-34 4.6E-35 2.4E-36 1.3E-37

Algorithm CQP on degenerate problem using Puixeux expansion It 0 1 2 3 4

p-feas 0.0E+00 0.0E+00 0.0E+00 0.0E+00 0.0E+00

d-feas 0.0E+00 0.0E+00 1.7E-18 8.1E-20 0.0E+00

com-slk 1.0E+00 2.2E-04 5.0E-08 5.6E-15 7.1E-29

obj 5.0E-01 1.1E-04 2.5E-08 2.8E-15 3.5E-29

Optimal solution value = 8.4205E-15

step 1.0E+00 1.0E+00 1.0E+00 1.0E+00

mu 1.0E-02 2.2E-06 1.1E-11 4.2E-22 6.0E-43

Summary

Summary Highly accurate solution of degenerate QP is not possible by standard higher order or curve fitting methods Degeneracy in QP may be overcome by using a “fractional” Puiseux expansion that is analytic (convex QP) Posed two strategies for incorporating the Puiseux expansion into interior-point framework Appears to work well for nonconvex QP, although not covered by existing theory Extends to classes of LCP Higher order Puiseux expansions improve efficiency by reducing the number of overall factorizations, although generally more system solves are required

References

Josef Stoer, Martin Wechs, and Shinji Mizuno High order infeasible-interior-point methods for solving sufficient linear complementarity problems Mathematics of Operations Research, 23, pp. 832-862, (1998) Yin Zhang On the convergence of a class of infeasible interior-point methods for the horizontal linear complementarity problem Siam J. Optimization, 1 (1994), pp. 208–227 Gongyun Zhao and Jie Sun On the rate of local convergence of high-order infeasible path-following algorithms for P∗ -linear complementarity problems Computational Optimization and Applications, 14, 293-307 (1999)