NASOQ

Quadratic programs (QP), minimizations of quadratic objectives subject to linear inequality and equality constraints, are at the heart of algorithms across scientific domains. Applications include fundamental tasks in geometry processing, simulation, engineering, animation and finance where the accurate, reliable, efficient, and scalable solution of QP problems is critical. However, available QP algorithms generally provide either accuracy or scalability – but not both. Some algorithms reliably solve QP problems to high accuracy but work only for smaller-scale QP problems due to their reliance on dense matrix methods. Alternately, many other QP solvers scale well via sparse, efficient algorithms but cannot reliably deliver solutions at requested accuracies. Towards addressing the need for accurate and efficient QP solvers at scale, we develop NASOQ, a new, full-space QP algorithm that provides accurate, efficient, and scalable solutions for QP problems. To enable NASOQ we construct a new row modification method and fast implementation of LDL factorization for indefinite systems. Together they enable efficient updates and accurate solutions of the iteratively modified KKT systems required for accurate QP solves. While QP methods have been previously tested on large synthetic benchmarks, to test and compare NASOQ’s suitability for real-world applications we collect here a new benchmark set comprising a wide range of graphics-related QPs across physical simulation, animation, and geometry processing tasks. We combine these problems with numerous pre-existing stress-test QP benchmarks to form, to our knowledge, the largest-scale test set of application-based QP problems currently available. Building off of our base NASOQ solver we then develop and test two NASOQ variants against best, state-of-the-art available QP libraries – both commercial and open-source. Our two NASOQ-based methods each solve respectively 98.8% and 99.5% of problems across a range of requested accuracies from 10−3 to 10−9 with average speedups ranging from 1.7× to 24.8× over fastest competing methods.

[1]  H. Huynh A Large-Scale Quadratic Programming Solver Based on Block-Lu Updates of the KKT System , 2008 .

[2]  Stefan Schaal,et al.  Quadratic programming for inverse dynamics with optimal distribution of contact forces , 2012, 2012 12th IEEE-RAS International Conference on Humanoid Robots (Humanoids 2012).

[3]  Doug L. James,et al.  Toward high-quality modal contact sound , 2011, ACM Trans. Graph..

[4]  Robert Bridson,et al.  Blended cured quasi-newton for distortion optimization , 2018, ACM Trans. Graph..

[5]  Stephen P. Boyd,et al.  Differentiable Convex Optimization Layers , 2019, NeurIPS.

[6]  Stephen J. Wright,et al.  Object-oriented software for quadratic programming , 2003, TOMS.

[7]  TIMOTHY A. DAVIS,et al.  Algorithm 1000 , 2019, ACM Transactions on Mathematical Software.

[8]  D K Smith,et al.  Numerical Optimization , 2001, J. Oper. Res. Soc..

[9]  Stephen J. Wright Primal-Dual Interior-Point Methods , 1997, Other Titles in Applied Mathematics.

[10]  Serge Gratton,et al.  A Note on GMRES Preconditioned by a Perturbed L D LT Decomposition with Static Pivoting , 2007, SIAM J. Sci. Comput..

[11]  Kenny Erleben Numerical methods for linear complementarity problems in physics-based animation , 2013, SIGGRAPH '13.

[12]  Joseph W. H. Liu The role of elimination trees in sparse factorization , 1990 .

[13]  J. Gondzio,et al.  Solving Nonlinear Financial Planning ProblemsWith 109 Decision Variables On MassivelyParallel Architectures , 2006 .

[14]  Jennifer A. Scott,et al.  Pivoting strategies for tough sparse indefinite systems , 2013, TOMS.

[15]  Shoaib Kamil,et al.  Sympiler: Transforming Sparse Matrix Codes by Decoupling Symbolic Analysis , 2017, SC17: International Conference for High Performance Computing, Networking, Storage and Analysis.

[16]  Shoaib Kamil,et al.  ParSy: Inspection and Transformation of Sparse Matrix Computations for Parallelism , 2018, SC18: International Conference for High Performance Computing, Networking, Storage and Analysis.

[17]  Dinesh K. Pai,et al.  Staggered projections for frictional contact in multibody systems , 2008, SIGGRAPH Asia '08.

[18]  Timothy A. Davis,et al.  Direct methods for sparse linear systems , 2006, Fundamentals of algorithms.

[19]  Michael A. Saunders,et al.  Inertia-Controlling Methods for General Quadratic Programming , 1991, SIAM Rev..

[20]  Alex Pothen,et al.  Interactively Cutting and Constraining Vertices in Meshes Using Augmented Matrices , 2016, ACM Trans. Graph..

[21]  Robert H. Halstead,et al.  Matrix Computations , 2011, Encyclopedia of Parallel Computing.

[22]  T. Tsuchiya,et al.  On the formulation and theory of the Newton interior-point method for nonlinear programming , 1996 .

[23]  Endong Wang,et al.  Intel Math Kernel Library , 2014 .

[24]  Stephen P. Boyd,et al.  ECOS: An SOCP solver for embedded systems , 2013, 2013 European Control Conference (ECC).

[25]  Elizabeth Wong Active-Set Methods for Quadratic Programming , 2011 .

[26]  Christian Kirches,et al.  qpOASES: a parametric active-set algorithm for quadratic programming , 2014, Mathematical Programming Computation.

