Releasing a modernized `QuadProg` for strictly convex quadratic programming

Announcements
1

Hej everyone,

I’m happy to announce that the modernized version of QuadProg is now officially released! This is an updated version of the quadprog solver initially written by Berwin A. Turlach in Fortran 77. It can be used to solve strictly convex quadratic programs of the form

\begin{aligned} \mathrm{minimize} & \quad \dfrac12 \mathbf{x}^T \mathbf{Px} - \mathbf{x}^T \mathbf{q} \\ \mathrm{subject~to} & \quad \mathbf{Ax} = \mathbf{b} \\ & \quad \mathbf{Cx} \geq \mathbf{d}. \end{aligned}

using an active set method. It is most efficient for small to moderate sized QP described using dense matrices. A specialized implementation for sparse constraint matrices is also available.

Updates to the original code include:

  • Sources have been translated from FORTRAN 77 fixed-form to Fortran 90 free-form.
  • All obsolescent features (goto, etc) have been removed. It is now 100% standard compliant (Fortran 2018).
  • It makes use of derived-type and easy to use interfaces. The qp_problem class is used to defined the quadratic program and solve to compute its solution.
  • Calls to blas functions and subroutines now replace some hand-crafted implementations for improved performances.
  • Calls to lapack subroutines now replace the original functionalities provided by linpack.
  • Utility solvers for non-negative least-squares (nnls) and bounded-variables least-squares (bvls) are provided.

Example

The following program solves the QP

\begin{aligned} \mathrm{minimize} & \quad \dfrac12 \left( x_1^2 + x_2^2 + x_3^2 \right) - 5 x_2 \\ \mathrm{subject~to} & \quad -4x_1 + 2x_2 \geq -8 \\ & \quad -3x_1 + x_2 - 2x_3 \geq -2 \\ & \quad x_3 \geq 0. \end{aligned} 
program example
    use quadprog
    implicit none
    integer, parameter :: dp = selected_real_kind(15, 307)
    ! Size of the problem.
    integer, parameter :: n = 3
    ! Quadratic cost.
    real(dp) :: P(n, n), q(n)
    ! Inequality constraints.
    real(dp) :: C(n, n), d(n)
    ! Convenience types.
    type(qp_problem) :: prob
    type(OptimizeResult) :: solution
    ! Miscellaneous
    integer :: i

    !> Setup the quadratic function..
    P = 0.0_dp ; q = [0.0_dp, 5.0_dp, 0.0_dp]
    do i = 1, n
        P(i, i) = 1.0_dp
    enddo
    
    !> Setup the inequality constraints.
    C(:, 1) = [-4.0_dp, 2.0_dp, 0.0_dp]
    C(:, 2) = [-3.0_dp, 1.0_dp, -2.0_dp]
    C(:, 3) = [0.0_dp, 0.0_dp, 1.0_dp]
    d = [-8.0_dp, 2.0_dp, 0.0_dp]

    !> Solve the inequality constrained QP.
    prob = qp_problem(P, q, C=C, d=d)
    solution = solve(prob)

    if (solution%success) then
        print *, "x   =", solution%x ! Solution of the QP.
        print *, "y   =", solution%y ! Lagrange multipliers.
        print *, "obj =", solution%obj ! Objective function.
    endif
end program

Two additional examples for simple Model Predictive Control (double integrator with actuator saturation) are included in the example folder.

Other features

  • Except for blas and lapack, the code base is entirely dependence-free and thus very easy to use.
  • For the moment, only an fpm-compatible build is provided. Whenever I’ll have time, I’ll setup cmake.
  • Non-negative least-squares can be solved as easily as x = nnls(A, b).
  • Bounded-variables least-squares can be solved equally easily with x = bvls(A, b, ub, lb).

What’s coming next?

In the next few weeks, I plan to add a couple of others solvers for standard quadratic programs, including lasso and elastic_net. I also plan to run some benchmarks to get some intuitions about the computational performances of this modernized code compared to other QP solvers (mostly available in Python, Julia or R). If anyone has experience (or just tips) for such benchmarks, I’m all ears!

11 Likes