tinyopt 700
Tinyopt, a lightweight, header only optimization library
Loading...
Searching...
No Matches
linear_solvers.h
1// Copyright 2026 Julien Michot.
2// SPDX-License-Identifier: Apache-2.0
3
4#pragma once
5
6#include <algorithm>
7#include <cassert>
8#include <cmath>
9#include <cstdint>
10#include <optional>
11
12#include <Eigen/SparseCholesky>
13
14#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_LLT)
15#include <Eigen/Cholesky>
16#endif
17#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_LU)
18#include <Eigen/LU>
19#include <Eigen/SparseLU>
20#endif
21#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_QR)
22#include <Eigen/QR>
23#include <Eigen/SparseQR>
24#endif
25#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_SVD)
26#include <Eigen/SVD>
27#endif
28
29#include <tinyopt/math/small_linear_solver.h>
30#include <tinyopt/types.h>
31
32#if defined(TINYOPT_ENABLE_SUITESPARSE)
33#include <tinyopt/math/suitesparse.h>
34#endif
35
36namespace tinyopt {
37
38enum class LinearSolverMethod : uint8_t {
39 LDLT,
40 LLT,
41 LU,
42 QR,
43 SVD,
44 SuiteSparse, // Cholmod from SuiteSparse
45};
46
47constexpr bool RequiresFullMatrix(LinearSolverMethod method) {
48 return method != LinearSolverMethod::LDLT && method != LinearSolverMethod::LLT &&
49 method != LinearSolverMethod::SuiteSparse;
50}
51
52template <typename Derived>
53void CompleteSymmetricMatrix(MatrixBase<Derived> &matrix) {
54 matrix.derived() = matrix.derived().template selfadjointView<Eigen::Upper>();
55}
56
57template <typename Scalar, int Options, typename StorageIndex>
58void CompleteSymmetricMatrix(SparseMatrix<Scalar, Options, StorageIndex> &matrix) {
59 SparseMatrix<Scalar, Options, StorageIndex> full =
60 matrix.template selfadjointView<Eigen::Upper>();
61 matrix.swap(full);
62}
63
64template <typename Derived, typename Derived2>
65std::optional<Vector<typename Derived::Scalar, Derived::RowsAtCompileTime>> SolveLDLT(
66 const MatrixBase<Derived> &A, const MatrixBase<Derived2> &b) {
67 using Scalar = typename Derived::Scalar;
68 using Result = Vector<Scalar, Derived::RowsAtCompileTime>;
69 const auto decomposition = A.template selfadjointView<Eigen::Upper>().ldlt();
70 if (decomposition.info() != Eigen::Success || !decomposition.isPositive()) return std::nullopt;
71 Result solution = decomposition.solve(b);
72 if (decomposition.info() != Eigen::Success) return std::nullopt;
73 // NOTE: solution.allFinite() is tested outside of this function in SolveLinearSystem() to avoid
74 // double-checking.
75 return solution;
76}
77
78template <typename Scalar, int RowsAtCompileTime = Dynamic>
79std::optional<Vector<Scalar, RowsAtCompileTime>> SolveLDLT(
80 const SparseMatrix<Scalar> &A, const Vector<Scalar, RowsAtCompileTime> &b) {
81 Eigen::SimplicialLDLT<SparseMatrix<Scalar>, Eigen::Upper> decomposition;
82 decomposition.compute(A);
83 if (decomposition.info() != Eigen::Success) return std::nullopt;
84 Vector<Scalar, RowsAtCompileTime> solution = decomposition.solve(b);
85 if (decomposition.info() != Eigen::Success) return std::nullopt;
86 // NOTE: solution.allFinite() is tested outside of this function in SolveLinearSystem() to avoid
87 // double-checking.
88 return solution;
89}
90
91template <typename Derived, typename Derived2>
92std::optional<Vector<typename Derived::Scalar, Derived::RowsAtCompileTime>> SolveLinearSystem(
93 const MatrixBase<Derived> &A, const MatrixBase<Derived2> &b, LinearSolverMethod method,
94 [[maybe_unused]] double svd_relative_threshold = 0.0) {
95 static_assert(Derived::RowsAtCompileTime == Dynamic || Derived::ColsAtCompileTime == Dynamic ||
96 Derived::RowsAtCompileTime == Derived::ColsAtCompileTime,
97 "SolveLinearSystem requires a square matrix");
98 assert(A.rows() == A.cols());
99 assert(A.rows() == b.rows());
100 if constexpr (Derived::RowsAtCompileTime == Derived::ColsAtCompileTime &&
101 Derived::RowsAtCompileTime >= 1 && Derived::RowsAtCompileTime <= 4) {
102 return SolveSmallDenseSystem(A, b);
103 } else {
104#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_LLT) || defined(TINYOPT_ENABLE_LINEAR_SOLVER_LU) || \
105 defined(TINYOPT_ENABLE_LINEAR_SOLVER_QR) || defined(TINYOPT_ENABLE_LINEAR_SOLVER_SVD)
106 using MatrixType = typename Derived::PlainObject;
107 using Result = Vector<typename Derived::Scalar, Derived::RowsAtCompileTime>;
108#endif
109#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_LU)
110 const auto epsilon = Eigen::NumTraits<typename Derived::Scalar>::epsilon();
111#endif
112
113 switch (method) {
114 case LinearSolverMethod::LDLT:
115#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_LDLT)
116 return SolveLDLT(A, b);
117#else
118 return std::nullopt;
119#endif
120 case LinearSolverMethod::LLT: {
121#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_LLT)
122 Eigen::LLT<MatrixType, Eigen::Upper> decomposition(A);
123 if (decomposition.info() != Eigen::Success) return std::nullopt;
124 Result solution = decomposition.solve(b);
125 if (decomposition.info() != Eigen::Success || !solution.allFinite()) return std::nullopt;
126 return solution;
127#else
128 return std::nullopt;
129#endif
130 }
131 case LinearSolverMethod::SuiteSparse:
132 return std::nullopt;
133#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_LU)
134 case LinearSolverMethod::LU: {
135 Eigen::PartialPivLU<MatrixType> decomposition(A);
136 if (decomposition.rcond() <= epsilon) return std::nullopt;
137 Result solution = decomposition.solve(b);
138 if (!solution.allFinite()) return std::nullopt;
139 return solution;
140 }
141#endif
142#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_QR)
143 case LinearSolverMethod::QR: {
144 Eigen::ColPivHouseholderQR<MatrixType> decomposition(A);
145 if (decomposition.rank() < A.cols()) return std::nullopt;
146 Result solution = decomposition.solve(b);
147 if (!solution.allFinite()) return std::nullopt;
148 return solution;
149 }
150#endif
151#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_SVD)
152 case LinearSolverMethod::SVD: {
153 if (!std::isfinite(svd_relative_threshold) || svd_relative_threshold < 0.0)
154 return std::nullopt;
155 Eigen::JacobiSVD<MatrixType> decomposition(A, Eigen::ComputeFullU | Eigen::ComputeFullV);
156 if (svd_relative_threshold > 0.0)
157 decomposition.setThreshold(static_cast<typename Derived::Scalar>(svd_relative_threshold));
158 Result solution = decomposition.solve(b);
159 if (!solution.allFinite()) return std::nullopt;
160 return solution;
161 }
162#endif
163 default:
164 return std::nullopt;
165 }
166 return std::nullopt;
167 }
168}
169
170template <typename Scalar, int RowsAtCompileTime = Dynamic>
171std::optional<Vector<Scalar, RowsAtCompileTime>> SolveLinearSystem(
172 const SparseMatrix<Scalar> &A, const Vector<Scalar, RowsAtCompileTime> &b,
173 LinearSolverMethod method, double = 0.0) {
174 assert(A.rows() == A.cols());
175 assert(A.rows() == b.rows());
176#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_LLT) || defined(TINYOPT_ENABLE_LINEAR_SOLVER_LU) || \
177 defined(TINYOPT_ENABLE_LINEAR_SOLVER_QR)
178 using Result = Vector<Scalar, RowsAtCompileTime>;
179#endif
180 switch (method) {
181 case LinearSolverMethod::LDLT:
182#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_LDLT)
183 return SolveLDLT(A, b);
184#else
185 return std::nullopt;
186#endif
187 case LinearSolverMethod::SuiteSparse:
188#if defined(TINYOPT_ENABLE_SUITESPARSE)
189 return suitesparse_detail::Solve<Scalar, RowsAtCompileTime>(A, b);
190#else
191 return std::nullopt;
192#endif
193 case LinearSolverMethod::LLT: {
194#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_LLT)
195 Eigen::SimplicialLLT<SparseMatrix<Scalar>, Eigen::Upper> decomposition;
196 decomposition.compute(A);
197 if (decomposition.info() != Eigen::Success) return std::nullopt;
198 Result solution = decomposition.solve(b);
199 if (decomposition.info() != Eigen::Success || !solution.allFinite()) return std::nullopt;
200 return solution;
201#else
202 return std::nullopt;
203#endif
204 }
205#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_LU)
206 case LinearSolverMethod::LU: {
207 Eigen::SparseLU<SparseMatrix<Scalar>> decomposition;
208 decomposition.compute(A);
209 if (decomposition.info() != Eigen::Success) return std::nullopt;
210 Result solution = decomposition.solve(b);
211 if (decomposition.info() != Eigen::Success || !solution.allFinite()) return std::nullopt;
212 return solution;
213 }
214#endif
215#if defined(TINYOPT_ENABLE_LINEAR_SOLVER_QR)
216 case LinearSolverMethod::QR: {
217 Eigen::SparseQR<SparseMatrix<Scalar>, Eigen::COLAMDOrdering<int>> decomposition;
218 decomposition.compute(A);
219 if (decomposition.info() != Eigen::Success || decomposition.rank() < A.cols())
220 return std::nullopt;
221 Result solution = decomposition.solve(b);
222 if (decomposition.info() != Eigen::Success || !solution.allFinite()) return std::nullopt;
223 return solution;
224 }
225#endif
226 case LinearSolverMethod::SVD:
227 return std::nullopt;
228 default:
229 return std::nullopt;
230 }
231 return std::nullopt;
232}
233
234} // namespace tinyopt