On the role of Linear Algebra in the development of Interior Point algorithms and software

1. Due Giorni di Algebra Lineare Numerica Bologna, 6-7 marzo, 2007 On the role of Linear Algebra in the development of Interior Point algorithms and software D. di Serafino,Second University of Naples, Italy daniela.diserafino@unina2.it joint work with S. Cafieri, École Polytechnique, Paris M. D’Apuzzo, V. De Simone, Second University of Naples G. Toraldo, University of Naples “Federico II”

2. Outline • Motivations • KKT system and IP methods • Iterative solution of KKT system • Constraint Preconditioner • Adaptive termination control • Numerical experiments On the role of Linear Algebra in the development of IP algorithms and software

3. Motivations KKT system symmetric indefinte Hsymmetric, Ddiagonal semidefinite positive,Jfull rank • ubiquitous in optimization algorithms • critical issue in an effective implementation of optimization methods Focus on IP methods On the role of Linear Algebra in the development of IP algorithms and software

4. Motivations Development of an IP-based software package for solving large-scale (convex) Quadratic Programming problems x, y, , tprimal and dual variables, s,,zslack variables PRQP - Potential Reduction software for Quadratic Programming problems On the role of Linear Algebra in the development of IP algorithms and software

5. KKT system and IP methods LINEAR ALGEBRA OPTIMIZATION type of problem different blocks formulation of the method in the KKT system large-scale problems sparse direct or iterative solvers local convexity of the problem KKT matrix inertia control accuracy requirements inexact solution of the KKT system adaptive stopping criteria increasing ill conditioning as the iterates approach the solutionpreconditioning techniques On the role of Linear Algebra in the development of IP algorithms and software

6. Infeasible primal-dual PR framework (Kojima, Mizuno & Todd, 1995) potential function(Tanabe,1988; Todd & Ye, 1990) barrier functions duality gap reduce the infeasibility at least at the same rate as the duality gap 0, 0 corresponding to the initial point w0 primal infeasibility dual infeasibility  = 0 feasible version On the role of Linear Algebra in the development of IP algorithms and software

7. PR basic steps • Given the current interior iterate w=(x,y,s,z,,v,t), apply a Newton step to the perturbed KKT conditions • Update w with suitably chosen On the role of Linear Algebra in the development of IP algorithms and software

8. Newton system reduction: KKT system diagonal E accounts for bound constraints equality + bound constraints: D pd inequality + bound constraints: Di pd ineq. + eq. + bound constraints: condensed system bound constraints only (J = I): On the role of Linear Algebra in the development of IP algorithms and software

9. Inherent ill conditioning approaching the solution, some entries ofDand Emay becomevery large, producing anincreasing ill conditioningin the matrix it is crucial to use preconditioning strategies On the role of Linear Algebra in the development of IP algorithms and software

10. KKT system solvers: direct vs iterative • Direct solvers • Widely used in well-established IP software (e.g. Mosek, PCx, Loqo, OOQP, KNITRO-Interior/Direct, IPOPT) • Ill conditioning not a severe problem (S. Wright, 1997; M. Wright, 1998) • Computational cost may become prohibitive for large-scale problems Iterative solvers • Increasing attention by IP community in the last 15 years (implemented in KNITRO-Interior/CG, HOPDM, PRQP) • Require suitable preconditioning techniques • Adaptive termination criteria On the role of Linear Algebra in the development of IP algorithms and software

11. Constraint Preconditioner (CP) • Indefinite preconditioner with the same structure as the KKT matrix • M “simple” approximation to H such that is spd on ker(Je). Common choice: M diagonal. • CP variants investigated by many researchers (Axelsson, 1979; Golub & Wathen, 1998; Luksan & Vlcek, 1998; Keller et al., 2000; Rozloznik & Simoncini, 2002; Durazzi & Ruggiero, 2003; Gondzio et al., 2004; Dollar et al., 2006-2007; di Serafino et al., 2007; Forsgren et al., 2007; …) On the role of Linear Algebra in the development of IP algorithms and software

12. CP: spectral properties • P-1K has at leastmunit eigenvalues • Ifrank(D)=p, then P-1K hasadditional m-p unit eigenvalues • The remaining eigenvalues are real positive (bounds are available) (Keller at al., 2000; Durazzi & Ruggiero, 2003; Bergamaschi et al., 2004; Dollar, 2007) • When the iterate approaches the solution, if qentries of Dget close to zero, then additional qeigenvalues tend to be clustered around 1 • We expect that the preconditioner increases its effectiveness as the IP method progresses On the role of Linear Algebra in the development of IP algorithms and software

