C++/CustomCppProblem/main.cpp

This example shows how to define optimization problems using ordinary C++ functions.

This example shows how to define optimization problems using ordinary C++ functions.

It solves a simple quadratic program of the form:

\[ \begin{aligned} & \underset{x}{\text{minimize}} && \tfrac12 \tp x Q x \\ & \text{subject to} && Ax \le b \\ \end{aligned} \]

Problem specification

The quadratic program described above can be formulated as a general NLP of the following form:

\[ \begin{aligned} & \underset{x}{\text{minimize}} && f(x) \\ & \text{subject to} && x \in C \subseteq \R^n \\ &&& g(x) \in D \subseteq \R^m. \\ \end{aligned} \]

The problem is specified by creating a class (Problem) that inherits from alpaqa's alpaqa::BoxConstrProblem class. It defines the problem-specific functions for the evaluation of the cost function \( f(x) = \tfrac12 \tp xQx \) (eval_f) and its gradient \( \nabla f(x) = Qx \) (eval_grad_f), as well as the constraint function \( g(x) = Ax \) (eval_g) and the function that evaluates the product of the gradient of the constraints and a given vector \( y \), \( \nabla g(x)\,y = \tp A y \) (eval_grad_g_prod).

If you have more efficient ways to combine evaluations of these functions and gradients, you can specify them as well, see the Problem formulations page and the alpaqa::TypeErasedProblem documentation for the full list of supported functions.

The alpaqa::BoxConstrProblem class exposes the two constraint sets, \( x \in C \) and \( g(x) \in D \), and provides the projection functions (eval_prox_grad_step, eval_proj_diff_g and eval_proj_multipliers) for you. The alpaqa::BoxConstrProblem constructor accepts the number of variables \( n = 2 \) and the number of constraints \( m = 1 \) as arguments.

Note

Alpaqa uses structural typing for problem definitions. This means that you just have to provide the supported functions with the correct names and arguments, and the library will pick them up automatically, you don't have to inherit from any abstract interfaces or override any virtual functions. See Problem formulations for more information.

See also

If you haven't already, be sure to go through the Problem formulations page.

Solver selection

The solver consists of three layers:

1. The outer augmented Lagrangian solver (alpaqa::ALMSolver) that handles the general constraints \( g(x) \in D \);
2. The inner PANOC solver (alpaqa::PANOCSolver) that is used to solve the ALM subproblems;
3. The L-BFGS direction (alpaqa::LBFGSDirection) that provides fast Newton-type directions to speed up PANOC.

You can try out different inner solvers (Inner solvers) and direction providers (Direction providers), you can even write your own.

Solver configuration

Each solver and direction class has a set of parameters, which are collected in a struct. Tuning these parameters can often significantly improve solver performance. You can also use them to limit the run time or the number of iterations, and to enable verbose output from the solvers.

Solver invocation

The solver is invoked by calling the solver object with an instance of the problem and an initial guess for the Lagrange multipliers and the decision variables. The solver will overwrite these guesses with the (approximate) solution, and returns a struct containing solver statistics, the most important of which is the status, which reports whether convergence to the desired tolerance was achieved.

#include <alpaqa/panoc-alm.hpp>
#include <alpaqa/problem/box-constr-problem.hpp>
#include <alpaqa/problem/problem-with-counters.hpp>
 
#include <iostream>
 
int main() {
    USING_ALPAQA_CONFIG(alpaqa::DefaultConfig);
 
    // Problem specification
    // minimize  ½ xᵀQx
    //  s.t.     Ax ≤ b
    struct Problem : alpaqa::BoxConstrProblem<config_t> {
        mat Q{n, n};
        mat A{m, n};
        vec b{m};
        mutable vec Qx{n};
 
        Problem() : alpaqa::BoxConstrProblem<config_t>{2, 1} {
            // Initialize problem matrices
            Q << 3, -1, -1, 3;
            A << 2, 1;
            b << -1;
 
            // Specify the bounds
            C.lowerbound = vec::Constant(n, -alpaqa::inf<config_t>);
            C.upperbound = vec::Constant(n, +alpaqa::inf<config_t>);
            D.lowerbound = vec::Constant(m, -alpaqa::inf<config_t>);
            D.upperbound = b;
        }
 
        // Evaluate the cost
        real_t eval_f(crvec x) const {
            Qx.noalias() = Q * x;
            return 0.5 * x.dot(Qx);
        }
        // Evaluat the gradient of the cost
        void eval_grad_f(crvec x, rvec gr) const { gr.noalias() = Q * x; }
        // Evaluate the constraints
        void eval_g(crvec x, rvec g) const { g.noalias() = A * x; }
        // Evaluate a matrix-vector product with the gradient of the constraints
        void eval_grad_g_prod(crvec x, crvec y, rvec gr) const {
            (void)x;
            gr.noalias() = A.transpose() * y;
        }
    };
 
    Problem problem;
 
    // Wrap the problem to count the function evaluations
    auto counted_problem = alpaqa::problem_with_counters_ref(problem);
 
    // Define the solvers to use
    using Direction   = alpaqa::LBFGSDirection<config_t>;
    using InnerSolver = alpaqa::PANOCSolver<Direction>;
    using OuterSolver = alpaqa::ALMSolver<InnerSolver>;
 
    // Settings for the outer augmented Lagrangian method
    OuterSolver::Params almparam;
    almparam.tolerance             = 1e-8; // tolerance
    almparam.dual_tolerance        = 1e-8;
    almparam.penalty_update_factor = 10; // penalty update factor
    almparam.max_iter              = 20;
    almparam.print_interval        = 1;
 
    // Settings for the inner PANOC solver
    InnerSolver::Params panocparam;
    panocparam.max_iter       = 500;
    panocparam.print_interval = 1;
    // Settings for the L-BFGS algorithm used by PANOC
    Direction::AcceleratorParams lbfgsparam;
    lbfgsparam.memory = 2;
 
    // Create an ALM solver using PANOC as inner solver
    OuterSolver solver{
        almparam,                 // params for outer solver
        {panocparam, lbfgsparam}, // inner solver
    };
 
    // Initial guess
    vec x(2);
    x << 2, 2; // decision variables
    vec y(1);
    y << 1; // Lagrange multipliers
 
    // Solve the problem
    auto stats = solver(counted_problem, x, y);
    // y and x have been overwritten by the solution
 
    // Print the results
    std::cout << '\n' << *counted_problem.evaluations << '\n';
    std::cout << "status: " << stats.status << '\n'
              << "f = " << problem.eval_f(x) << '\n'
              << "inner iterations: " << stats.inner.iterations << '\n'
              << "outer iterations: " << stats.outer_iterations << '\n'
              << "ε = " << stats.ε << '\n'
              << "δ = " << stats.δ << '\n'
              << "elapsed time:     "
              << std::chrono::duration<double>{stats.elapsed_time}.count()
              << " s" << '\n'
              << "x = " << x.transpose() << '\n'
              << "y = " << y.transpose() << '\n'
              << "avg τ = " << (stats.inner.sum_τ / stats.inner.count_τ) << '\n'
              << "L-BFGS rejected = " << stats.inner.lbfgs_rejected << '\n'
              << "L-BFGS failures = " << stats.inner.lbfgs_failures << '\n'
              << "Line search failures = " << stats.inner.linesearch_failures
              << '\n'
              << std::endl;
}