refactor: Refactor Options and Optimize function

Author: julien-michot

benchmarks/dense.cpp

@@ -20,16 +20,16 @@
 
 using namespace tinyopt;
 using namespace tinyopt::benchmark;
-using namespace tinyopt::nlls::lm;
+using namespace tinyopt::lm;
 using namespace tinyopt::losses;
 
 static const bool enable_log = false;
 
 TEST_CASE("Float", "[benchmark][fixed][scalar]") {
   auto loss = [](const auto &x) { return x * x - 2.0f; };
   Options options = CreateOptions(enable_log);
-  options.solver.use_ldlt = false;
-  options.solver.log.print_failure = true;
+  options.hessian.use_ldlt = false;
+  options.log.print_failure = true;
   BENCHMARK("√2") {
     float x = Vec1::Random()[0];
     if (enable_log) TINYOPT_LOG("x:{:.12e}", x);
@@ -40,7 +40,7 @@ TEST_CASE("Float", "[benchmark][fixed][scalar]") {
 TEST_CASE("Double", "[benchmark][fixed][scalar]") {
   auto loss = [](const auto &x) { return x * x - 2.0; };
   Options options = CreateOptions(enable_log);
-  options.solver.use_ldlt = false;
+  options.hessian.use_ldlt = false;
   static StatCounter<double> counter;
   BENCHMARK("√2") {
     double x = Vec1::Random()[0];  // 0.480009157900 fails to converge

benchmarks/densef.cpp

@@ -18,7 +18,7 @@
 
 using namespace tinyopt;
 using namespace tinyopt::benchmark;
-using namespace tinyopt::nlls::lm;
+using namespace tinyopt::lm;
 using namespace tinyopt::losses;
 
 TEMPLATE_TEST_CASE("Dense", "[benchmark][fixed][dense][float]", Vec3f, Vec6f, VecXf) {

benchmarks/options.h

@@ -8,7 +8,7 @@
 namespace tinyopt::benchmark {
 
 inline auto CreateOptions(bool enable_log = false) {
-  tinyopt::nlls::lm::Options options;
+  Options options;
   options.max_iters = 10;
 
   options.min_error = 0;  // Ceres does not seem to be using this
@@ -21,8 +21,8 @@ inline auto CreateOptions(bool enable_log = false) {
 
   // No log?
   options.log.enable = enable_log;
-  options.solver.log.enable = enable_log;
-  options.save.H = false;
+  options.log.enable = enable_log;
+  options.hessian.save_last = false;
   return options;
 }
 

benchmarks/sparse.cpp

@@ -14,12 +14,13 @@
 
 #include <tinyopt/tinyopt.h>
 #include "options.h"
+#include "utils.h"
 
 using namespace tinyopt;
 using namespace tinyopt::benchmark;
-using namespace tinyopt::nlls::lm;
+using namespace tinyopt::lm; 
 
-auto simple_loss = [](const auto &x, auto &grad, SparseMatrix<double> &H) {
+auto simple_loss = [](const auto &x, auto &grad, SparseMat &H) {
   const VecX res = 10 * x.array() - 2;
   // Update the gradient and Hessian approx.
   if constexpr (!traits::is_nullptr_v<decltype(grad)>) {
@@ -39,7 +40,7 @@ auto simple_loss = [](const auto &x, auto &grad, SparseMatrix<double> &H) {
     } else if constexpr (0) {  // yet another way, using a dense jacobian
       H = (J.transpose() * J).sparseView();
     } else {  // yet another way, using a sparse jacobian
-      SparseMatrix<double> Js(res.rows(), x.size());
+      SparseMat Js(res.rows(), x.size());
       for (int i = 0; i < x.size(); ++i) Js.coeffRef(i, i) = 10;
       H = Js.transpose() * Js;
     }

docs/API.md

@@ -222,13 +222,13 @@ Optimize(x, loss); // ok, let's optimize
 ### Sparse Systems
 
 Ok so you have quite large and sparse systems? Just tell `Tinyopt` about it by simply
-defining you loss to have a `SparseMatrix<double> &H` type instead of a Dense Matrix or `auto`.
+defining you loss to have a `SparseMat &H` type instead of a Dense Matrix or `auto`.
 
 *NOTE* that automatic differentation is not supported for sparse Hessian matrices but is for first order solvers.
 In that case, simply use this loss signature `auto loss = []<typename T>(const auto &x, SparseMatrix<T> &gradient)`.
 
 ```cpp
-auto loss = [](const auto &x, auto &grad, SparseMatrix<double> &H) {
+auto loss = [](const auto &x, auto &grad, SparseMat &H) {
   // Define your residuals
   const VecX res = 10 * x.array() - 2; // the residuals
   // Update the full gradient matrix, using the Jacobian J of the residuals w.r.t 'x'
@@ -246,7 +246,7 @@ auto loss = [](const auto &x, auto &grad, SparseMatrix<double> &H) {
 
 ```
 
-As an alternative, you can use the `Optimizer<SparseMatrix<MatX>>` class instead of the `Optimize` function.
+As an alternative, you can use the `gn::Optimizer<SparseMatX>` class instead of the `Optimize` function.
 
 There are many ways to fill `H` in Eigen, have a look at `tests/sparse.cpp` for some examples.
 

include/tinyopt/diff/num_diff.h

@@ -189,7 +189,7 @@ auto EstimateNumJac(const X_t &x, const Func &f,
  * double norm = acc_loss(x, g, H);
  *
  * The returned function can be passed to an optimizer, e.g.
- * auto optimizer = Optimizer<SolverGD<Vec2>>();
+ * auto optimizer = Optimizer_<SolverGD<Vec2>>();
  * optimizer(x, acc_loss);
  *
  * @endcode
@@ -275,7 +275,7 @@ auto CreateNumDiffFunc1(X_t &, const ResidualsFunc &residuals,
  * double norm = acc_loss(x, g, H);
  *
  * The returned function can be passed to an optimizer, e.g.
- * auto optimizer = Optimizer<SolverLM<Mat2>>();
+ * auto optimizer = Optimizer_<SolverLM<Mat2>>();
  * optimizer(x, acc_loss);
  *
  * @endcode

include/tinyopt/optimize.h

@@ -4,15 +4,70 @@
 #pragma once
 
 #include <tinyopt/optimizers/optimizer.h>
+#include <tinyopt/optimizers/options.h>
+
+#include <tinyopt/optimizers/optimizers.h>
 
 namespace tinyopt {
 
 /// Simplest interface to optimize `x` and minimize residuals (loss function).
 /// Internally call the optimizer and run the optimization.
-template <typename Optimizer, typename X_t, typename Res_t>
-inline auto Optimize(X_t &x, const Res_t &func, const typename Optimizer::Options &options = {}) {
-  Optimizer optimizer(options);
-  return optimizer(x, func);
+template <typename X_t, typename Res_t>
+inline Output Optimize(X_t &x, const Res_t &func, const Options &options = {}) {
+  // Detect Scalar, supporting at most one nesting level
+  using Scalar = std::conditional_t<
+      std::is_scalar_v<typename traits::params_trait<X_t>::Scalar>,
+      typename traits::params_trait<X_t>::Scalar,
+      typename traits::params_trait<typename traits::params_trait<X_t>::Scalar>::Scalar>;
+  static_assert(std::is_scalar_v<Scalar>);
+  constexpr Index Dims = traits::params_trait<X_t>::Dims;
+
+  // Detect Hessian Type, if it's dense or sparse
+  constexpr bool isDense =
+      std::is_invocable_v<Res_t, const X_t &> ||
+      std::is_invocable_v<Res_t, const X_t &, Vector<Scalar, Dims> &> ||
+      std::is_invocable_v<Res_t, const X_t &, Vector<Scalar, Dims> &, Matrix<Scalar, Dims, Dims> &>;
+
+  using Hessian_t = std::conditional_t<isDense, Matrix<Scalar, Dims, Dims>, SparseMatrix<Scalar>>;
+  using Gradient_t = std::conditional_t<isDense, Vector<Scalar, Dims>, SparseMatrix<Scalar>>;
+
+  constexpr bool secondOrderValid =
+    !std::is_invocable_v<Res_t, const X_t &, Vector<Scalar, Dims> &>;
+
+  // Check if this is an unconstrained first order problem
+  constexpr bool firstOrderAllowed = !secondOrderValid;
+  
+  switch (options.solver_type) {
+    // Second order methods
+    case Options::Solver::GaussNewton:
+    if constexpr (secondOrderValid) {
+      gn::Optimizer<Hessian_t> optimizer(options);
+      return optimizer(x, func);
+    }
+    case Options::Solver::LevenbergMarquardt: 
+    if constexpr (secondOrderValid) {
+      lm::Optimizer<Hessian_t> optimizer(options);
+      return optimizer(x, func);
+    }
+    // First order methods
+    case Options::Solver::GradientDescent: 
+    if constexpr (firstOrderAllowed) {
+      if constexpr (std::is_invocable_v<Res_t, const X_t &>) {
+        const auto res = Res_t(x);
+        if constexpr (traits::is_scalar_v<typename std::decay_t<decltype(res)>>) {
+          gd::Optimizer<Gradient_t> optimizer(options);
+          return optimizer(x, func);
+        } else {
+          throw std::invalid_argument("Error: cost function must return a scalar for Gradient Descent");
+        }
+      } else {
+        gd::Optimizer<Gradient_t> optimizer(options);
+        return optimizer(x, func);
+      }
+    }
+    default: 
+      throw std::invalid_argument("Error: Unknown solver type");
+  }
 }
 
 }  // namespace tinyopt

include/tinyopt/optimizers/gd.h

@@ -4,43 +4,24 @@
 #pragma once
 
 #include <tinyopt/math.h>
-#include <tinyopt/optimize.h>
 
+#include <tinyopt/optimizers/optimizer.h>
 #include <tinyopt/solvers/gd.h>
 
 /// Gradient Descent specific solver, optimizer and their options
 namespace tinyopt::gd {
 
-/// Gradient Descent Optimization Options
-struct Options : Options1 {
-  Options(const Options1 options = {}) : Options1{options} {}
-  gd::SolverOptions solver;
-};
-
 /// Gradient Descent Solver
 template <typename Gradient_t>
 using Solver = solvers::SolverGD<Gradient_t>;
 
-/// Gradient Descent Optimizater type
+/// Gradient Descent Sparse Solver
+template <typename Hessian_t = SparseMat>
+using SparseSolver = solvers::SolverGD<Hessian_t>;
+
+/// Gradient Descent Optimizer
 template <typename Gradient_t>
-using Optimizer = optimizers::Optimizer<Solver<Gradient_t>, Options>;
-
-/// Gradient Descent Optimize function
-template <typename X_t, typename Res_t>
-inline auto Optimize(X_t &x, const Res_t &func, const Options &options = Options()) {
-  using Scalar = std::conditional_t<
-      std::is_scalar_v<typename traits::params_trait<X_t>::Scalar>,
-      typename traits::params_trait<X_t>::Scalar,
-      typename traits::params_trait<typename traits::params_trait<X_t>::Scalar>::Scalar>;
-  static_assert(std::is_scalar_v<Scalar>);
-  constexpr Index Dims = traits::params_trait<X_t>::Dims;
-  // Detect Hessian Type, if it's dense or sparse
-  constexpr bool isDense = std::is_invocable_v<Res_t, const X_t &> ||
-                           std::is_invocable_v<Res_t, const X_t &, Vector<Scalar, Dims> &>;
-  using Gradient_t = std::conditional_t<isDense, Vector<Scalar, Dims>, SparseMatrix<Scalar>>;
-
-  static_assert(Solver<Gradient_t>::FirstOrder);
-  return tinyopt::Optimize<Optimizer<Gradient_t>>(x, func, options);
-}
+using Optimizer = Optimizer_<solvers::SolverGD<Gradient_t>>;
+
 
 }  // namespace tinyopt::gd

include/tinyopt/optimizers/gn.h

@@ -4,47 +4,23 @@
 #pragma once
 
 #include <tinyopt/math.h>
-#include <tinyopt/optimize.h>
 
+#include <tinyopt/optimizers/optimizer.h>
 #include <tinyopt/solvers/gn.h>
 
 /// Gauss-Newton specific solver, optimizer and their options
-namespace tinyopt::nlls::gn {
-
-/// Gauss-Newton Optimization Options
-struct Options : Options2 {
-  Options(const Options2 options = {}) : Options2{options} {}
-  gn::SolverOptions solver;
-};
+namespace tinyopt::gn {
 
 /// Gauss-Newton Solver
 template <typename Hessian_t>
 using Solver = solvers::SolverGN<Hessian_t>;
 
 /// Gauss-Newton Sparse Solver
-template <typename Hessian_t = SparseMatrix<double>>
+template <typename Hessian_t = SparseMat>
 using SparseSolver = solvers::SolverGN<Hessian_t>;
 
-/// Gauss-Newton Optimizater type
+/// Gauss-Newton Optimizer
 template <typename Hessian_t>
-using Optimizer = optimizers::Optimizer<Solver<Hessian_t>, Options>;
-
-/// Gauss-Newton Optimize function
-template <typename X_t, typename Res_t>
-inline auto Optimize(X_t &x, const Res_t &func, const Options &options = Options()) {
-  // Detect Scalar, supporting at most one nesting level
-  using Scalar = std::conditional_t<
-      std::is_scalar_v<typename traits::params_trait<X_t>::Scalar>,
-      typename traits::params_trait<X_t>::Scalar,
-      typename traits::params_trait<typename traits::params_trait<X_t>::Scalar>::Scalar>;
-  static_assert(std::is_scalar_v<Scalar>);
-  constexpr Index Dims = traits::params_trait<X_t>::Dims;
-  // Detect Hessian Type, if it's dense or sparse
-  constexpr bool isDense =
-      std::is_invocable_v<Res_t, const X_t &> ||
-      std::is_invocable_v<Res_t, const X_t &, Vector<Scalar, Dims> &, Matrix<Scalar, Dims, Dims> &>;
-  using Hessian_t = std::conditional_t<isDense, Matrix<Scalar, Dims, Dims>, SparseMatrix<Scalar>>;
-  return tinyopt::Optimize<Optimizer<Hessian_t>>(x, func, options);
-}
-
-}  // namespace tinyopt::nlls::gn
+using Optimizer = Optimizer_<solvers::SolverGN<Hessian_t>>;
+
+}  // namespace tinyopt::gn

include/tinyopt/optimizers/lm.h

@@ -4,49 +4,24 @@
 #pragma once
 
 #include <tinyopt/math.h>
-#include <tinyopt/optimize.h>
 
-#include <tinyopt/optimizers/options.h>
+#include <tinyopt/optimizers/optimizer.h>
 #include <tinyopt/solvers/lm.h>
 #include <type_traits>
 
 /// Levenberg-Marquardt specific solver, optimizer and their options
-namespace tinyopt::nlls::lm {
-
-/// Levenberg-Marquardt Optimization Options
-struct Options : Options2 {
-  Options(const Options2 options = {}) : Options2{options} {}
-  lm::SolverOptions solver;
-};
+namespace tinyopt::lm {
 
 /// Levenberg-Marquardt Solver
 template <typename Hessian_t>
 using Solver = solvers::SolverLM<Hessian_t>;
 
 /// Levenberg-Marquardt Sparse Solver
-template <typename Hessian_t = SparseMatrix<double>>
+template <typename Hessian_t = SparseMat>
 using SparseSolver = solvers::SolverLM<Hessian_t>;
 
-/// Levenberg-Marquardt Optimizater type
+/// Levenberg-Marquardt Optimizer
 template <typename Hessian_t>
-using Optimizer = optimizers::Optimizer<Solver<Hessian_t>, Options>;
-
-/// Levenberg-Marquardt Optimize function
-template <typename X_t, typename Res_t>
-inline auto Optimize(X_t &x, const Res_t &func, const Options &options = Options()) {
-  // Detect Scalar, supporting at most one nesting level
-  using Scalar = std::conditional_t<
-      std::is_scalar_v<typename traits::params_trait<X_t>::Scalar>,
-      typename traits::params_trait<X_t>::Scalar,
-      typename traits::params_trait<typename traits::params_trait<X_t>::Scalar>::Scalar>;
-  static_assert(std::is_scalar_v<Scalar>);
-  constexpr Index Dims = traits::params_trait<X_t>::Dims;
-  // Detect Hessian Type, if it's dense or sparse
-  constexpr bool isDense =
-      std::is_invocable_v<Res_t, const X_t &> ||
-      std::is_invocable_v<Res_t, const X_t &, Vector<Scalar, Dims> &, Matrix<Scalar, Dims, Dims> &>;
-  using Hessian_t = std::conditional_t<isDense, Matrix<Scalar, Dims, Dims>, SparseMatrix<Scalar>>;
-  return tinyopt::Optimize<Optimizer<Hessian_t>>(x, func, options);
-}
-
-}  // namespace tinyopt::nlls::lm
+using Optimizer = Optimizer_<solvers::SolverLM<Hessian_t>>;
+
+}  // namespace tinyopt::lm