Modernizing `quadprog` for convex quadratic programming

1

Working on control theory problems, I’ve recently needed to be able to solve convex quadratic programs (QP) of the form

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

where \mathbf{P} is a dense symmetric positive definite matrix of moderate size (say 512 x 512). The matrices \mathbf{A} and \mathbf{C} are also dense. One of the best algorithms for such small-to-moderate scale QP is the Godfarb & Idnani algorithm. A reference implementation is the one by Berwin Turlach (here) which is used under the hood by quadprog in R and other languages.

The original implementation being written in F77, I’ve converted it to F90 and made very small changes here and there to the code to make it completely stand-alone (it does not even need blas, linpack or lapack). You can find the code on this Github repo (still work-in-progress). The legacy solver has a pretty verbose interface that reads like so

interface
      pure module subroutine qpgen2(dmat, dvec, fddmat, n, sol, lagr, crval, amat, bvec, fdamat, q, &
                                    meq, iact, nact, iter, work, ierr)
         integer, intent(in)     :: fddmat, n
      !! Dimensions of the symmetric positive definite matrix Dmat.
         integer, intent(in)     :: fdamat, q
      !! Dimensions of the constraint matrix Amat
         integer, intent(in)     :: meq
      !! Number of equality constraints.
         integer, intent(out)    :: iact(*), nact
      !! Indices and number of active constraints at the optimum.
         integer, intent(out)    :: iter(*)
      !! Number of iterations.
         integer, intent(inout)  :: ierr
      !! Information flag.
         real(dp), intent(inout) :: dmat(fddmat, *), dvec(*)
      !! Sym. pos. def. matrix and vector defining the quadratic cost.
         real(dp), intent(out)   :: lagr(*), sol(*)
      !! Lagrange multipliers and solution vector.
         real(dp), intent(inout) :: amat(fdamat, *), bvec(*)
      !! Matrix and vector defining the (in-)equality constraints.
         real(dp), intent(inout) :: work(*)
      !! Workspace.
         real(dp), intent(out)   :: crval
      !! Cost function at the optimum.
      end subroutine
   end interface

In the process of modernizing the code and making it more user-friendly, I wrapped this qpgen2 subroutine into another function with the following interface

interface
    !!  #### Description
    !!
    !!  ...
    !!
    !!  #### Syntax
    !!
    !!  - To solve an unconstrained quadratic program:
    !!
    !!  ```fortran
    !!  call solve_qp(P, q, x)
    !!  ```
    !!
    !!  #### Arguments
    !!
      module function solve_qp(P, q, Aeq, beq, C, d, y, obj, &
                               factorized, overwrite_p, info) result(x)
         real(dp), intent(inout), contiguous, target  :: P(:, :)
      !! n x n Symmetric positive-definite matrix defining the quadratic form.
         real(dp), intent(in)                         :: q(:)
      !! Vector defining the linear term of the quadratic cost.
         real(dp), optional, intent(in)               :: Aeq(:, :)
      !! Matrix defining the linear equality constraints Aeq @ x = beq.
         real(dp), optional, intent(in)               :: beq(:)
      !! Right-hand side vector of the linear equality constraints Aeq @ x = beq.
         real(dp), optional, intent(in)               :: C(:, :)
      !! Matrix defining the linear inequality constraints C @ x >= d.
         real(dp), optional, intent(in)               :: d(:)
      !! Right-hand side vector of the linear inequality constraints C @ x >= d.
         real(dp), optional, allocatable, target, intent(out) :: y(:)
      !! Vector of Lagrange multipliers at the optimum.
         real(dp), optional, intent(out), target      :: obj
      !! Cost function at the optimum.
         logical, optional, intent(inout)            :: factorized
      !! Whether P has already been factorized using ??? (default .false.)
         logical, optional, intent(in)               :: overwrite_p
      !! Whether P can be overwritten (default .false.)
         integer, intent(out)              :: info
      !! Information flag.

         real(dp), allocatable :: x(:)
      !! Vector corresponding to the minimizer.
      end function
   end interface

which is much more flexible on the user-side and also quite a bit more natural. In my set of unit tests, everything works as expected if I use directly qpgen2, no matter whether I use gfortran or ifx to compile the code with fpm. However, if I use solve_qp, this is where the sh*t hits the fan and I may need help. Solving an unconstrained QP with the following call