13. CP+CG applied to KKT system behaves as CG applied to a reduced system with matrix using as preconditioner (Z basis of ker(Je)) CP: use of Conjugate Gradient (CG) method Starting guess such that Preconditioned Projected CG • No breakdown • Convergence in at most n-m+p iterations, p=rank(D) (Gould, Hribar, Nocedal, 2001; Rozloznik, Simoncini, 2002; Cafieri, D’Apuzzo, De Simone, di Serafino, 2007; Dollar, 2007) On the role of Linear Algebra in the development of IP algorithms and software

14. Building CP Explicit LBLT factorization of P(B block diagonal with 1x1 or 2x2 blocks) LBLT factorization of the Schur complement of -M Implicit Schilders factorization of P (Dollar et al., 2006) May still account for a large part of the computational cost! “Requiring a factorization of P may still be considered a disadvantage, and methods which avoid this are urgently needed ” (Gould, Orban, Toint, Acta Numerica, 2005) On the role of Linear Algebra in the development of IP algorithms and software

15. Reducing the cost of building CP Always set the (2,2)-block to 0 (Dollar, Gould, Schilders, Wathen, 2007) Approximate J by dropping away entries below a prescribed tolerance and outside a fixed band (caseD=0) (Bergamaschi, Gondzio, Venturin, Zilli, 2006) Approximate the Schur complement via incomplete Cholesky factorization (caseD=0) (Benzi, Simoncini, 2006) Reuse CP for some IP iterations (Cafieri, D’Apuzzo, De Simone, di Serafino, 2007b) On the role of Linear Algebra in the development of IP algorithms and software

16. Reusing CP The idea is not new (condensed system: Carpenter & Shanno, 1993; Karmakar & Ramakrishnan, 1991) Reusing CPuse an approximate CP CG cannot be used, apply SQMR computed at a previous outer step On the role of Linear Algebra in the development of IP algorithms and software

17. Reusing CP: spectral properties(work in progress) has • an eigenvalue at 1 with multiplicity at least 2m-p-j • at most 2j eigenvalues with nonzero imaginary part The eigenvalues λsatify On the role of Linear Algebra in the development of IP algorithms and software

18. Reusing CP CP is reused until its effectiveness deteriorates the number of inner iterations must not increase too much the number of outer steps for which CP is reused must be bounded by taking into account the increasing ill conditioning and the accuracy requirements while(PR stopping criterion not satisfied)do … k = k+1 factorize Pk apply to the KKT system SQMR with Pk j = k; l = 0 while(iit(k) ≤α∙ iit(j)and l ≤ β∙lmax)do apply to the KKT system SQMR with P j k = k+1; l = l+1 endwhile lmax = l iit(k) = # inner iter. at outer step k a = 2,b = 1(by numerical experiments) … endwhile On the role of Linear Algebra in the development of IP algorithms and software

19. Termination control • Adaptive stopping criteria, according to the quality of the IP iterate(to reduce the number of inner iterations) • Basic idea: • at early outer iterations low accuracy in solving the KKT system • as the iterates approach the optimal solution the accuracy requirement grows up Regard the IP method as an inexact Newton method: condition for convergence perturbation of the KKT cond. of the original pb. residual of the KKT system KKT cond. of the original pb. On the role of Linear Algebra in the development of IP algorithms and software

20. Termination control Inexact PR convergence Cafieri, D’Apuzzo, De Simone, di Serafino, Toraldo, 2007c polynomial convergence On the role of Linear Algebra in the development of IP algorithms and software

21. Termination control: theory + comput. study CG/SQMR stopping criterion: • Theory • Computational study • Further reduce the number of inner iterations • Avoid slowdown in decreasing the infeasibility (Cafieri, D’Apuzzo, De Simone, di Serafino, 2007d) On the role of Linear Algebra in the development of IP algorithms and software

22. Software architecture of PRQP SIF, AMPL interfaces Pre/post-processing utilities PRQP core linear algebra kernels CG SQMR CP MA27 (LBLT) ICFS BLAS On the role of Linear Algebra in the development of IP algorithms and software

23. Testing details Selected CUTEr problems • PR stopping criteria: δ = relative duality gap, σ= relative infeasibility # PR iterations ≤ 80 • Computational environment • 2.53 GHz Pentium IV, 1.256 GB RAM, 128 KB L1 cache • Linux Ubuntu 4.1.2, g77 3.4.6 and gcc 4.1.3 compilers • sparsity of Hessian and constr. matrices ≥ 99% • * = diagonal Hessian On the role of Linear Algebra in the development of IP algorithms and software

24. PRQP: selectedresults * = diagonal Hessian On the role of Linear Algebra in the development of IP algorithms and software

25. PRQP: selected results * = diagonal Hessian On the role of Linear Algebra in the development of IP algorithms and software

26. Future work • Devising and experimenting new strategies for CP approximation • Applying and analysing CPs in the context of IP methods for nonconvex (quadratic) programming • Developing and experimenting inertia-controlling iterative solvers Thanks for your attention! On the role of Linear Algebra in the development of IP algorithms and software