@@ -6,140 +6,166 @@ The library is designed for efficiency, modern C++ compliance, and easy
66integration into existing projects. It is distributed under a permissive
77license, making it suitable for commercial use.
88
9- ### Example Usage: Minimizing the Rosenbrock Function
9+ Key Features:
10+ - ** Unconstrained Optimization** : Supports algorithms such as ** Gradient Descent, Conjugate Gradient, Newton's Method, BFGS, L-BFGS, and Nelder-Mead** .
11+ - ** Constrained Optimization** : Provides support for ** L-BFGS-B** and ** Augmented Lagrangian methods** .
12+ - ** Expression Templates** : Enable efficient evaluation of function expressions without unnecessary allocations.
13+ - ** First- and Second-Order Methods** : Support for both ** gradient-based** and ** Hessian-based** optimizations.
14+ - ** Differentiation Utilities** : Tools to check gradient and Hessian correctness.
15+ - ** Header-Only** : Easily integrate into any C++17 project with no dependencies.
1016
11- Let minimize the classic Rosenbrock function using the BFGS solver.
17+
18+ ## Quick Start
19+
20+ ### ** Minimizing a Function (BFGS Solver)**
1221
1322``` cpp
14- /* *
15- * @brief Alias for a 2D function supporting second-order differentiability in cppoptlib.
16- *
17- * This defines a function template specialized for a 2-dimensional input vector,
18- * which supports evaluation of the function value and its derivatives.
19- *
20- * The differentiability level is determined by the Differentiability enum:
21- * - Differentiability::First: Supports first-order derivatives (i.e., the gradient)
22- * without computing the Hessian. This is useful for optimization methods that do not
23- * require second-order information, saving computational effort.
24- * - Differentiability::Second: Supports second-order derivatives, meaning that both the
25- * gradient and Hessian are computed. This level is needed for methods that rely on curvature
26- * information.
27- *
28- * In this alias, Differentiability::Second is used, so both the gradient and Hessian are
29- * assumed to be implemented.
30- */
31- using Function2D = cppoptlib::function::Function<double , 2 , cppoptlib::function::Differentiability::Second>;
32-
33- /* *
34- * @brief Implementation of a quadratic function with optional gradient and Hessian computation.
35- *
36- * The function is defined as:
37- *
38- * f(x) = 5 * x₀² + 100 * x₁² + 5
39- *
40- * Its gradient is given by:
41- *
42- * ∇f(x) = [10 * x₀, 200 * x₁]ᵀ
43- *
44- * And its Hessian is:
45- *
46- * ∇²f(x) = [ [10, 0],
47- * [ 0, 200] ]
48- *
49- * This implementation computes the gradient and Hessian only if the corresponding pointers
50- * are provided by the caller.
51- */
52- class Function : public Function2D {
23+ #include < iostream>
24+ #include " cppoptlib/function.h"
25+ #include " cppoptlib/solver/bfgs.h"
26+
27+ using namespace cppoptlib ::function;
28+
29+ class ExampleFunction : public FunctionXd<ExampleFunction > {
5330public:
5431 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
5532
56- /**
57- * @brief Evaluates the function and optionally its gradient and Hessian.
58- *
59- * @param x Input vector.
60- * @param gradient (Optional) Pointer to store the computed gradient.
61- * @param hessian (Optional) Pointer to store the computed Hessian.
62- * @return The function value f(x).
63- * /
64- scalar_t operator()(const vector_t &x, vector_t * gradient = nullptr,
65- matrix_t * hessian = nullptr) const override {
66-
67- // Even for functions declared as Differentiability::First, the gradient is not always required.
68- // To save computation, we only calculate and store the gradient if a non-null pointer is provided.
33+ /**
34+ * The function is defined as:
35+ *
36+ * f(x) = 5 * x₀² + 100 * x₁² + 5
37+ *
38+ * Its gradient is given by:
39+ *
40+ * ∇f(x) = [10 * x₀, 200 * x₁]ᵀ
41+ *
42+ * And its Hessian is:
43+ *
44+ * ∇²f(x) = [ [10, 0],
45+ * [ 0, 200] ]
46+ */
47+ ScalarType operator()(const VectorType &x, VectorType * gradient = nullptr,
48+ MatrixType * hessian = nullptr) const override {
6949 if (gradient) {
70- gradient->resize(x.size());
71- // The gradient components:
72- // ∂f/∂x₀ = 2 * 5 * x₀ = 10 * x₀
73- // ∂f/∂x₁ = 2 * 100 * x₁ = 200 * x₁
50+ * gradient = VectorType::Zero(2);
7451 (* gradient)[ 0] = 2 * 5 * x[ 0] ;
7552 (* gradient)[ 1] = 2 * 100 * x[ 1] ;
7653 }
7754
78- // If the Hessian is requested, compute and store it.
79- // (This is applicable only for Differentiability::Second)
8055 if (hessian) {
81- hessian->resize(x.size(), x.size());
82- // Set the Hessian components:
83- // ∂²f/∂x₀² = 10, ∂²f/∂x₁² = 200, and the off-diagonals are 0.
56+ *hessian = MatrixType::Zero(2, 2);
8457 (*hessian)(0, 0) = 10;
8558 (*hessian)(0, 1) = 0;
8659 (*hessian)(1, 0) = 0;
8760 (*hessian)(1, 1) = 200;
8861 }
8962
90- // Return the function value: f(x) = 5* x₀² + 100* x₁² + 5.
9163 return 5 * x[0 ] * x[0 ] + 100 * x[1 ] * x[1 ] + 5 ;
9264 }
9365};
9466
95-
96- int main(int argc, char const * argv[ ] ) {
97- // Create an instance of the Rosenbrock function.
98- Rosenbrock f;
99-
100- // Initial guess for the solution.
67+ int main () {
68+ ExampleFunction f;
10169 Eigen::VectorXd x(2);
10270 x << -1, 2;
103- std::cout << "Initial point: " << x.transpose() << std::endl;
104-
105- // Evaluate
106- Eigen::VectorXd gradient(2);
107- double value = f(x, &gradient);
108- std::cout << "Function value at initial point: " << value << std::endl;
109- std::cout << "Gradient at initial point: " << gradient << std::endl;
11071
111- // Minimize the Rosenbrock function using the BFGS solver.
112- using Solver = cppoptlib::solver::Bfgs<Rosenbrock>;
72+ using Solver = cppoptlib::solver::Bfgs<ExampleFunction>;
11373 auto init_state = f.GetState(x);
11474 Solver solver;
11575 auto [solution_state, solver_progress] = solver.Minimize(f, init_state);
11676
117- // Display the results of the optimization.
11877 std::cout << "Optimal solution: " << solution_state.x.transpose() << std::endl;
119- std::cout << "Optimal function value: " << f(solution_state.x) << std::endl;
78+ return 0;
79+ }
80+ ```
12081
121- std::cout << "Number of iterations: " << solver_progress.num_iterations << std::endl;
122- std::cout << "Solver status: " << solver_progress.status << std::endl;
82+ ### ** Solving a Constrained Optimization Problem**
12383
124- return 0;
84+ CppNumericalSolvers allows solving constrained optimization problems using the
85+ ** Augmented Lagrangian** method. Below, we solve the problem:
86+
87+ ```
88+ min f(x) = x_0 + x_1
89+ ```
90+
91+ subject to the constraint:
92+
93+ ```
94+ x_0^2 + x_1^2 = 2
95+ ```
96+
97+ ``` cpp
98+ #include < iostream>
99+ #include " cppoptlib/function.h"
100+ #include " cppoptlib/solver/augmented_lagrangian.h"
101+ #include " cppoptlib/solver/lbfgs.h"
102+
103+ using namespace cppoptlib ::function;
104+ using namespace cppoptlib::solver;
105+
106+ class SumObjective : public FunctionXd<SumObjective > {
107+ public:
108+ ScalarType operator()(const VectorType &x, VectorType * gradient = nullptr) const {
109+ if (gradient) * gradient = VectorType::Ones(2);
110+ return x.sum();
111+ }
112+ };
113+
114+ class Circle : public FunctionXd<Circle > {
115+ public:
116+ ScalarType operator()(const VectorType &x, VectorType * gradient = nullptr) const {
117+ if (gradient) * gradient = 2 * x;
118+ return x.squaredNorm();
119+ }
120+ };
121+
122+ int main() {
123+ SumObjective::VectorType x(2);
124+ x << 2, 10;
125+
126+ FunctionExpr objective = SumObjective();
127+ FunctionExpr circle = Circle();
128+
129+ ConstrainedOptimizationProblem prob(
130+ objective,
131+ // 1. Equality constraint: circle(x) - 2 == 0, forcing the solution onto
132+ // the circle's boundary.
133+ {FunctionExpr(circle - 2)},
134+ // 2. Inequality constraint: 2 - circle(x) >= 0, ensuring the solution
135+ // remains inside the circle.
136+ {FunctionExpr(2 - circle)});
137+
138+ Lbfgs<FunctionExpr2d1 > inner_solver;
139+ AugmentedLagrangian solver(prob, inner_solver);
140+
141+ AugmentedLagrangeState<double > l_state(x, /* num_eq=* /1,
142+ /* num_ineq=* /1,
143+ /* penalty=* /1.0);
144+ auto [ solution, solver_state] = solver.Minimize(prob, l_state);
145+
146+ std::cout << "Optimal x: " << solution.x.transpose() << std::endl;
147+ return 0;
125148}
149+
126150```
127151
128- You can easily adapt this code for your specific optimization problem by
129- defining your objective function and selecting an appropriate solver from
130- CppNumericalSolvers. Dive into the implementation for more details on available
131- solvers and advanced usage.
132152
153+ ## Use in Your Project
154+
155+ Either use Bazel or CMake:
133156
134- See the examples for a constrained problem.
157+ ```sh
158+ git clone https://github.com/PatWie/CppNumericalSolvers.git
159+ ```
160+
161+ Compile with:
162+
163+ ``` sh
164+ g++ -std=c++17 -Ipath/to/CppNumericalSolvers/include myfile.cpp -o myprogram
165+ ```
135166
136- ### References
137167
138- **L-BFGS-B**: A LIMITED MEMORY ALGORITHM FOR BOUND CONSTRAINED OPTIMIZATION
139- _Richard H. Byrd, Peihuang Lu, Jorge Nocedal and Ciyou Zhu_
140168
141- **L-BFGS**: Numerical Optimization, 2nd ed. New York: Springer
142- _J. Nocedal and S. J. Wright_
143169
144170### Citing this implementation
145171
0 commit comments