This library implements basic optimization algorithms, including Newton's Method for unconstrained problems, SQP method for constrained problems and Gauss-Newton algorithm for non-linear least square problems. It uses Eigen as matrix wrapper and take modern C++ function object as input.
IGsolver::NR_solver
uses Newton's Method to solve unconstrained problems, and implements LM method to resolve singularity. You can customize the shifting matrix of LM method.
IGsolver::GN_solver
uses Gauss-Newton algorithm to solve non-linear least sqaure problems, which also implements LM method.
If the hessian of optimization problem is not positive definite, Cholesky decomposition may fail. To fix this problem, change SOLVER
defined in NR_solver.cpp
to sparse matrix solvers that support non-positive definite matrices, like LU solver or LDLT solver.
IGsolver::SQP_solver
uses Sequential Quadratic Programming (SQP) method to solve constrained problems. It takes augmented Lagrangian as merit function, and LU solver to decomposite SQP matrix.
git clone https://github.com/milkpku/IGsolver
cd IGsolver
mkdir build
cd build
cmake ..
After compilation, you can get library IGsolver.lib
, link it to your own project.
#include <IGsolver/NR_Solver.h>
using namespace IGsolver::NR;
/* f = 0.5 * x1^2 * (x_1^2/6 + 1) + x2 * arctan(x2) - 0.5 * ln(x2^2 + 1)
* f' = [x1^3/3 + x1; arctan(x2)]
* f'' = diag{x1^2 + 1, 1/(1+x2^2)}
*/
Fun_grad_hessian f_grad = [](const dVec& X, double& eval, dVec& grad, SpMat& hessian)
{
double x1 = X(0);
double x2 = X(1);
eval = 0.5 * x1 * x1 * (x1 * x1 / 6 + 1) + x2 * atan(x2) - 0.5 * log(x2 * x2 + 1);
grad.resize(2);
grad(0) = x1 * x1 * x1 / 3 + x1;
grad(1) = atan(x2);
hessian.resize(2, 2);
hessian.insert(0, 0) = x1 * x1 + 1;
hessian.insert(1, 1) = 1 / (1 + x2 * x2);
};
Fun_eval f_eval = [](const dVec& X, const dVec& dX, const double e, double& e_next)
{
double x1 = X(0) + dX(0);
double x2 = X(1) + dX(1);
e_next = 0.5 * x1 * x1 * (x1 * x1 / 6 + 1) + x2 * atan(x2) - 0.5 * log(x2 * x2 + 1);
};
Fun_iter f_iter = [](const int n_iter, const int cut_cnt, const dVec& X, const double& eval, const dVec& dX, const dVec& grad)
{
printf("[iter %d] e = %g, grad = %g, dX = %g, cut_cnt = %d\n",
n_iter, eval, grad.norm(), dX.norm(), cut_cnt);
};
NR_Config config;
dVec sol(2);
sol << 1, 0.7;
NR_solver(sol, f_eval, f_grad, f_iter, config);
#include <IGsolver/SQP_Solver.h>
using namespace IGsolver::SQP;
/* min x1^2 + x2^2
* s.t. x1^2 - x2 -1 = 0
*/
Fun_eval func_eval = [](const dVec& X, double& eval, dVec& c)
{
eval = X.squaredNorm();
c.resize(1);
c(0) = X(0) * X(0) - X(1) - 1;
};
Fun_grad_hess_Jc func_grad_hess_Jc = [](const dVec& X, double& eval,
dVec& grad, SpMat& hessian, dVec& c, SpMat& Jc)
{
eval = X.squaredNorm();
grad = 2 * X;
hessian.resize(2, 2);
hessian.insert(0, 0) = 2;
hessian.insert(1, 1) = 2;
c.resize(1);
c(0) = X(0) * X(0) - X(1) - 1;
Jc.resize(1, 2);
Jc.insert(0, 0) = 2 * X(0);
Jc.insert(0, 1) = -1;
};
Fun_iter func_iter = [](const int n_iter, const int cut_cnt, const dVec& X,
const double& eval, const dVec& dX, const dVec& grad_res, const dVec& c)
{
printf("[iter %d] (x1, x2) = (%g, %g), f(x) = %g, gnorm = %g, c = %g\n",
n_iter, X(0), X(1), eval, grad_res.norm(), c(0));
};
SQP_Config config;
dVec sol = dVec::Random(2);
SQP_solver(sol, func_eval, func_grad_hess_Jc, func_iter, config);
Eigen's sparse matrix solver is slow when matrix becomes big, so please consider other sparse matrix solvers like PARDISO or SuiteSparse. Fortunately, Eigen has convenient API to connect PARDISO and SuiteSparse.
Just redefine SOLVER
in src/{NR|SQP|GN}_solver.cpp
as recommended in comment. For example, adding
#include <Eigen/CholmodSupport>
#define SOLVER Eigen::CholmodSupernodalLLT<SpMat>
in the beginning of src/NR_solver.cpp
or src/GN_solver
replaces Eigen's solver with SuiteSparse Cholesky solver, which speeds up the decomposition for several magnitudes. Also, adding
#include <Eigen/UmfPackSupport>
#define SOLVER Eigen::UmfPackLU<SpMat>
in the beginning of src/SQP_solver.cpp
replaces Eigen's solver with SuiteSparse Umfpack LU solver. For more information of thrid-party sparse matrix solvers, please refer to Eigen's manual of Solving Sparse System.
You can customize config settings in {NR|SQP|GN}_config
, which are self-documented.