CppNumericalSolvers: header-only optimisation for C++17

Header-only, one dependency, and a wide Eigen range

CppNumericalSolvers is a header-only C++17 optimization library, and the integration story is unusually short. The only dependency is Eigen3, and the README makes a point of the version range: it is compatible with Eigen 3.4 through Eigen 5.x, and the versions referenced are minimums, not pins. If a project already uses Eigen 5, the library builds against it unchanged.

That sentence matters more than it looks. Eigen has a history of breaking source compatibility across major versions, so a library that pins to a narrow range forces a choice, and one that accepts 3.4 to 5.x removes a whole category of dependency conflict. The cost is that the library has to compile against an API surface that changed across those versions, which is a maintenance burden the author carries.

Header-only means there is nothing to build and nothing to link before the compiler finds it. The include layout follows a namespace, cppoptlib, with the function definitions under cppoptlib/function.h and each solver under cppoptlib/solver/ with a header named after it. The examples in the README include the two headers they need and nothing else.

The other end of the integration story is the build system support, and it is broad. There is CMake integration through FetchContent and through find_package after an install, a pkg-config file so a plain g++ command works, and Bazel support with a MODULE.bazel and a git_override. The repository also carries a .bcr directory, which is the Bazel Central Registry, so the library is consumable from the registry as well as from a checkout, and a WORKSPACE file alongside the module file for older Bazel.

The licence is MIT, with the LICENSE file at the repository root, and the project also ships a CITATION.cff, which is a small sign of a project that expects to be cited in academic work.

The benchmark, and what a self-run benchmark is worth

The performance claim is the first thing in the README, and it is specific. On a 376 problem benchmark against Nocedal's Fortran L-BFGS, libLBFGS, LBFGSpp and LBFGS-Lite, the library has the highest reliability at 95 per cent converged, the most first place wins at 230 out of 376, which the README describes as twice the next library, and the lowest mean nfev of any solver tested, where nfev is the number of function evaluations.

Those are the three numbers that matter for a solver, and they are the right three to report. Reliability first, because a solver that converges more often is more useful than a solver that wins races. Win count second, because a per-problem win rate is more informative than an average. And function evaluations last, because they are the proxy for work done, and they matter for the case where the function evaluation is expensive, which is the usual reason to use a gradient based method in the first place.

The comparison set is sensible. Nocedal's Fortran implementation is the reference everyone cites, libLBFGS is the widely used portable C implementation, and LBFGSpp and LBFGS-Lite are the two common C++ ports. Comparing against those four says something real about the class of solver rather than about an easy target.

The caveat is the obvious one and it should be stated plainly: this is the author's own benchmark, hosted on the author's own pages, which the repository description points at as the homepage. Nobody has independently reproduced those numbers, and a benchmark's design, its problem set and its convergence criteria are all choices the benchmark's author made. The claims are specific enough to be checkable, and the benchmark page is public, so a reader can go and look. But the phrase highest reliability in a README is a marketing sentence until you have read the benchmark.

For a team choosing a solver, the practical takeaway is that this library's L-BFGS implementation is at least competitive with the standard portable implementations, and the honest way to use that information is to run the same problem on your own function before you commit.

Defining an objective, and the differentiability mode that drives everything

An objective is a class that inherits from FunctionCRTP and implements an operator() that takes the point, a pointer to a gradient and, optionally, a pointer to a Hessian. The quick start is the shortest complete example:

#include "cppoptlib/function.h"
#include "cppoptlib/solver/lbfgs.h"

// f(x) = 5*x0^2 + 100*x1^2 + 5
class Quadratic : public cppoptlib::function::FunctionCRTP<
    Quadratic, double, cppoptlib::function::DifferentiabilityMode::First, 2> {
 public:
  ScalarType operator()(const VectorType &x, VectorType *grad) const {
    if (grad) *grad = Eigen::Vector2d(10 * x[0], 200 * x[1]);
    return 5 * x[0] * x[0] + 100 * x[1] * x[1] + 5;
  }
};

