example 1

#include "pressio/solvers_nonlinear.hpp"
#include <Eigen/Core>
#include <Eigen/Sparse>
#include <optional>

struct MyProblem
{
  using state_type    = Eigen::VectorXd;
  using residual_type = state_type;
  using jacobian_type = Eigen::SparseMatrix<double>;

  state_type createState() const {
    state_type a(2);
    a.setZero();
    return a;
  }

  residual_type createResidual() const {
    residual_type a(2);
    a.setZero();
    return a;
  }

  jacobian_type createJacobian() const {
    jacobian_type a(2, 2);
    a.setZero();
    return a;
  }

  void residualAndJacobian(const state_type& x,
			   residual_type& res,
			   std::optional<jacobian_type*> Jin) const
  {
    res(0) = x(0)*x(0) + x(1)*x(1) - 4.0;
    res(1) = x(0)*x(1) - 1.0;

    if (Jin){
     auto & jac = *Jin.value();
      jac.coeffRef(0, 0) = 2*x(0);
      jac.coeffRef(0, 1) = 2*x(1);
      // Have incorrect entries so that line search is required
      jac.coeffRef(1, 0) = 0.1*x(1);
      jac.coeffRef(1, 1) = 0.1*x(0);
    }
  }
};

class MyLinearSolver{
  using A_t = typename MyProblem::jacobian_type;
  using b_t = typename MyProblem::residual_type;
  using x_t = typename MyProblem::state_type;

public:
  void solve(A_t const & A, b_t const & b, x_t & x)
  {
    Eigen::BiCGSTAB< A_t > solver;
    solver.compute(A);
    x = solver.solve(b);
  }
};

int main()
{
  PRESSIOLOG_INITIALIZE(pressiolog::LogLevel::sparse);

  namespace pnonls = pressio::nlsol;

  MyProblem problem;
  MyLinearSolver linSolver;
  auto solver = pnonls::create_newton_solver(problem, linSolver);
  const auto updateMethod = pnonls::Update::BacktrackStrictlyDecreasingObjective;
  solver.setUpdateCriterion(updateMethod);

  auto x = problem.createState();
  x(0) = 0.3; x(1) = 0.4;
  solver.solve(x);

  // check solution
  std::cout << "Computed solution: \n "
	    << "["
	    << std::setprecision(14)
	    << x(0) << " "
	    << x(1) << " " << "] \n"
	    << "Expected solution: \n "
	    << "[1.93185165, 0.51763809] \n";

  std::cout << "PASSED\n";
  PRESSIOLOG_FINALIZE();
  return 0;
}

example 2

Solve one variant of the Rosenbrock problem via Gauss-Newton and Levenberg-Marquardt.

#include "pressio/solvers_linear.hpp"
#include "pressio/solvers_nonlinear.hpp"
#include <Eigen/Core>
#include <optional>

struct MyRosenbrockSystem
{
  using scalar_type   = double;
  using state_type    = Eigen::VectorXd;
  using residual_type = state_type;
  using jacobian_type = Eigen::MatrixXd;

  state_type createState() const{
    auto s = state_type(4);
    s.setZero();
    return s;
  };

  residual_type createResidual() const{
    auto r = residual_type(6);
    r.setZero();
    return r;
  }

  jacobian_type createJacobian() const{
    auto j = jacobian_type(6, 4);
    j.setZero();
    return j;
  }

  void residualAndJacobian(const state_type& x,
			   residual_type& res,
			   std::optional<jacobian_type*> Jin) const
  {
    auto x1 = x(0);
    auto x2 = x(1);
    auto x3 = x(2);
    auto x4 = x(3);

    res(0) = 10.*(x4 - x3*x3);
    res(1) = 10.*(x3 - x2*x2);
    res(2) = 10.*(x2 - x1*x1);
    res(3) = (1.-x1);
    res(4) = (1.-x2);
    res(5) = (1.-x3);

    if (Jin){
      auto & J = *Jin.value();
      J.setZero();
      J(0,2) = -20.*x3; J(0,3) = 10.; // first row
      J(1,1) = -20.*x2; J(1,2) = 10.; // second row
      J(2,0) = -20.*x1; J(2,1) = 10.; // etc
      J(3,0) = -1.;
      J(4,1) = -1.;
      J(5,2) = -1.;
    }
  }
};

int main()
{
  PRESSIOLOG_INITIALIZE(pressiolog::LogLevel::sparse);

  namespace pls    = pressio::linsol;
  namespace pnonls = pressio::nlsol;

  MyRosenbrockSystem problem;

  using hessian_t = Eigen::MatrixXd;
  using linear_solver_tag = pls::direct::HouseholderQR;
  using linear_solver_t   = pls::Solver<linear_solver_tag, hessian_t>;
  linear_solver_t linearSolver;

  auto x0 = problem.createState();
  x0(0) = 0.0; x0(1) = 1.1; x0(2) = 1.2; x0(3) = 1.5;

  {
    auto solver = pnonls::create_gauss_newton_solver(problem, linearSolver);
    solver.setStopTolerance(1e-5);

    auto x = pressio::ops::clone(x0);
    solver.solve(x);

    // check solution
    std::cout << "Computed solution: \n "
	      << "["
	      << std::setprecision(14)
	      << x(0) << " "
	      << x(1) << " "
	      << x(2) << " "
	      << x(3) << " " << "] \n"
	      << "exact minimum at: \n "
	      << "[1.0, 1.0, 1.0, 1.0] \n";
  }

  {
    auto solver = pnonls::create_levenberg_marquardt_solver(problem, linearSolver);
    solver.setStopTolerance(1e-5);
    solver.setUpdateCriterion(pnonls::Update::LMSchedule2);

    auto x = pressio::ops::clone(x0);
    solver.solve(x);

    // check solution
    std::cout << "Computed solution: \n "
	      << "["
	      << std::setprecision(14)
	      << x(0) << " "
	      << x(1) << " "
	      << x(2) << " "
	      << x(3) << " " << "] \n"
	      << "exact minimum at: \n "
	      << "[1.0, 1.0, 1.0, 1.0] \n";
  }

  std::cout << "PASSED\n";
  PRESSIOLOG_FINALIZE();
  return 0;
}