Add QR, LU, SVD linear solvers

Author: julien-michot

benchmarks/dense.cpp

@@ -28,7 +28,7 @@ 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.linear_solver = tinyopt::solvers::Options2::LinearSolver::Inverse;
   options.solver.log.print_failure = true;
   BENCHMARK("√2") {
     float x = Vec1::Random()[0];
@@ -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.solver.linear_solver = tinyopt::solvers::Options2::LinearSolver::Inverse;
   static StatCounter<double> counter;
   BENCHMARK("√2") {
     double x = Vec1::Random()[0];  // 0.480009157900 fails to converge

cmake/Options.cmake

@@ -28,7 +28,13 @@ option(TINYOPT_BUILD_PACKAGES "Build packages" OFF)
 option(TINYOPT_BUILD_DOCS "Build documentation" OFF)
 
 
-# Adding Definitions
+# Linear Solvers
+option(TINYOPT_BUILD_SOLVER_LDLT "Build Cholesky/LDLT linear solver" ON)
+option(TINYOPT_BUILD_SOLVER_QR "Build (col-piv) QR linear solver" ON)
+option(TINYOPT_BUILD_SOLVER_LU "Build (partial) LU linear solver" ON)
+option(TINYOPT_BUILD_SOLVER_SVD "Build SVD linear solver" ON)
+
+# Adding Definitions # TODO: use target_compile_definitions()
 if (NOT TINYOPT_ENABLE_FORMATTERS)
   add_definitions(-DTINYOPT_NO_FORMATTERS=1)
 endif ()
@@ -38,3 +44,16 @@ endif ()
 if (TINYOPT_DISABLE_NUMDIFF)
   add_definitions(-DTINYOPT_DISABLE_NUMDIFF=1)
 endif ()
+# Linear solvers
+if (TINYOPT_BUILD_SOLVER_LDLT)
+  add_definitions(-DTINYOPT_BUILD_SOLVER_LDLT=1)
+endif ()
+if (TINYOPT_BUILD_SOLVER_QR)
+  add_definitions(-DTINYOPT_BUILD_SOLVER_QR=1)
+endif ()
+if (TINYOPT_BUILD_SOLVER_LU)
+  add_definitions(-DTINYOPT_BUILD_SOLVER_LU=1)
+endif ()
+if (TINYOPT_BUILD_SOLVER_SVD)
+  add_definitions(-DTINYOPT_BUILD_SOLVER_SVD=1)
+endif ()

include/tinyopt/math.h

@@ -6,6 +6,7 @@
 #include <optional>
 
 #include <tinyopt/types.h>
+#include <Eigen/src/SVD/JacobiSVD.h>
 
 namespace tinyopt {
 
@@ -276,6 +277,110 @@ std::optional<Vector<Scalar, RowsAtCompileTime>> SolveLDLT(
   return X;
 }
 
+/**
+ * @brief Solves the linear system A * X = B for X using LU decomposition.
+ *
+ * @tparam Derived The type of the matrix A.
+ * @tparam Derived2 The type of the vector B.
+ *
+ * @param A The coefficient matrix A.
+ * @param b The right-hand side vector B.
+ *
+ * @return An `std::optional` containing the solution vector X if the system is solvable, or
+ * `std::nullopt` otherwise.
+ */
+template <typename Derived, typename Derived2>
+std::optional<Vector<typename Derived::Scalar, Derived::RowsAtCompileTime>> SolveLU(
+    const MatrixBase<Derived> &A, const MatrixBase<Derived2> &b) {
+  auto lu = A.partialPivLu();
+  if (lu.determinant() != 0) {
+    return lu.solve(b);
+  }
+  return std::nullopt;
+}
+
+/**
+ * @brief Solves the linear system A * X = B for X using QR decomposition.
+ *
+ * @tparam Derived The type of the matrix A.
+ * @tparam Derived2 The type of the vector B.
+ *
+ * @param A The coefficient matrix A.
+ * @param b The right-hand side vector B.
+ *
+ * @return An `std::optional` containing the solution vector X.
+ */
+template <typename Derived, typename Derived2>
+std::optional<Vector<typename Derived::Scalar, Derived::RowsAtCompileTime>> SolveQR(
+    const MatrixBase<Derived> &A, const MatrixBase<Derived2> &b) {
+  return A.colPivHouseholderQr().solve(b);
+}
+
+/**
+ * @brief Solves the linear system A * X = B for X using SVD decomposition.
+ *
+ * @tparam Derived The type of the matrix A.
+ * @tparam Derived2 The type of the vector B.
+ *
+ * @param A The coefficient matrix A.
+ * @param b The right-hand side vector B.
+ *
+ * @return An `std::optional` containing the solution vector X.
+ */
+template <typename Derived, typename Derived2>
+std::optional<Vector<typename Derived::Scalar, Derived::RowsAtCompileTime>> SolveSVD(
+    const MatrixBase<Derived> &A, const MatrixBase<Derived2> &b) {
+  return A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
+}
+
+
+/**
+ * @brief Solves the sparse linear system A * X = B for X using LU decomposition.
+ *
+ * @tparam Scalar The scalar type.
+ * @tparam RowsAtCompileTime The compile-time number of rows of vector B.
+ *
+ * @param A The sparse coefficient matrix A.
+ * @param b The right-hand side vector B.
+ *
+ * @return An `std::optional` containing the solution vector X if the system is solvable, or
+ * `std::nullopt` otherwise.
+ */
+template <typename Scalar, int RowsAtCompileTime = Dynamic>
+std::optional<Vector<Scalar, RowsAtCompileTime>> SolveLU(
+    const SparseMatrix<Scalar> &A, const Vector<Scalar, RowsAtCompileTime> &b) {
+  Eigen::SparseLU<SparseMatrix<Scalar>> solver;
+  solver.compute(A);
+  if (solver.info() != Eigen::Success) return std::nullopt;
+  auto X = solver.solve(b);
+  if (solver.info() != Eigen::Success) return std::nullopt;
+  return X;
+}
+
+/**
+ * @brief Solves the sparse linear system A * X = B for X using QR decomposition.
+ *
+ * @tparam Scalar The scalar type.
+ * @tparam RowsAtCompileTime The compile-time number of rows of vector B.
+ *
+ * @param A The sparse coefficient matrix A.
+ * @param b The right-hand side vector B.
+ *
+ * @return An `std::optional` containing the solution vector X if the system is solvable, or
+ * `std::nullopt` otherwise.
+ */
+template <typename Scalar, int RowsAtCompileTime = Dynamic>
+std::optional<Vector<Scalar, RowsAtCompileTime>> SolveQR(
+    const SparseMatrix<Scalar> &A, const Vector<Scalar, RowsAtCompileTime> &b) {
+  Eigen::SparseQR<SparseMatrix<Scalar>, Eigen::COLAMDOrdering<int>> solver;
+  solver.compute(A);
+  if (solver.info() != Eigen::Success) return std::nullopt;
+  auto X = solver.solve(b);
+  if (solver.info() != Eigen::Success) return std::nullopt;
+  return X;
+}
+
+
 /// Integer square root function for positive integers
 /// Will return `N` for negative or 0 values
 constexpr inline int SQRT(int N) {

include/tinyopt/solvers/gn.h

@@ -40,9 +40,10 @@ class SolverGN
   using Options = nlls::gn::SolverOptions;
 
   explicit SolverGN(const Options &options = {}) : Base(options), options_{options} {
-    // Sparse matrix must use LDLT
+    // Inverse is not supported for sparse matrices
     if constexpr (traits::is_sparse_matrix_v<H_t>) {
-      if (!options.use_ldlt) TINYOPT_LOG("Warning: LDLT must be used with Sparse Matrices");
+      if (options.linear_solver == Options::LinearSolver::Inverse)
+        TINYOPT_LOG("Warning: Inverse is not supported with Sparse Matrices");
     }
   }
 
@@ -149,7 +150,7 @@ class SolverGN
 
     // Fill the lower part if H if needed
     {
-      if (!options_.H_is_full && !options_.use_ldlt) {
+      if (!options_.H_is_full && options_.linear_solver != Options::LinearSolver::LDLT) {
         H_.template triangularView<Lower>() = H_.template triangularView<Upper>().transpose();
       }
     }
@@ -160,17 +161,49 @@ class SolverGN
   inline std::optional<Vector<Scalar, Dims>> Solve() const override {
     if (!this->cost().isValid()) return std::nullopt;
 
+    std::optional<Vector<Scalar, Dims>> dx_;
+
     // Solver linear system
-    if (options_.use_ldlt || traits::is_sparse_matrix_v<H_t>) {
-      const auto dx_ = tinyopt::SolveLDLT(H_, -grad_);
-      if (dx_) return dx_;                                    // Hopefully not a copy...
-    } else if constexpr (!traits::is_sparse_matrix_v<H_t>) {  // Use default inverse
-      if constexpr (Dims == 1) {
-        if (H_(0, 0) > FloatEpsilon<Scalar>()) return -H_.inverse() * grad_;
-        return Vector<Scalar, Dims>::Zero(grad_.size());
-      } else
-        return -H_.inverse() * grad_;
+    switch (options_.linear_solver) {
+#if TINYOPT_BUILD_SOLVER_LDLT
+      case Options::LinearSolver::LDLT:
+        dx_ = tinyopt::SolveLDLT(H_, grad_);
+        break;
+#endif
+#if TINYOPT_BUILD_SOLVER_LU
+      case Options::LinearSolver::LU:
+        dx_ = tinyopt::SolveLU(H_, grad_);
+        break;
+#endif
+#if TINYOPT_BUILD_SOLVER_QR
+      case Options::LinearSolver::QR:
+        dx_ = tinyopt::SolveQR(H_, grad_);
+        break;
+#endif
+#if TINYOPT_BUILD_SOLVER_SVD
+      case Options::LinearSolver::SVD:
+        dx_ = tinyopt::SolveSVD(H_, grad_);
+        break;
+#endif
+      case Options::LinearSolver::Inverse:
+        if constexpr (!traits::is_sparse_matrix_v<H_t>) {  // Use default inverse
+          if constexpr (Dims == 1) {
+            if (H_(0, 0) > FloatEpsilon<Scalar>()) dx_ = H_.inverse() * grad_;
+          } else {
+            dx_ = H_.inverse() * grad_;
+          }
+        }
+        break;
+      default:
+        TINYOPT_LOG("❌ Unsupported linear solver");
+        break;
+    }
+
+    if (dx_) {
+      dx_.value() = (-dx_.value()).eval();
+      return dx_;
     }
+
     // Log on failure
     if (options_.log.enable && options_.log.print_failure) {
       TINYOPT_LOG("❌ Failed solve linear system");

include/tinyopt/solvers/lm.h

@@ -52,9 +52,10 @@ class SolverLM : public tinyopt::solvers::SolverGN<Hessian_t> {
   using Options = nlls::lm::SolverOptions;
 
   explicit SolverLM(const Options &options = {}) : Base(options), options_{options} {
-    // Sparse matrix must use LDLT
+    // Inverse is not supported for sparse matrices
     if constexpr (traits::is_sparse_matrix_v<H_t>) {
-      if (!options.use_ldlt) TINYOPT_LOG("Warning: LDLT must be used with Sparse Matrices");
+      if (options.linear_solver == Options::LinearSolver::Inverse)
+        TINYOPT_LOG("Warning: Inverse is not supported with Sparse Matrices");
     }
     reset();
   }
@@ -104,7 +105,7 @@ class SolverLM : public tinyopt::solvers::SolverGN<Hessian_t> {
 
       // Fill the lower part if H if needed
       if constexpr (!traits::is_sparse_matrix_v<H_t>) {
-        if (!options_.H_is_full && !options_.use_ldlt) {
+        if (!options_.H_is_full && options_.linear_solver != Options::LinearSolver::LDLT) {
           this->H_.template triangularView<Lower>() =
               this->H_.template triangularView<Upper>().transpose();
         }

include/tinyopt/solvers/options.h

@@ -19,7 +19,7 @@ struct Options1 {
    * @name Cost scaling options (mostly for NLLS solvers really)
    * @{
    */
-   struct CostScaling {
+  struct CostScaling {
     bool use_squared_norm = true;  ///< Use squared norm instead of norm (faster)
     bool downscale_by_2 = false;   ///< Rescale the cost by 0.5
     /// Normalize the final error by the number of residuals (after use_squared_norm)
@@ -29,7 +29,7 @@ struct Options1 {
   /** @} */
 
   struct {
-    bool enable = true;  // Enable solver logging
+    bool enable = true;          // Enable solver logging
     bool print_failure = false;  // Log when a failure to solve the linear system happens
   } log;
 };
@@ -46,8 +46,20 @@ struct Options2 : Options1 {
    * @{
    */
 
-  bool use_ldlt = true;   ///< If not, will use H.inverse() without any checks on invertibility
-                          ///< except for Dims==1
+  enum class LinearSolver {
+    LDLT,
+    LU,
+    QR,
+    SVD,
+    Inverse
+  };  ///< Linear solver to use for the linear system
+  LinearSolver linear_solver =
+#ifdef TINYOPT_BUILD_SOLVER_LDLT
+      LinearSolver::LDLT;
+#else
+      LinearSolver::Inverse;
+#endif
+
   bool H_is_full = true;  ///< Specify if H is only Upper triangularly or fully filled
 
   /** @} */

tests/CMakeLists.txt

@@ -23,6 +23,10 @@ add_executable(tinyopt_test_solvers solvers.cpp)
 target_link_libraries(tinyopt_test_solvers PRIVATE ${THIRDPARTY_TEST_LIBS} tinyopt)
 add_test_target(tinyopt_test_solvers)
 
+add_executable(tinyopt_test_linear_solvers linear_solvers.cpp)
+target_link_libraries(tinyopt_test_linear_solvers PRIVATE ${THIRDPARTY_TEST_LIBS} tinyopt)
+add_test_target(tinyopt_test_linear_solvers)
+
 add_executable(tinyopt_test_optimizers optimizers.cpp)
 target_link_libraries(tinyopt_test_optimizers PRIVATE ${THIRDPARTY_TEST_LIBS} tinyopt)
 add_test_target(tinyopt_test_optimizers)