[27]  Jorge J. Moré,et al.  Benchmarking optimization software with performance profiles , 2001, Math. Program..

[28]  David I. W. Levin,et al.  Eulerian-on-lagrangian cloth simulation , 2018, ACM Trans. Graph..

[29]  Michael A. Saunders,et al.  User's Guide for SQOPT Version 7.5: Software for Large-Scale Linear and Quadratic Programming , 2016 .

[30]  Stephen P. Boyd,et al.  CVXGEN: a code generator for embedded convex optimization , 2011, Optimization and Engineering.

[31]  Timothy A. Davis,et al.  Row Modifications of a Sparse Cholesky Factorization , 2005, SIAM J. Matrix Anal. Appl..

[32]  James F. O'Brien,et al.  Updated sparse cholesky factors for corotational elastodynamics , 2012, TOGS.

[33]  M. Mahdavi,et al.  ARTICLE IN PRESS Available online at www.sciencedirect.com , 2007 .

[34]  Donald Goldfarb,et al.  A numerically stable dual method for solving strictly convex quadratic programs , 1983, Math. Program..

[35]  Stephen P. Boyd,et al.  OSQP: an operator splitting solver for quadratic programs , 2017, 2018 UKACC 12th International Conference on Control (CONTROL).

[36]  William W. Hager,et al.  Updating the Inverse of a Matrix , 1989, SIAM Rev..

[37]  Yanran Ding,et al.  qpSWIFT: A Real-Time Sparse Quadratic Program Solver for Robotic Applications , 2019, IEEE Robotics and Automation Letters.

[38]  GingoldYotam,et al.  Interactive Design and Stability Analysis of Decorative Joinery for Furniture , 2017 .

[39]  C. Mészáros,et al.  A repository of convex quadratic programming problems , 1999 .

[40]  Lorenz T. Biegler,et al.  On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming , 2006, Math. Program..

[41]  Olga Sorkine-Hornung,et al.  Ink-and-ray: Bas-relief meshes for adding global illumination effects to hand-drawn characters , 2014, TOGS.

[42]  Daniele Panozzo,et al.  libigl: prototyping geometry processing research in C++ , 2017, SIGGRAPH ASIA.

[43]  Eitan Grinspun,et al.  Reflections on simultaneous impact , 2012, ACM Trans. Graph..

[44]  Gustaf Söderlind,et al.  Digital filters in adaptive time-stepping , 2003, TOMS.

[45]  Nicholas I. M. Gould,et al.  GALAHAD, a library of thread-safe Fortran 90 packages for large-scale nonlinear optimization , 2003, TOMS.

[46]  M. Powell On the quadratic programming algorithm of Goldfarb and Idnani , 1985 .

[47]  Stephen P. Boyd,et al.  Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers , 2011, Found. Trends Mach. Learn..

[48]  Gene H. Golub,et al.  Numerical solution of saddle point problems , 2005, Acta Numerica.

[49]  Alec Jacobson,et al.  Seamless reconstruction of part-based high-relief models from hand-drawn images , 2018, Expressive.

[50]  Nancy Wilkins-Diehr,et al.  XSEDE: Accelerating Scientific Discovery , 2014, Computing in Science & Engineering.

[51]  Yousef Saad,et al.  Iterative methods for sparse linear systems , 2003 .

[52]  YANQING CHEN,et al.  Algorithm 8 xx : CHOLMOD , supernodal sparse Cholesky factorization and update / downdate ∗ , 2006 .

[53]  Marc Alexa,et al.  Localized solutions of sparse linear systems for geometry processing , 2017, ACM Trans. Graph..

[54]  Olaf Schenk,et al.  Two-level dynamic scheduling in PARDISO: Improved scalability on shared memory multiprocessing systems , 2002, Parallel Comput..

[55]  Timothy A. Davis,et al.  Dynamic Supernodes in Sparse Cholesky Update/Downdate and Triangular Solves , 2009, TOMS.

[56]  Pascal Hénon,et al.  PaStiX: a high-performance parallel direct solver for sparse symmetric positive definite systems , 2002, Parallel Comput..

[57]  Nicholas I. M. Gould,et al.  On the Solution of Equality Constrained Quadratic Programming Problems Arising in Optimization , 2001, SIAM J. Sci. Comput..

[58]  J. Zico Kolter,et al.  OptNet: Differentiable Optimization as a Layer in Neural Networks , 2017, ICML.

[59]  Mihai Anitescu,et al.  Using Krylov subspace and spectral methods for solving complementarity problems in many‐body contact dynamics simulation , 2013 .

[60]  Doug L. James,et al.  Real-time reduced large-deformation models and distributed contact for computer graphics and haptics , 2007 .

[61]  Olga Sorkine-Hornung,et al.  Bounded biharmonic weights for real-time deformation , 2011, Commun. ACM.

[62]  Timothy A. Davis,et al.  Modifying a Sparse Cholesky Factorization , 1999, SIAM J. Matrix Anal. Appl..

[63]  R. Fletcher Practical Methods of Optimization , 1988 .

[64]  Iain S. Duff,et al.  MA57---a code for the solution of sparse symmetric definite and indefinite systems , 2004, TOMS.