CppNumericalSolvers

Solvers accept a Progress state at construction to determine when Minimize should return. You can use built-in presets or customize them.

1. Default Stopping: DefaultStoppingSolverProgress

Best for well-conditioned problems. It terminates when:

• Gradient norm: $|g|\infty < 1e-5 \times \max(1, |x|\infty)$
• Plateau: The objective moves by less than $1e-6$ (relative to $\max(1, |f|)$) over the last 3 iterations.

auto stop = cppoptlib::solver::DefaultStoppingSolverProgress<Fn, State>();
cppoptlib::solver::Lbfgs<Fn> solver(stop);

2. Conservative Stopping: ConservativeStoppingSolverProgress

Use this when the objective has flat regions (e.g., degenerate saddles or valleys) where the default might mistake a plateau for a minimum. It uses much tighter tolerances:

• Gradient norm: 5e-6
• Plateau window: 5 iterations
• Plateau delta: 1e-10

auto stop = cppoptlib::solver::ConservativeStoppingSolverProgress<Fn, State>();
cppoptlib::solver::Lbfgs<Fn> solver(stop);

3. Customizing via Per-field Overrides

Both presets return a plain struct. You can modify individual fields to create a custom stopping condition without copying the whole preset.

auto stop = cppoptlib::solver::DefaultStoppingSolverProgress<Fn, State>();
stop.num_iterations = 500;
stop.gradient_norm = 1e-7;
cppoptlib::solver::Lbfgs<Fn> solver(stop);

To use CppNumericalSolvers, define your objective function by inheriting from cppoptlib::function::FunctionCRTP. You must specify the scalar type, the differentiability mode (e.g., First or Second), and optionally the dimension. Implement the operator() to return the function value and optionally compute the gradient or Hessian.

Example of a quadratic function $f(x) = 5x_0^2 + 100x_1^2 + 5$:

#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;
  }
};

int main() {
  Quadratic f;
  Eigen::Vector2d x0(-10, 2);
  cppoptlib::solver::Lbfgs<Quadratic> solver;
  auto [solution, state] = solver.Minimize(f, cppoptlib::function::FunctionState(x0));
  // solution.x ≈ (0, 0), solution.value ≈ 5
}

To solve problems with equality or inequality constraints, use the Augmented Lagrangian method. You define the objective and constraints as cppoptlib::function::FunctionExpr objects, wrap them in a cppoptlib::function::ConstrainedOptimizationProblem, and then use an AugmentedLagrangian solver with an inner solver (like Lbfgs).

Example: Minimize $x_0 + x_1$ subject to $x_0^2 + x_1^2 = 2$

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

// ... (Sum and CircleNorm implementations) ...

int main() {
  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));
  // sol.x ≈ (-1, -1), f* ≈ -2
}

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

class Sum : public cppoptlib::function::FunctionXd<Sum> {
 public:
  ScalarType operator()(const VectorType &x, VectorType *g) const {
    if (g) *g = VectorType::Ones(x.size());
    return x.sum();
  }
};

class CircleNorm : public cppoptlib::function::FunctionXd<CircleNorm> {
 public:
  ScalarType operator()(const VectorType &x, VectorType *g) const {
    if (g) *g = 2 * x;
    return x.squaredNorm();
  }
};

int main() {
  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));
  // sol.x ≈ (-1, -1), f* ≈ -2
}

You can build complex objective functions from reusable parts using expression templates. This allows you to combine simple functions (like a squared error term and a regularization term) into a single cppoptlib::function::FunctionExpr without writing boilerplate code for the combined derivatives.

Example: Ridge Regression $F(x) = ||Ax - y||^2 + ext{lambda} imes ||x||^2$

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

class SquaredError : public cppoptlib::function::FunctionCRTP<...
// ... (implementation of SquaredError and L2Reg) ...

int main() {
  Eigen::MatrixXd A(3, 2);  A << 1,2, 3,4, 5,6;
  Eigen::VectorXd y(3);     y << 7, 8, 9;
  double lambda = 0.1;

  // 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));
  std::cout << "x* = " << sol.x.transpose() << ", f* = " << sol.value << "\n";
}

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

class SquaredError : public cppoptlib::function::FunctionCRTP<
    SquaredError, double, cppoptlib::function::DifferentiabilityMode::Second> {
  const Eigen::MatrixXd &A;
  const Eigen::VectorXd &y;
 public:
  SquaredError(const Eigen::MatrixXd &A, const Eigen::VectorXd &y) : A(A), y(y) {}
  int GetDimension() const { return A.cols(); }
  ScalarType operator()(const VectorType &x, VectorType *grad, MatrixType *hess) const {
    Eigen::VectorXd r = A * x - y;
    if (grad) *grad = 2 * A.transpose() * r;
    if (hess) *hess = 2 * A.transpose() * A;
    return r.squaredNorm();
  }
};

class L2Reg : public cppoptlib::function::FunctionCRTP<
    L2Reg, double, cppoptlib::function::DifferentiabilityMode::Second> {
  int dim;
 public:
  explicit L2Reg(int d) : dim(d) {}
  int GetDimension() const { return dim; }
  ScalarType operator()(const VectorType &x, VectorType *grad, MatrixType *hess) const {
    if (grad) *grad = 2 * x;
    if (hess) { hess->setIdentity(dim, dim); *hess *= 2; }
    return x.squaredNorm();
  }
};

int main() {
  Eigen::MatrixXd A(3, 2);  A << 1,2, 3,4, 5,6;
  Eigen::VectorXd y(3);     y << 7, 8, 9;
  double lambda = 0.1;

  // 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));
  std::cout << "x* = " << sol.x.transpose() << ", f* = " << sol.value << "\n";
}