x = solve_qp(P, q, obj=obj, info=info)

causes not problem whatsoever, no matter the compiler. However, if I want to solve a constrained problem, and using the following calling sequence

      x = solve_qp(P=P, q=q, Aeq=C, beq=b, y=y, obj=obj, info=info)

only ifx can compile the code and produce the expected output. If I use gfortran, it encounters an internal compile error

 33%] Compiling...
test/problems.f90:158:67:

  158 |       x = solve_qp(P=P, q=q, Aeq=C, beq=b, y=y, obj=obj, info=info)
      |                                                                   1
internal compiler error: Segmentation fault
0x7fc99922551f ???
	./signal/../sysdeps/unix/sysv/linux/x86_64/libc_sigaction.c:0
0x7fc99920cd8f __libc_start_call_main
	../sysdeps/nptl/libc_start_call_main.h:58
0x7fc99920ce3f __libc_start_main_impl
	../csu/libc-start.c:392
Please submit a full bug report, with preprocessed source (by using -freport-bug).
Please include the complete backtrace with any bug report.

I’ve never really had to deal with internal compiler errors so far and so I was wondering if anyone could give me some insights on how to generate a proper bug report for gfortran given that I ain’t sure how to reproduce this issue with a MWE (since I don’t know exactly where it comes from).

2 Likes
2

I suspect the problem coming from the fact that info is a non-optional argument but you placed it at the end of the procedure signature. My advise, try again placing all non-optional first and then all optional arguments.

Non-related to your post: forall is a deprecated feature of the language, I know it is neat for one-liners, but I would advise against its inclusion in new code.

3

Nope, didn’t change anything. That was my first idea as well. As for forall, yep, it was just out of convenience. I eventually plan to make the code as much F2018 compliant as possible with the exception of some go to in qpgen2 which I don’t want to mess with.

4

Humm, I tweaked it a bit here Compiler Explorer it seems to be coming from your choice of splitting the procedure signature within the interface declaration and having a module procedure for the actual implementation. Moving your whole signature below, the ICE goes away.

5

Awesome ! I’ll make the changes on Monday as I get back to the lab’. I’ll also add some derived types to have an even simpler interface.

On a different note, at some point, I think we, as a community, would need to start a discussion about the future of stdlib, what is beyond its scope, etc. Basically, if we think of stdlib as the Fortran equivalent of numpy, what about scipy then? Given this topic, I’m obviously thinking about where an optimization module would fit.

1 Like
6

Yep, I’m personally convinced that stdlib ought to be a complete scientific toolbox for the community. It should cover all those nice tools and algorithms one can find in numpy, scipy, scikit and more.

I wanted to suggest, if you were willing to do this modernization in the scope of stdlib but refrained doing so because it should be your decision :slight_smile: . If you share the same vision, then that’s awesome!!

2 Likes
7

I’m all for stdlib becoming a complete scientific toolbox, given that Fortran is so heavily used in science. I’m currently stretched a little thin, but I’ve been meaning to (and still plan to) dive a little deeper into stdlib, to understand how it works, and to see if/how I could contribute. Would love nothing more than to have a single lib as a toolbox for all regular analyses (i.e. everything I don’t need specialised models for). It would also make it less complicated to re-introduce it to university classrooms for a number of reasons.

4 Likes
8

This warms my heart :grin:

1 Like
9

Removing the go to’s is fairly straightforward but takes a little thought. I made the mods and they seem to work with ifx. The check program said everything worked. I’ll PM my updated solve.QP file.

1 Like
10

Definitely ! For the moment I needed to have rapidly a slightly modernized stand-alone version of this algorithm for a project I’m working one but in the long term (e.g. couple of months) that is definitely something I had in mind. One caveat though is the licence of the original code (LGPL if I recall correctly) which might be incompatible with that of stdlib. I’ll be in the lab in Paris tomorrow with no students, no meeting, no nothing. I’ll take some time to write a longer post reflecting on the past year as an active stdlib user and, from this (admittedly biased) point of view, how I’d like to see the project evolve in the near and longer future.

1 Like
11

Yes, stdlib should be the full toolbox, kind of like scipy.

As to compiler errors, create an MRE (minimal reproducible example) and submit it as a bug report.

1 Like