peterzs / igsolver

IGsolver is a C++ library implements several basic optimization algorithms

IGsolver

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.

Unconstrained problem

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.

Constrained problem

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.

Installation

  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.

Usage

Newton-Raphson solver

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

SQP solver

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

Advanced Usage

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.