# Releasing a modernized \`QuadProg\` for strictly convex quadratic programming

**URL:** <https://fortran-lang.discourse.group/t/releasing-a-modernized-quadprog-for-strictly-convex-quadratic-programming/9587>\
**Category:** Announcements\
**Created:** [April 14, 2025, 1:40pm UTC](https://fortran-lang.discourse.group/t/releasing-a-modernized-quadprog-for-strictly-convex-quadratic-programming/9587 "2025-04-14T13:40:09Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\
**Post date:** [April 14, 2025, 1:40pm UTC](https://fortran-lang.discourse.group/t/releasing-a-modernized-quadprog-for-strictly-convex-quadratic-programming/9587/1 "2025-04-14T13:40:10Z")

</div>

Hej everyone,

I’m happy to announce that the modernized version of [`QuadProg`](https://github.com/loiseaujc/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](https://en.wikipedia.org/wiki/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`](https://www.netlib.org/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}

```fortran
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!