The template parameters are the API's central idea. The first is the deriving type, the second is the scalar type, the third is the differentiability mode, and the fourth is optional, a compile time dimension. The differentiability mode is what selects the solver family, and it is the one parameter to get right, because declaring a mode that does not match the derivatives you implement produces wrong answers rather than a compile error.

The gradient pointer is nullable, and that is not an accident. The operator is called in two situations, once to evaluate the function value and once to evaluate value and gradient together, and the same operator serves both, so the pointer is how the caller says which one it wants. The Ridge Regression example in the README shows the same pattern with a Hessian pointer, where the gradient and the Hessian are computed from the same residual.

The dimension template parameter is worth mentioning because it removes a class of error. With a compile time dimension of 2, the gradient is constructed as an Eigen::Vector2d, so a dimension mismatch is a compile error rather than a runtime one. The README also shows the dynamic dimension form, where the class inherits from FunctionCRTP with only the differentiability mode given, which is what you use when the problem size is not known until runtime.

The solve call itself returns a pair, the solution and the state, so the final value and the iteration count are both available without the caller having to reconstruct them.

Nine solvers, and the one that decides your problem shape

The solver table is the map of the library, and it is worth reading as a taxonomy of optimisation problems rather than as a list.

The first order solvers are gradient descent, conjugate gradient, L-BFGS and BFGS. Newton and trust-region Newton are second order, which means they use the Hessian, and Nelder-Mead is zeroth order, meaning it needs no derivatives at all. Two solvers handle constraints: L-BFGS-B handles box constraints, and the augmented Lagrangian method handles equality and inequality constraints with any order.

That last row is the most interesting one, because augmented Lagrangian is not a bound constrained method in disguise. It solves an unconstrained problem with a penalty term and a multiplier, which means it can express equality and inequality constraints, at the cost of tuning. The README's example solves the problem of minimising x0 plus x1 subject to x0 squared plus x1 squared equalling 2, which has the solution of negative one, negative one.

The mechanics of that example show how a constrained problem is assembled. The objective and the constraint are each built as a FunctionExpr, the constraint being the norm expression minus 2, and the two are then handed to a ConstrainedOptimizationProblem:

  cppoptlib::function::FunctionExpr objective = Sum();
  cppoptlib::function::FunctionExpr constraint = cppoptlib::function::FunctionExpr(CircleNorm()) - 2.0;

  cppoptlib::function::ConstrainedOptimizationProblem problem(objective, {constraint});

  cppoptlib::solver::Lbfgs<decltype(problem)::ObjectiveFunctionType> inner;
  cppoptlib::solver::AugmentedLagrangian solver(problem, inner);

  Eigen::VectorXd x0(2);  x0 << 5, -3;
  auto [sol, state] = solver.Minimize(problem,
      cppoptlib::solver::AugmentedLagrangeState<double>(x0, 1, 0, 10.0));

Two things are visible there. The augmented Lagrangian solver wraps an inner unconstrained solver, so the inner solver is a template parameter, and you can put L-BFGS inside it. And the initial state is a distinct type, AugmentedLagrangeState, carrying the starting point along with the penalty and multiplier parameters.

For a team picking a solver, the practical guidance is short. If you have gradients, L-BFGS is the default and the one the benchmark favours. If you have Hessians and the problem is small and well conditioned, Newton converges in a handful of iterations. If you have no derivatives at all, Nelder-Mead is the only option in the table. If the constraints are simple bounds, L-BFGS-B. If they are general, the augmented Lagrangian route, and you accept the extra parameters.

Expression templates for objectives you build from parts

The expression template feature is the one that makes the library pleasant for the kind of objective that appears in machine learning work, where the thing you want to minimise is a sum of a data term and a regularisation term.

