一个轻量级仅头文件的 C++17 库，提供针对（无）约束非线性函数及表达式模板的数值优化方法https://github.com/PatWie/CppNumericalSolvers
CppNumericalSolvers 库 include 目录下的文件及其功能说明
根目录文件
文件名 | 功能说明 |
|---|---|
function.h | (主函数头文件) 定义目标函数的核心头文件。它通常包含了其他 |
function_base.h | (函数基类) 定义了所有目标函数的基类接口，如 |
function_expressions.h | (函数表达式) 实现了强大的函数表达式模板功能。它允许你通过运算符重载（如 |
function_penalty.h | (罚函数) 提供了用于约束优化的罚函数（Penalty Method）相关实现。这是一种将约束问题转化为无约束问题的方法。 |
function_problem.h | (优化问题定义) 定义了一个完整的优化问题结构，特别是对于有约束问题，它将目标函数和一系列约束函数（等式和不等式）封装在一起。 |
linesearch 目录 (线搜索算法)
这个目录包含了用于在一维方向上寻找最优步长的算法。
文件名 | 功能说明 |
|---|---|
armijo.h | (Armijo 线搜索) 实现了 Armijo 回溯线搜索算法，一种简单而常用的非精确线搜索方法，用于确保每一步迭代都能充分减小目标函数值。 |
more_thuente.h | (Moré-Thuente 线搜索) 实现了 Moré-Thuente 线搜索算法，这是一种更复杂的线搜索方法，试图同时满足 Wolfe 条件，通常能提供更好的收敛性。 |
solver 目录 (求解器)
这个目录包含了库中实现的核心优化算法。
文件名 | 功能说明 |
|---|---|
augmented_lagrangian.h | (增广拉格朗日法) 实现了增广拉格朗日求解器，这是一种强大的处理带等式和不等式约束优化问题的方法。 |
bfgs.h | (BFGS 求解器) 实现了 BFGS (Broyden–Fletcher–Goldfarb–Shanno) 拟牛顿法求解器，一种非常流行和高效的无约束优化算法。 |
conjugated_gradient_descent.h | (共轭梯度下降法) 实现了共轭梯度法求解器，适用于大规模无约束优化问题，因为它不需要存储海森矩阵。 |
gradient_descent.h | (梯度下降法) 实现了基本的梯度下降法求解器，沿着梯度的反方向进行迭代。 |
lbfgs.h | (L-BFGS 求解器) 实现了 L-BFGS (有限内存 BFGS) 求解器，是 BFGS 的内存优化版本，特别适用于变量维度非常高的问题。 |
lbfgsb.h | (L-BFGS-B 求解器) 实现了 L-BFGS-B 求解器，用于处理带箱形约束（Bound-constrained）的优化问题，即变量有上下界限制。 |
nelder_mead.h | (Nelder-Mead 求解器) 实现了 Nelder-Mead 单纯形法，这是一种直接搜索方法，不需要计算梯度，适用于目标函数不可微或计算梯度困难的情况。 |
newton_descent.h | (牛顿法) 实现了牛顿下降法求解器，它使用目标函数的二阶导数（海森矩阵），收敛速度快但计算成本较高。 |
progress.h | (进度监控) 提供了用于监控优化过程的工具，例如 README 中提到的 |
solver.h | (求解器基类) 定义了所有求解器的通用接口和基类，规定了求解器的基本行为，如 |
utils 目录 (工具)
这个目录包含了一些辅助功能和工具。
文件名 | 功能说明 |
|---|---|
derivatives.h | (导数工具) 提供了数值计算导数的工具，主要通过有限差分法来近似计算梯度和海森矩阵。这对于验证用户提供的解析导数的正确性非常有用。 |
示例代码解析
1. constrained_simple.cc
这份代码是一个使用 CppNumericalSolvers 库解决带有一个等式约束和一个不等式约束的优化问题的完整示例。它演示了如何定义目标函数、约束函数，如何将它们组合成一个约束问题，以及如何使用 AugmentedLagrangian (增广拉格朗日) 求解器来找到最优解。
// Copyright 2025, https://github.com/PatWie/CppNumericalSolvers
// 版权所有 2025, https://github.com/PatWie/CppNumericalSolvers
#include<iostream>// 引入标准输入输出流库，用于打印结果#include"cppoptlib/function.h"// 引入库中定义函数所需的核心头文件
#include"cppoptlib/solver/augmented_lagrangian.h"// 引入增广拉格朗日求解器
#include"cppoptlib/solver/lbfgs.h"// 引入 L-BFGS 求解器，用作增广拉格朗日法的内层求解器// Define a 2D function type with first-order differentiability.
// 定义一个二维、一阶可微的函数类型别名。
// 这简化了后面定义函数类的代码。
template<classF>
using Function2d = cppoptlib::function::FunctionCRTP<F,double, cppoptlib::function::DifferentiabilityMode::First,2>;// Alias for an "any function" with the above properties.
// 为具有上述属性的“任意函数”（通过类型擦除）创建别名。
// 这在将具体函数存入通用容器或作为参数传递时很有用。
using FunctionExpr2d = cppoptlib::function::FunctionExpr<
double, cppoptlib::function::DifferentiabilityMode::First,2>;//
// QuadraticObjective2: f(x) = (x[0]-1)^2 + (x[1]-2)^2
// 二次目标函数2：f(x) = (x[0]-1)^2 + (x[1]-2)^2
//
// The unconstrained optimum is (1,2) with f(x)=0.
// 无约束最优解是 (1,2)，此时 f(x)=0。
// However, the constraints (below) force a different solution.
// 然而，下面的约束会强制得到一个不同的解。
//
classQuadraticObjective2:publicFunction2d<QuadraticObjective2>{
public:EIGEN_MAKE_ALIGNED_OPERATOR_NEW // Eigen 内存对齐宏// 重载 operator()，实现函数求值和梯度计算ScalarType operator()(const VectorType &x,VectorType *gradient =nullptr)const{
if(gradient){// 如果 gradient 指针非空，则计算梯度
// Gradient: 2*(x - [1,2])
// 梯度公式：2*(x - [1,2])VectorType ref(2);ref <<1,2;
*gradient =2*(x - ref);
}
// 返回函数值
return(x(0)-1)*(x(0)-1)+(x(1)-2)*(x(1)-2);
}
};//
// EqualityConstraint2: g(x) = x[0] - 0.5 = 0
// 等式约束2：g(x) = x[0] - 0.5 = 0
//
// This forces x[0] to be exactly 0.5.
// 这强制 x[0] 必须等于 0.5。
// Thus, the optimal solution must have x[0] = 0.5.
// 因此，最优解必须满足 x[0] = 0.5。
//
classEqualityConstraint2:publicFunction2d<EqualityConstraint2>{
public:EIGEN_MAKE_ALIGNED_OPERATOR_NEW// 实现约束函数的求值和梯度计算ScalarType operator()(const VectorType &x,VectorType *gradient =nullptr)const{
if(gradient){// 如果请求梯度
// Gradient is [1, 0].
// 梯度是 [1, 0]。
*gradient =(VectorType(2)<<1,0).finished();
}
// 返回约束函数的值
returnx(0)-0.5;
}
};//
// InequalityConstraint3: h(x) = 2 - (x[0] + x[1]) >= 0
// 不等式约束3：h(x) = 2 - (x[0] + x[1]) >= 0
//
// This requires x[0] + x[1] <= 2. With x[0]=0.5 from the equality constraint,
// we have x[1] <= 1.5.
// 这要求 x[0] + x[1] <= 2。结合等式约束 x[0]=0.5，我们得到 x[1] <= 1.5。
//
classInequalityConstraint3:publicFunction2d<InequalityConstraint3>{
public:EIGEN_MAKE_ALIGNED_OPERATOR_NEW// 实现约束函数的求值和梯度计算ScalarType operator()(const VectorType &x,VectorType *gradient =nullptr)const{
if(gradient){// 如果请求梯度
// Gradient is [-1, -1].
// 梯度是 [-1, -1]。
*gradient =(VectorType(2)<<-1,-1).finished();
}
// 返回约束函数的值
return2-(x(0)+x(1));
}
};//
// Main demo: Solve the constrained problem with multiple constraints
// 主演示程序：解决带有多个约束的约束问题
//
intmain(){
// Initial guess for x. We choose a starting point that is not optimal.
// x 的初始猜测值。我们选择一个非最优的起始点。QuadraticObjective2::VectorType x(2);x <<1,1;// Define the objective function.
// 定义目标函数。
// 将具体的目标函数类包装成通用的 FunctionExpr 类型。cppoptlib::function::FunctionExpr objective =QuadraticObjective2();// Define the constraint functions.
// 定义约束函数。
// 同样将具体的约束函数类包装成通用的 FunctionExpr 类型。cppoptlib::function::FunctionExpr eq =EqualityConstraint2();cppoptlib::function::FunctionExpr ineq =InequalityConstraint3();// Construct the constrained optimization problem.
// 构造约束优化问题。
//
// Problem Statement:
// 问题描述：
// minimize f(x) = (x[0]-1)^2 + (x[1]-2)^2
// 最小化 f(x) = (x[0]-1)^2 + (x[1]-2)^2
// subject to:
// 约束条件：
// equality constraint: x[0] - 0.5 == 0 (forces x[0]=0.5)
// 等式约束：x[0] - 0.5 == 0 (强制 x[0]=0.5)
// inequality constraint: 2 - (x[0]+x[1]) >= 0 (forces x[0]+x[1] <= 2)
// 不等式约束：2 - (x[0]+x[1]) >= 0 (强制 x[0]+x[1] <= 2)
//
// Expected Outcome:
// 预期结果：
// The unconstrained optimum is (1,2) with f(x)=0, but it is infeasible
// because 1+2=3 > 2. With x[0] forced to 0.5, the inequality requires x[1]
// <= 1.5. The best feasible choice is x = (0.5, 1.5), giving:
// 无约束最优解是 (1,2)，f(x)=0，但这个解是不可行的，
// 因为 1+2=3 > 2。在 x[0] 被强制为 0.5 的情况下，不等式要求 x[1] <= 1.5。
// 最佳的可行选择是 x = (0.5, 1.5)，此时：
// f(x) = (0.5-1)^2 + (1.5-2)^2 = 0.25 + 0.25 = 0.5.
//
// 将目标函数、等式约束列表和不等式约束列表打包成一个 ConstrainedOptimizationProblem 对象。cppoptlib::function::ConstrainedOptimizationProblem prob(objective,
/* equality constraints */{eq},// 等式约束列表
/* inequality constraints */{ineq});// 不等式约束列表// Set up an inner LBFGS solver for unconstrained subproblems.
// 设置一个内层的 L-BFGS 求解器，用于解决无约束子问题。cppoptlib::solver::Lbfgs<FunctionExpr2d> unconstrained_solver;// Create the augmented Lagrangian solver that handles the constraints.
// 创建处理约束的增广拉格朗日求解器。
// 它接收约束问题 `prob` 和一个无约束求解器 `unconstrained_solver`作为参数。cppoptlib::solver::AugmentedLagrangian solver(prob, unconstrained_solver);// Initialize the augmented Lagrange state.
// 初始化增广拉格朗日的状态。
// 参数：初始点x，等式约束数量(1)，不等式约束数量(1)，初始惩罚因子(1.0)cppoptlib::solver::AugmentedLagrangeState<double,2>l_state(x,1,1,1.0);// Run the solver.
// 运行求解器。
// `Minimize` 函数返回一个包含最终解和求解器状态的元组。
auto[solution, solver_state]= solver.Minimize(prob, l_state);// Output the results.
// 输出结果。std::cout <<"Optimal f(x): "<<objective(solution.x)<< std::endl;std::cout <<"Optimal x: "<< solution.x.transpose()<< std::endl;std::cout <<"Iterations: "<< solver_state.num_iterations << std::endl;std::cout <<"Solver status: "<< solver_state.status << std::endl;// Expected Output:
// 预期输出：
// Optimal x should be close to (0.5, 1.5).
// 最优解 x 应该接近 (0.5, 1.5)。
// The optimal function value f(x) should be approximately 0.5.
// 最优函数值 f(x) 应该约等于 0.5。
// Both constraints are active: x[0]=0.5 (equality) and x[0]+x[1]=2
// (inequality).
// 两个约束都是激活的：x[0]=0.5 (等式) 和 x[0]+x[1]=2 (不等式)。
return0;
}2. constrained_simple2.cc
这份代码是使用 CppNumericalSolvers 库解决另一个约束优化问题的示例。这个例子非常经典：在一个圆形区域上，寻找一个点，使得该点坐标的和最小。它巧妙地展示了如何使用同一个基础函数（Circle）来构造等式和不等式约束。
// Copyright 2025, https://github.com/PatWie/CppNumericalSolvers
// 版权所有 2025, https://github.com/PatWie/CppNumericalSolvers
#include<iostream>// 引入标准输入输出流库，用于打印结果#include"cppoptlib/function.h"// 引入库中定义函数所需的核心头文件
#include"cppoptlib/solver/augmented_lagrangian.h"// 引入增广拉格朗日求解器
#include"cppoptlib/solver/lbfgs.h"// 引入 L-BFGS 求解器//
// SumObjective: f(x) = x[0] + x[1] (to be minimized)
// 求和目标函数：f(x) = x[0] + x[1] (待最小化)
//
classSumObjective:public cppoptlib::function::FunctionXd<SumObjective>{
public:EIGEN_MAKE_ALIGNED_OPERATOR_NEW // Eigen 内存对齐宏// Return the sum of the components.
// 返回各分量之和。ScalarType operator()(const VectorType &x,VectorType *gradient =nullptr)const{
if(gradient){// 如果 gradient 指针非空，则计算梯度
// The gradient is a vector of ones.
// 梯度是一个全为 1 的向量。
*gradient =VectorType::Ones(2);
}
// 返回 x 各分量之和
return x.sum();
}
};//
// Circle: c(x) = x[0]^2 + x[1]^2
// 圆函数：c(x) = x[0]^2 + x[1]^2
// This function will be used to form both an equality constraint (forcing
// the solution to lie on the circle) and an inequality constraint (ensuring
// the solution remains within the circle).
// 这个函数将被用来构造一个等式约束（强制解落在圆上）和一个不等式约束（确保解保持在圆内）。
//
classCircle:public cppoptlib::function::FunctionXd<Circle>{
public:EIGEN_MAKE_ALIGNED_OPERATOR_NEW// Compute the squared norm.
// 计算平方范数。ScalarType operator()(const VectorType &x,VectorType *gradient =nullptr)const{
if(gradient){// 如果请求梯度
// 梯度是 2*x
*gradient =2* x;
}
// 返回 x 的平方范数 (x[0]^2 + x[1]^2)
return x.squaredNorm();
}
};//
// Main demo: solve the constrained problem
// 主演示程序：解决约束问题
//
intmain(){
// Initial guess for x.
// x 的初始猜测值。SumObjective::VectorType x(2);x <<2,10;// Define the objective: f(x) = x[0] + x[1].
// 定义目标函数：f(x) = x[0] + x[1]。
// 将具体的 SumObjective 类包装成通用的 FunctionExpr 类型。cppoptlib::function::FunctionExpr objective =SumObjective();// Define the circle function: c(x) = x[0]^2 + x[1]^2.
// 定义圆函数：c(x) = x[0]^2 + x[1]^2。
// 将具体的 Circle 类包装成通用的 FunctionExpr 类型。cppoptlib::function::FunctionExpr circle =Circle();// Build the constrained optimization problem.
// 构造约束优化问题。
// We impose two constraints:
// 我们施加两个约束：
// 1. Equality constraint: circle(x) - 2 == 0, forcing the solution onto the
// circle's boundary.
// 1. 等式约束: circle(x) - 2 == 0，强制解落在圆的边界上。
// 2. Inequality constraint: 2 - circle(x) >= 0, ensuring the solution remains
// inside the circle.
// 2. 不等式约束: 2 - circle(x) >= 0，确保解保持在圆的内部。
//
// 这里使用了表达式模板来动态地创建约束函数。
// `circle - 2` 创建了一个新的函数表达式，其值为 circle(x) - 2。
// `2 - circle` 创建了另一个新的函数表达式，其值为 2 - circle(x)。cppoptlib::function::ConstrainedOptimizationProblem prob(objective,
/* equality constraints */
{cppoptlib::function::FunctionExpr(circle -2)},// 等式约束列表
/* inequality constraints */
{cppoptlib::function::FunctionExpr(2- circle)});// 不等式约束列表
// 注意：这个例子中的两个约束实际上是冗余的。等式约束已经包含了不等式约束。
// 这里主要是为了演示如何同时处理两种类型的约束。// Set up an inner solver (LBFGS) for the unconstrained subproblems.
// 设置一个内层的 L-BFGS 求解器，用于解决无约束子问题。cppoptlib::solver::Lbfgs<cppoptlib::function::FunctionExprXd> inner_solver;// Create the augmented Lagrangian solver using the inner solver.
// 使用内层求解器创建增广拉格朗日求解器。cppoptlib::solver::AugmentedLagrangian solver(prob, inner_solver);// Initialize the augmented Lagrange state.
// 初始化增广拉格朗日的状态。
// 参数：初始点x，等式约束数量(1)，不等式约束数量(1)，初始惩罚因子(1.0)cppoptlib::solver::AugmentedLagrangeState<double>l_state(x,1,1,1.0);// Run the solver.
// 运行求解器。
auto[solution, solver_state]= solver.Minimize(prob, l_state);// Output the results.
// 输出结果。std::cout <<"Optimal f(x): "<<objective(solution.x)<< std::endl;std::cout <<"Optimal x: "<< solution.x.transpose()<< std::endl;std::cout <<"Iterations: "<< solver_state.num_iterations << std::endl;std::cout <<"Solver status: "<< solver_state.status << std::endl;// 预期结果：
// 为了最小化 x1+x2，同时满足 x1^2+x2^2=2，
// 几何上是在圆上找一条 y=-x+c 的直线，使其与圆相切且截距 c 最小。
// 这发生在圆的第三象限，即 x1=x2 的点。
// x1^2 + x1^2 = 2 => 2*x1^2 = 2 => x1 = -1。
// 所以最优解是 (-1, -1)，最优值是 -2。return0;
}3. debug.cc
这份代码是一份全面的单元测试和功能演示，旨在展示 CppNumericalSolvers 库中与函数定义、函数表达式和约束问题构建相关的核心功能。它不进行任何实际的优化求解，而是侧重于验证以下几点：
如何定义不同可微性等级（零阶、一阶、二阶）的函数。
如何使用
FunctionExpr类型擦除包装器来统一处理这些不同类型的函数。如何使用重载的运算符（
+,-,*）通过表达式模板来动态构建复杂的函数。如何使用库提供的工具来构建约束问题的惩罚项和增广拉格朗日函数。
#include<Eigen/Dense>// 引入 Eigen 密集矩阵和向量库
#include<cmath>// 引入 C++ 标准数学库
#include<iostream>// 引入标准输入输出流库
#include<memory>// 引入智能指针
#include<stdexcept>// 引入标准异常
#include<type_traits>// 引入类型萃取工具#include"cppoptlib/function.h"// 引入库中定义函数所需的核心头文件usingnamespace cppoptlib::function;// 使用 cppoptlib::function 命名空间// LinearFunction: supports first-order information.
// 线性函数：支持一阶信息（梯度）。
// 继承自 FunctionXd，这是一个预定义的别名，代表一阶可微、double类型的函数。
structLinearFunction:publicFunctionXd<LinearFunction>{Eigen::VectorXd m;// 斜率向量
double b;// 截距LinearFunction(const Eigen::VectorXd &m_,double b_):m(m_),b(b_){}// Simplified operator(): only function value and gradient.
// 简化的 operator()：只提供函数值和梯度。
doubleoperator()(const VectorType &x, VectorType *grad =nullptr)const{
double fx = m.dot(x)+ b;// f(x) = m^T * x + b
if(grad){// 如果请求梯度
*grad = m;// 梯度就是 m
}
return fx;
}
};// ConstantFunction: supports no derivative information.
// 常数函数：不支持导数信息。
// 直接继承自 FunctionCRTP，并指定可微性模式为 None。
structConstantFunction
:public cppoptlib::function::FunctionCRTP<ConstantFunction,double,DifferentiabilityMode::None>{
double c;// 常数值
ConstantFunction(double c_):c(c_){}// 最简化的 operator()：只接受 x，返回函数值。
doubleoperator()(const VectorType &x)const{
(void)x;// 明确表示 x 未被使用，避免编译器警告。
return c;
}
};// QuadraticFunction: supports second-order information.
// 二次函数：支持二阶信息（梯度和海森矩阵）。
// 继承自 FunctionCRTP，并指定可微性模式为 Second。
structQuadraticFunction
:public cppoptlib::function::FunctionCRTP<QuadraticFunction,double,DifferentiabilityMode::Second>{Eigen::MatrixXd A;// 二次项矩阵Eigen::VectorXd b;// 线性项向量
double c;// 常数项QuadraticFunction(const Eigen::MatrixXd &A_,const Eigen::VectorXd &b_,
double c_)
:A(A_),b(b_),c(c_){}// 完整的 operator()：提供函数值、梯度和海森矩阵。
doubleoperator()(const VectorType &x, VectorType *grad =nullptr,MatrixType *hess =nullptr)const{
double fx =0.5* x.transpose()* A * x + b.dot(x)+ c;// f(x) = 0.5*x^T*A*x + b^T*x + c
if(grad){// 如果请求梯度
*grad = A * x + b;// 梯度是 A*x + b
}
if(hess){// 如果请求海森矩阵
*hess = A;// 海森矩阵就是 A
}
return fx;
}
};
// 二次函数2，与上面类似，但使用静态维度 (2D)
structQuadraticFunction2
:public cppoptlib::function::FunctionCRTP<
QuadraticFunction2,double, DifferentiabilityMode::Second,2>{Eigen::Matrix2d A;// 使用固定大小的 Eigen 类型Eigen::Vector2d b;
double c;QuadraticFunction2(const Eigen::Matrix2d &A_,const Eigen::Vector2d &b_,
double c_)
:A(A_),b(b_),c(c_){}doubleoperator()(const VectorType &x, VectorType *grad =nullptr,MatrixType *hess =nullptr)const{
double fx =0.5* x.transpose()* A * x + b.dot(x)+ c;
if(grad){
*grad = A * x + b;
}
if(hess){
*hess = A;
}
return fx;
}
};//-----------------------------------------------------------------
// Example Usage
// 示例用法
intmain(){
using VectorType = Eigen::VectorXd;// 使用动态大小的向量类型
using MatrixType = Eigen::MatrixXd;// 使用动态大小的矩阵类型// Test LinearFunction via FunctionExpr (First-order).
// 通过 FunctionExpr 测试线性函数 (一阶)。
{VectorType m(2);m <<2,-1;
double b_val =0.5;LinearFunction lin(m, b_val);// 创建一个具体的线性函数对象FunctionExpr anyLin(lin);// 将其包装进类型擦除的 FunctionExpr 中VectorType x(2);x <<1,2;VectorType grad(2);
double f_lin =anyLin(x,&grad);// 调用包装器，就像调用原始函数一样std::cout <<"Linear function value: "<< f_lin <<"\n";std::cout <<"Gradient: "<< grad.transpose()<<"\n\n";
}// Test ConstantFunction via FunctionExpr (None).
// 通过 FunctionExpr 测试常数函数 (零阶)。
{ConstantFunction cf(3.14);// 创建一个常数函数对象FunctionExpr anyConst(cf);// 包装它VectorType x(2);x <<1,2;
double f_const =anyConst(x);// 调用std::cout <<"Constant function value: "<< f_const <<"\n\n";// 测试表达式模板：对一个 FunctionExpr 对象取负FunctionExpr negAnyConst =-anyConst;f_const =negAnyConst(x);std::cout <<"Constant function value: "<< f_const <<"\n\n";
}// Test QuadraticFunction via FunctionExpr (Second-order).
// 通过 FunctionExpr 测试二次函数 (二阶)。
{MatrixType A(2,2);A <<3,1,1,2;VectorType b(2);b <<1,-1;
double c_val =0.5;QuadraticFunction quad(A, b, c_val);// 创建二次函数对象FunctionExpr anyQuad(quad);// 包装VectorType x(2);x <<1,2;VectorType grad(2);MatrixType hess(2,2);
double f_quad =anyQuad(x,&grad,&hess);// 调用，获取值、梯度和海森矩阵std::cout <<"Quadratic function value: "<< f_quad <<"\n";std::cout <<"Gradient: "<< grad.transpose()<<"\n";std::cout <<"Hessian:\n"<< hess <<"\n\n";
}
// Test QuadraticFunction2 via FunctionExpr (Second-order).
// 测试静态维度的二次函数。
{
using VectorType = Eigen::Vector2d;
using MatrixType = Eigen::Matrix2d;MatrixType A;A <<3,1,1,2;VectorType b;b <<1,-1;
double c_val =0.5;QuadraticFunction2 quad(A, b, c_val);FunctionExpr anyQuad(quad);VectorType x(2);x <<1,2;VectorType grad(2);MatrixType hess(2,2);
double f_quad =anyQuad(x,&grad,&hess);std::cout <<"Quadratic function value: "<< f_quad <<"\n";std::cout <<"Gradient: "<< grad.transpose()<<"\n";std::cout <<"Hessian:\n"<< hess <<"\n\n";
}// Test expression templates:
// 测试表达式模板：
// Expression: (quad + lin - lin * 2.0)
// 表达式: (quad + lin - lin * 2.0)
{MatrixType A(2,2);A <<3,1,1,2;VectorType b(2);b <<1,-1;
double c_val =0.5;QuadraticFunction quad(A, b, c_val);VectorType m(2);m <<2,-1;
double b_val =0.5;LinearFunction lin(m, b_val);// Expression: quad + lin (Min(Differentiability(quad),
// Differentiability(lin)) = First)
// 表达式: quad + lin (最终可微性取两者中的较低者，即一阶)
auto expr1 = quad + lin;
// Expression: lin * 2.0
// 表达式: lin * 2.0
auto expr2 = lin *2.0;
// Expression: (quad + lin) - (lin * 2.0)
// 表达式: (quad + lin) - (lin * 2.0)
auto expr = expr1 - expr2;FunctionExpr anyExpr(expr);// 将最终的复杂表达式包装起来VectorType x(2);x <<1,2;VectorType grad(2);
double f_expr =anyExpr(x,&grad);// 调用复杂表达式std::cout <<"Expression (quad + lin - lin*2.0) value: "<< f_expr <<"\n";std::cout <<"Expression gradient: "<< grad.transpose()<<"\n";// 测试两个函数相乘的表达式FunctionExpr prodFunc(quad * quad);VectorType grad2(2);
quad(x,&grad2);std::cout <<"Expression gradient: "<< grad2.transpose()<<"\n";// 这行打印的是 quad 的梯度
double f2_expr =prodFunc(x,&grad2);std::cout <<"Expression (quad * quad) value: "<< f2_expr <<"\n";std::cout <<"Expression gradient: "<< grad2.transpose()<<"\n";// 这行打印的是 (quad*quad) 的梯度
}// --- Equality Constraint Example ---
// --- 等式约束示例 ---
// Constraint: 2*x0 + x1 - 3 == 0
// 约束: 2*x0 + x1 - 3 == 0
// Penalty function: 0.5 * [2*x0 + x1 - 3]^2
// 惩罚函数: 0.5 * [2*x0 + x1 - 3]^2
{
// Define a linear function: f(x) = 2*x0 + 1*x1 - 3.
// 定义一个线性函数: f(x) = 2*x0 + 1*x1 - 3。VectorType m(2);m <<2,1;
double b_val =-3;
// 使用库函数将线性约束 c(x)=0 转化为二次惩罚函数 P(x)=0.5*c(x)^2FunctionExpr anyEqPenalty =
quadraticEqualityPenalty(LinearFunction(m, b_val));// Evaluate at x = (1, 2)
// 在点 x = (1, 2) 处求值VectorType x(2);x <<1,2;VectorType grad(2);
double eqPenaltyValue =anyEqPenalty(x,&grad);// 调用惩罚函数std::cout <<"Equality constraint penalty (2*x0 + x1 - 3 == 0) at x=(1,2): "
<< eqPenaltyValue <<"\n";std::cout <<"Gradient: "<< grad.transpose()<<"\n\n";
}// --- 构造增广拉格朗日函数示例 ---
{
using TScalar =double;
using VectorType = Eigen::VectorXd;
using MatrixType = Eigen::MatrixXd;// Define the objective function as a quadratic:
// 将目标函数定义为二次函数：MatrixType A(2,2);A <<3,1,1,2;VectorType b(2);b <<1,-1;TScalar c_val =0.5;FunctionExpr objective =QuadraticFunction(A, b, c_val);// Define two constraints.
// 定义两个约束。
// Constraint 0: Equality constraint: 2*x0 + x1 - 3 == 0.
// 约束 0: 等式约束: 2*x0 + x1 - 3 == 0。VectorType m1(2);m1 <<2,1;TScalar b1 =-3;FunctionExpr eqConstraint =LinearFunction(m1, b1);// Constraint 1: Inequality constraint: x0 - 1 <= 0.
// 约束 1: 不等式约束: x0 - 1 <= 0。
// (We interpret "≤" as an InequalityLe constraint.)
// (我们将其解释为“小于等于”的不等式约束。)VectorType m2(2);m2 <<1,0;TScalar b2 =-1;FunctionExpr ineqConstraint =LinearFunction(m2, b2);// Build the optimization problem.
// 构造优化问题。
// 注意：这里的约束格式与库的标准 c(x)>=0 不完全一致，但主要用于演示构建过程。
// eqConstraint 表示 c(x)=0，-eqConstraint 也表示 c(x)=0。ConstrainedOptimizationProblem optimization_problem(objective,{FunctionExpr(-eqConstraint)},{ineqConstraint});// 创建拉格朗日乘子和惩罚因子的状态LagrangeMultiplierState<double>l_state(1,1);l_state.equality_multipliers[0]=1;l_state.inequality_multipliers[0]=1;l_state =LagrangeMultiplierState<double>({1},{1});// 另一种初始化方式PenaltyState<double>p_state(1);// 使用库函数构建完整的惩罚函数和增广拉格朗日函数FunctionExpr unconstrainedPenalty =
ToPenalty(optimization_problem, p_state);FunctionExpr unconstrainedAugmentedLagrangian =
ToAugmentedLagrangian(optimization_problem, l_state, p_state);VectorType x(2);x <<1,3;// 分别构建拉格朗日部分和惩罚部分FunctionExpr penalty_part =FormPenaltyPart(optimization_problem, p_state);FunctionExpr l_part =FormLagrangianPart(optimization_problem, l_state);// Evaluate at the initial state.
// 在初始点处求值，以验证构建是否正确。VectorType grad(2);TScalar objValue;objValue =eqConstraint(x,&grad);std::cout <<"eqConstraint at x0: "<< objValue <<"\n";objValue =ineqConstraint(x,&grad);std::cout <<"ineqConstraint at x0: "<< objValue <<"\n";objValue =unconstrainedAugmentedLagrangian(x,&grad);std::cout <<"Augmented Lagrangian value at x0: "<< objValue <<"\n";objValue =l_part(x,&grad);std::cout <<"Lpart value at x0: "<< objValue <<"\n";objValue =penalty_part(x,&grad);std::cout <<"Penalty value at x0: "<< objValue <<"\n";objValue =unconstrainedPenalty(x,&grad);std::cout <<"obj + penalty value at x0: "<< objValue <<"\n";std::cout <<"Gradient: "<< grad.transpose()<<"\n";
}return0;
}4. linear_regression.cc
这份代码非常巧妙，它用两种完全不同的方法来解决同一个带箱形约束的线性回归问题，并对比它们的结果，从而展示了 CppNumericalSolvers 库的灵活性。
- 方法一 (L-BFGS-B)
: 直接使用专门为箱形约束设计的
Lbfgsb求解器。 - 方法二 (Augmented Lagrangian)
: 将箱形约束显式地表达为四个独立的不等式约束，然后使用通用的
AugmentedLagrangian求解器来解决。
这个问题本身是一个简单的最小二乘问题，目标是找到参数 beta1 和 beta2 (在代码中是 x[0] 和 x[1])，使得模型能最好地拟合两个数据点。
// Copyright 2025, https://github.com/PatWie/CppNumericalSolvers
// 版权所有 2025, https://github.com/PatWie/CppNumericalSolvers#include<iostream>// 引入标准输入输出流库
#include<limits>// 引入数值极限#include"Eigen/Core"// 引入 Eigen 核心库
#include"cppoptlib/function.h"// 引入库中定义函数所需的核心头文件
#include"cppoptlib/solver/augmented_lagrangian.h"// 引入增广拉格朗日求解器
#include"cppoptlib/solver/lbfgs.h"// 引入 L-BFGS 求解器
#include"cppoptlib/solver/lbfgsb.h"// 引入 L-BFGS-B 求解器usingnamespace cppoptlib::function;// 使用 cppoptlib::function 命名空间// 定义一个线性回归的目标函数（最小二乘法）
// 继承自 FunctionXf，这是一个预定义的别名，代表一阶可微、float类型的函数。
classLinearRegression:publicFunctionXf<LinearRegression>{
public:EIGEN_MAKE_ALIGNED_OPERATOR_NEW // Eigen 内存对齐宏// 重载 operator()，实现函数求值和梯度计算ScalarType operator()(const VectorType &x,VectorType *gradient =nullptr)const{
// Compute residuals for the two observations:
// 计算两个观测值的残差：
// Observation 1: beta1 + 2*beta2 - 4
// 观测点 1: x[0] + 2*x[1] - 4
// Observation 2: 3*beta1 + beta2 - 5
// 观测点 2: 3*x[0] + x[1] - 5ScalarType r1 = x[0]+2* x[1]-4;ScalarType r2 =3* x[0]+ x[1]-5;// Compute the objective function: sum of squared errors
// 计算目标函数：误差平方和ScalarType f = r1 * r1 + r2 * r2;// Compute the gradient if requested.
// 如果请求梯度，则进行计算。
// The gradient is: grad f = 2 * [ r1 + 3*r2, 2*r1 + r2 ]
// 梯度公式：grad f = 2 * [ r1 + 3*r2, 2*r1 + r2 ]
if(gradient){gradient->resize(x.size());
(*gradient)[0]=2*(r1 +3* r2);
(*gradient)[1]=2*(2* r1 + r2);
}return f;
}
};// 定义一个通用的边界约束函数
// 它可以表示 x[i] - lower_bound >= 0 这种形式的约束
classBoundConstraint:publicFunctionXf<BoundConstraint>{
public:
int index;// 0 or 1 // 变量的索引 (0 或 1)ScalarType lower_bound;// 边界值BoundConstraint(int i, ScalarType bound):index(i),lower_bound(bound){}// 实现约束函数的求值和梯度计算ScalarType operator()(const VectorType &x, VectorType *grad =nullptr)const{
if(grad){// 如果请求梯度grad->setZero();// 先将梯度向量置零
(*grad)[index]=1.0;// ∇(x[i] - lower_bound)
// 对应 x[i] 的偏导数是 1，其他是 0
}
// 返回约束函数的值
return x[index]- lower_bound;
}
};intmain(){
// 创建一个 L-BFGS-B 求解器实例cppoptlib::solver::Lbfgsb<FunctionExprXf> solver;// optimal solution is suppose to be [1, 1.6] under the box constraints
// 在箱形约束下，最优解应该是 [1, 1.6]
// 将具体的线性回归类包装成通用的 FunctionExpr 类型FunctionExpr f =LinearRegression();// Either use L-BFG-B
// ---------------------------
// 方法一：直接使用 L-BFGS-BEigen::VectorXf x(2);x <<-1,2;// 设置初始猜测值Eigen::VectorXf lb(2);lb <<0,1;// 设置下界 [0, 1]Eigen::VectorXf ub(2);ub <<1,2;// 设置上界 [1, 2]solver.SetBounds(lb, ub);// 将边界信息传递给 L-BFGS-B 求解器constauto initial_state = cppoptlib::function::FunctionState(x);// 创建初始状态
auto[solution, solver_state]= solver.Minimize(f, initial_state);// 运行求解器std::cout <<"argmin "<< solution.x.transpose()<< std::endl;// 打印 L-BFGS-B 找到的解// Or model it as a augmented Lagrangian
// ------------------------------------------
// 方法二：将其建模为增广拉格朗日问题// 将箱形约束 0 <= x[0] <= 1 和 1 <= x[1] <= 2 拆分为四个不等式约束 c(x) >= 0
// 1. x[0] >= 0 => x[0] - 0 >= 0FunctionExpr lb0 =BoundConstraint(0,0.0f);
// 2. x[1] >= 1 => x[1] - 1 >= 0FunctionExpr lb1 =BoundConstraint(1,1.0f);
// 3. x[0] <= 1 => 1 - x[0] >= 0
// 这里通过 -1 * (x[0] - 1) 来实现FunctionExpr ub0 =-1*BoundConstraint(0,1.0f);
// 4. x[1] <= 2 => 2 - x[1] >= 0
// 这里通过 -1 * (x[1] - 2) 来实现FunctionExpr ub1 =-1*BoundConstraint(1,2.0f);// 构造约束优化问题ConstrainedOptimizationProblem prob(f,
/* equality constraints */{},// 没有等式约束
/* inequality constraints */
{lb0,/* x[0] - 0 >= 0 */lb1,/* x[1] - 1 >= 0 */ub0,/* 1 - x[0] >= 0 */ub1,/* 2 - x[1] >= 0 */
});
// 创建一个 L-BFGS 求解器，用作增广拉格朗日法的内层求解器cppoptlib::solver::Lbfgs<FunctionExprXf> unconstrained_solver;
// 创建增广拉格朗日求解器cppoptlib::solver::AugmentedLagrangian aug_solver(prob, unconstrained_solver);
// 初始化增广拉格朗日的状态
// 参数：初始点x，等式约束数量(0)，不等式约束数量(4)，初始惩罚因子(1.0)cppoptlib::solver::AugmentedLagrangeState l_state(x,0,4,1.0f);// Run the agumented solver.
// 运行增广拉格朗日求解器。
auto[aug_solution, aug_solver_state]= aug_solver.Minimize(prob, l_state);
// 打印增广拉格朗日法找到的解std::cout <<"argmin "<< aug_solution.x.transpose()<< std::endl;return0;
}5. simple.cc
这份代码是一个全面的求解器测试平台。它的目的是在一个简单、明确的二次函数上，测试和演示 CppNumericalSolvers 库中提供的几乎所有无约束优化求解器。
通过注释和取消注释不同的 using Solver = ... 行，用户可以轻松地切换不同的求解器（如梯度下降、共轭梯度、牛顿法、BFGS、L-BFGS等），并观察它们的性能、收敛路径和最终结果。代码还演示了如何使用库中的工具来验证用户自定义的梯度和海森矩阵是否正确。
// Copyright 2020, https://github.com/PatWie/CppNumericalSolvers
// 版权所有 2020, https://github.com/PatWie/CppNumericalSolvers#include<iostream>// 引入标准输入输出流库
#include<limits>// 引入数值极限#include"Eigen/Core"// 引入 Eigen 核心库
#include"cppoptlib/function.h"// 引入库中定义函数所需的核心头文件
#include"cppoptlib/solver/bfgs.h"// 引入 BFGS 求解器
#include"cppoptlib/solver/conjugated_gradient_descent.h"// 引入共轭梯度下降求解器
#include"cppoptlib/solver/gradient_descent.h"// 引入梯度下降求解器
#include"cppoptlib/solver/lbfgs.h"// 引入 L-BFGS 求解器
#ifEIGEN_VERSION_AT_LEAST(3,4,0)// 条件编译：只有当 Eigen 版本大于等于 3.4.0 时才包含 L-BFGS-B
#include"cppoptlib/solver/lbfgsb.h"// 引入 L-BFGS-B 求解器
#endif
#include"cppoptlib/solver/nelder_mead.h"// 引入 Nelder-Mead 求解器
#include"cppoptlib/solver/newton_descent.h"// 引入牛顿下降求解器
#include"cppoptlib/utils/derivatives.h"// 引入用于验证导数的工具// 定义一个类型别名，代表一个二阶可微、double类型的函数基类。
template<classF>
using FunctionXd = cppoptlib::function::FunctionCRTP<F,double, cppoptlib::function::DifferentiabilityMode::Second>;
// 定义一个类型别名，代表一个通用的二阶可微、double类型的函数表达式。
using FunctionExprXd2 = cppoptlib::function::FunctionExpr<
double, cppoptlib::function::DifferentiabilityMode::Second>;// 定义一个具体的测试函数：f(x) = 5*x[0]^2 + 100*x[1]^2 + 5
classFunction:publicFunctionXd<Function>{
public:EIGEN_MAKE_ALIGNED_OPERATOR_NEW // Eigen 内存对齐宏
// 以下几行是被注释掉的，但它们明确地定义了函数的属性，
// 实际上这些属性已经通过继承 FunctionXd<Function> 隐式地定义了。
// static constexpr int Dimension = Eigen::Dynamic;
// static constexpr cppoptlib::function::DifferentiabilityMode
// Differentiability = cppoptlib::function::DifferentiabilityMode::Second;
// using ScalarType = double;// 重载 operator()，实现函数求值、梯度和海森矩阵的计算。ScalarType operator()(const VectorType &x, VectorType *gradient =nullptr,MatrixType *hessian =nullptr)const{
if(gradient){// 如果请求梯度gradient->resize(x.size());
(*gradient)[0]=2*5* x[0];// 梯度分量1: ∂f/∂x0 = 10*x0
(*gradient)[1]=2*100* x[1];// 梯度分量2: ∂f/∂x1 = 200*x1
}if(hessian){// 如果请求海森矩阵hessian->resize(x.size(), x.size());
(*hessian)(0,0)=10;
(*hessian)(0,1)=0;
(*hessian)(1,0)=0;
(*hessian)(1,1)=200;
}// 返回函数值
return5* x[0]* x[0]+100* x[1]* x[1]+5;
}
};intmain(){
// 通过取消注释不同的行来选择要使用的求解器。
// using Solver = cppoptlib::solver::GradientDescent<FunctionExprXd2>;
// using Solver = cppoptlib::solver::ConjugatedGradientDescent<FunctionExprXd2>;
// using Solver = cppoptlib::solver::NewtonDescent<FunctionExprXd2>;
// using Solver = cppoptlib::solver::Bfgs<FunctionExprXd2>;
using Solver = cppoptlib::solver::Lbfgs<FunctionExprXd2>;// 当前选择的是 L-BFGS
// using Solver = cppoptlib::solver::Lbfgsb<FunctionExprXd2>;
// using Solver = cppoptlib::solver::NelderMead<FunctionExprXd2>;constexprauto dim =2;// 定义问题维度
// 将具体的 Function 类包装成通用的 FunctionExpr 类型FunctionExprXd2 f =Function();Function::VectorType x(dim);// 创建一个二维向量x <<-10,2;// 设置初始猜测值// ---- 初步测试和验证 ----Function::VectorType gradient = Function::VectorType::Zero(2);
constdouble value =f(x,&gradient);// 计算初始点的函数值和梯度std::cout << value << std::endl;// 打印初始函数值std::cout << gradient << std::endl;// 打印初始梯度// 使用库工具验证解析梯度和海森矩阵是否正确
// IsGradientCorrect 会比较解析梯度和数值梯度std::cout << cppoptlib::utils::IsGradientCorrect(f, x)<< std::endl;// 应该输出 1 (true)
// IsHessianCorrect 会比较解析海森矩阵和数值海森矩阵std::cout << cppoptlib::utils::IsHessianCorrect(f, x)<< std::endl;// 应该输出 1 (true)// ---- 运行优化求解器 ----
// 创建初始状态
constauto initial_state = cppoptlib::function::FunctionState(x);std::cout <<"init "<< initial_state.x.transpose()<< std::endl;// 打印初始点Solver solver;// 创建所选求解器的实例// 设置一个回调函数，用于在每次迭代时打印详细进度solver.SetCallback(cppoptlib::solver::PrintProgressCallback<Solver::FunctionType,Solver::StateType>(std::cout));
// 运行最小化过程
auto[solution, solver_state]= solver.Minimize(f, initial_state);// ---- 打印最终结果 ----std::cout <<"argmin "<< solution.x.transpose()<< std::endl;// 打印找到的最优解 xstd::cout <<"f in argmin "<<f(solution.x)<< std::endl;// 打印最优解处的函数值std::cout <<"iterations "<< solver_state.num_iterations << std::endl;// 打印总迭代次数std::cout <<"status "<< solver_state.status << std::endl;// 打印求解器的最终状态（终止原因）return0;
}6. CMakeLists.txt
这份 CMakeLists.txt 文件是一个用于自动化编译 C++ 项目的脚本。它的主要作用是定义如何编译几个示例程序，并确保它们能找到所需的依赖库（如 Eigen3）。它通过定义一个可复用的 function 来简化多个相似示例的构建过程。
# Function to build examples
# 定义一个名为 build_example 的函数，用于构建示例程序
# 这个函数接受一个参数 `name`，代表示例程序的名称 (例如 "simple")
function(build_example name)
# add_executable 指令用于告诉 CMake 创建一个可执行文件。
# 第一个参数 ${name} 是可执行文件的目标名称 (target name)。
# 第二个参数 ${name}.cc 是用于编译该可执行文件的源文件名。
# 例如，如果 name 是 "simple"，这条命令就是 add_executable(simple simple.cc)。
add_executable(${name}${name}.cc)# target_compile_options 指令用于为指定的目标设置编译选项。
# ${name} 是目标名称。
# PRIVATE 表示这些选项只对这个目标本身生效，不会传递给链接到它的其他目标。
# -std=c++17: 要求编译器使用 C++17 标准。
# -Wall: 开启所有常用的编译器警告 (Warnings all)。
# -Wextra: 开启一些额外的、不被 -Wall 包含的警告。
target_compile_options(${name}PRIVATE -std=c++17 -Wall -Wextra)# target_link_libraries 指令用于为指定的目标链接所需的库。
# ${name} 是目标名称。
# PRIVATE 表示链接关系只对这个目标本身生效。
# CppNumericalSolvers: 链接到本 CppNumericalSolvers 项目本身（假设在主 CMakeLists.txt 中已定义）。
# Eigen3::Eigen: 链接到通过 find_package 找到的 Eigen3 库。
target_link_libraries(${name}PRIVATE CppNumericalSolvers Eigen3::Eigen)
endfunction()# Find Eigen and GoogleTest
# 寻找 Eigen 和 GoogleTest 库
# find_package 指令用于在系统中查找并加载一个外部库的配置。
# Eigen3 是要查找的库的名称。
# REQUIRED 表示如果找不到这个库，CMake 将会报错并停止构建过程。
find_package(Eigen3 REQUIRED)
# 查找 GoogleTest 库（用于单元测试），如果找不到则报错。
# 注意：尽管这里找到了 GTest，但在下面的示例构建中并没有使用它。
# 这可能意味着 GTest 是用于其他部分（如测试目标）的，或者这是一个不完整的脚本片段。
find_package(GTest REQUIRED)# Build examples
# 构建示例程序
# 调用之前定义的 build_example 函数来构建名为 "simple" 的示例。
# 这会执行：
# add_executable(simple simple.cc)
# target_compile_options(simple PRIVATE ...)
# target_link_libraries(simple PRIVATE ...)
build_example(simple)# 调用 build_example 函数来构建名为 "constrained_simple" 的示例。
# 这会执行：
# add_executable(constrained_simple constrained_simple.cc)
# target_compile_options(constrained_simple PRIVATE ...)
# target_link_libraries(constrained_simple PRIVATE ...)
build_example(constrained_simple)这份 CMakeLists.txt 脚本定义了一个自动化的构建流程，其核心逻辑可以概括如下：
1. 依赖声明 (Dependency Declaration)
脚本首先声明了项目的硬性依赖：
- Eigen3
: 一个用于线性代数（矩阵、向量）的 C++ 模板库。这是必需的 (REQUIRED)。
- GTest
: Google 的 C++ 测试框架。这也是必需的 (REQUIRED)，但在此脚本片段中未被示例目标直接使用。
2. 构建蓝图定义 (Build Blueprint Definition)
通过 function(build_example name)，脚本定义了一个名为 build_example 的构建模板或蓝图。这个蓝图规定了任何一个“示例程序”的构建规则：
- 输入 (Input)
:
一个源文件，其名称为
name.cc。
- 输出 (Output)
:
一个可执行文件，其名称为
name。
- 编译规则 (Compilation Rules)
:
必须使用 C++17 标准进行编译。
必须开启
-Wall和-Wextra警告选项，以确保代码质量。
- 链接规则 (Linking Rules)
:
必须链接到
CppNumericalSolvers库本身（提供优化算法）。必须链接到
Eigen3库（提供数学计算支持）。
3. 构建目标实例化 (Build Target Instantiation)
最后，脚本通过调用这个蓝图两次，实例化了两个具体的构建目标：
- 目标
simple:
基于源文件
simple.cc构建。遵循
build_example蓝图的所有编译和链接规则。
- 目标
constrained_simple:
基于源文件
constrained_simple.cc构建。同样遵循
build_example蓝图的所有规则。
总结:
该 CMake 脚本通过一个可复用的函数，高效地定义了两个示例程序（simple 和 constrained_simple）的构建过程。每个示例都会被编译成一个独立的可执行文件，使用 C++17 标准，并链接到 CppNumericalSolvers 和 Eigen3 库。这体现了良好的 CMake 实践，即通过函数来避免重复代码，使构建脚本更简洁、更易于维护。
)