The Ridge Regression example is the canonical one. The objective is the squared error of A times x minus y, plus lambda times the squared norm of x, and the README composes it from two separate classes:

  // Compose: F(x) = ||Ax-y||^2 + 0.1 * ||x||^2
  cppoptlib::function::FunctionExpr objective(SquaredError(A, y) + lambda * L2Reg(A.cols()));

  Eigen::VectorXd x0 = Eigen::VectorXd::Zero(A.cols());
  cppoptlib::solver::Lbfgs<decltype(objective)> solver;
  auto [sol, state] = solver.Minimize(objective, cppoptlib::function::FunctionState(x0));

What that buys is the absence of boilerplate. Without expression templates you would write a third class that combines the two, holds both, and forwards the gradient and Hessian calls, and then a fourth when you change the regulariser. With them, the objective is an expression tree built at run time, and the solver is templated on the type of that expression, which is why the solver line says Lbfgs<decltype(objective)>.

The pieces in the example are ordinary classes. SquaredError holds a reference to the matrix A and the vector y, computes the residual A times x minus y, and derives the gradient as 2 A transpose r and the Hessian as 2 A transpose A. L2Reg holds a dimension, derives the gradient as 2x and the Hessian as twice the identity. Both declare DifferentiabilityMode::Second, since they provide Hessians.

One detail is easy to miss and matters in practice: both classes hold references to the matrix and the vector, so the caller is responsible for keeping them alive for the duration of the solve. That is a performance friendly choice, since copying a large matrix for every objective would be wasteful, and it is also a lifetime rule the reader has to honour.

The build story, and the parts of the repository that are not the library

The integration section offers four routes, and the choice depends on your build system rather than on the library.

CMake through FetchContent is the shortest, declaring the repository and tag and making it available, then linking your target against CppNumericalSolvers and Eigen3::Eigen. CMake through find_package after an install is the more conventional route, where the install is a cmake configure with a prefix and a cmake --install, and the consuming project calls find_package for both CppNumericalSolvers and Eigen3 and links the exported target. pkg-config is there for a plain compiler invocation, with the flags obtained by pkg-config and the source compiled with -std=c++17.

Bazel is the fourth, and the details are specific enough to be worth quoting. The module is added with a bazel_dep at a version, and a git_override points at a commit if you want to track master rather than a release. The dependency is then on the include target. The README also explains the Eigen version handling, since cppoptlib declares eigen 3.4.0.bcr.3 as its minimum, and the way to build against Eigen 5 is to request a higher version in your own MODULE.bazel, because bzlmod resolves to the highest requested version.

Outside the library itself, the repository contains a few things that tell you how it is developed. There is a Dockerfile.test, which is how the test suite runs in a controlled environment, a format.sh script, a generator.bzl, and the .bcr directory for the Bazel Central Registry. The release history shows a v2.0.0 tagged as cppoptlib 2.0.0 on 2026-05-11 and a 1.1.0 on 2026-04-27, so the project moved to 2.0 recently, and the last push was on 2026-07-24.

The name of the package is worth flagging for anyone wiring it into a larger build, because the library is CppNumericalSolvers while the namespace, the CMake target and the Bazel module are all cppoptlib, and a mismatch there is the kind of thing that costs an afternoon.

Editorial conclusion

Adopt CppNumericalSolvers when you are solving a nonlinear least squares or general unconstrained problem in C++ and want a solver with no build system friction beyond Eigen, and consider it for constrained problems where the augmented Lagrangian route is acceptable. It is the wrong tool if you need a solver this library does not implement, since the list is nine and does not include a proper SQP or interior point method, and the wrong tool if you want the performance claims to be neutral evidence, because the benchmark is the author's own. Before integrating, confirm the differentiability mode you declare matches the derivatives you actually implement, since the mode selects the solver family and a mismatched mode is a silent source of wrong answers.

Frequently asked questions

Official sources

1. License: MIT
2. PatWie/CppNumericalSolvers on GitHub
3. Project website
4. README
5. Releases