Sleipnir C++ API
Loading...
Searching...
No Matches
hessian.hpp
1// Copyright (c) Sleipnir contributors
2
3#pragma once
4
5#include <utility>
6
7#include <Eigen/SparseCore>
8#include <gch/small_vector.hpp>
9
10#include "sleipnir/autodiff/adjoint_expression_graph.hpp"
11#include "sleipnir/autodiff/variable.hpp"
12#include "sleipnir/autodiff/variable_matrix.hpp"
13#include "sleipnir/util/assert.hpp"
14#include "sleipnir/util/concepts.hpp"
15#include "sleipnir/util/symbol_exports.hpp"
16
17namespace slp {
18
28template <int UpLo>
29 requires(UpLo == Eigen::Lower) || (UpLo == (Eigen::Lower | Eigen::Upper))
30class SLEIPNIR_DLLEXPORT Hessian {
31 public:
38 Hessian(Variable variable, Variable wrt)
39 : Hessian{std::move(variable), VariableMatrix{std::move(wrt)}} {}
40
49 : m_variables{detail::AdjointExpressionGraph{variable}
50 .generate_gradient_tree(wrt)},
51 m_wrt{wrt} {
52 slp_assert(m_wrt.cols() == 1);
53
54 // Initialize column each expression's adjoint occupies in the Jacobian
55 for (size_t col = 0; col < m_wrt.size(); ++col) {
56 m_wrt[col].expr->col = col;
57 }
58
59 for (auto& variable : m_variables) {
60 m_graphs.emplace_back(variable);
61 }
62
63 // Reset col to -1
64 for (auto& node : m_wrt) {
65 node.expr->col = -1;
66 }
67
68 for (int row = 0; row < m_variables.rows(); ++row) {
69 if (m_variables[row].expr == nullptr) {
70 continue;
71 }
72
73 if (m_variables[row].type() == ExpressionType::LINEAR) {
74 // If the row is linear, compute its gradient once here and cache its
75 // triplets. Constant rows are ignored because their gradients have no
76 // nonzero triplets.
77 m_graphs[row].append_gradient_triplets(m_cached_triplets, row, m_wrt);
78 } else if (m_variables[row].type() > ExpressionType::LINEAR) {
79 // If the row is quadratic or nonlinear, add it to the list of nonlinear
80 // rows to be recomputed in Value().
81 m_nonlinear_rows.emplace_back(row);
82 }
83 }
84
85 if (m_nonlinear_rows.empty()) {
86 m_H.setFromTriplets(m_cached_triplets.begin(), m_cached_triplets.end());
87 if constexpr (UpLo == Eigen::Lower) {
88 m_H = m_H.triangularView<Eigen::Lower>();
89 }
90 }
91 }
92
102 VariableMatrix result{VariableMatrix::empty, m_variables.rows(),
103 m_wrt.rows()};
104
105 for (int row = 0; row < m_variables.rows(); ++row) {
106 auto grad = m_graphs[row].generate_gradient_tree(m_wrt);
107 for (int col = 0; col < m_wrt.rows(); ++col) {
108 if (grad[col].expr != nullptr) {
109 result[row, col] = std::move(grad[col]);
110 } else {
111 result[row, col] = Variable{0.0};
112 }
113 }
114 }
115
116 return result;
117 }
118
124 const Eigen::SparseMatrix<double>& value() {
125 if (m_nonlinear_rows.empty()) {
126 return m_H;
127 }
128
129 for (auto& graph : m_graphs) {
130 graph.update_values();
131 }
132
133 // Copy the cached triplets so triplets added for the nonlinear rows are
134 // thrown away at the end of the function
135 auto triplets = m_cached_triplets;
136
137 // Compute each nonlinear row of the Hessian
138 for (int row : m_nonlinear_rows) {
139 m_graphs[row].append_gradient_triplets(triplets, row, m_wrt);
140 }
141
142 m_H.setFromTriplets(triplets.begin(), triplets.end());
143 if constexpr (UpLo == Eigen::Lower) {
144 m_H = m_H.triangularView<Eigen::Lower>();
145 }
146
147 return m_H;
148 }
149
150 private:
151 VariableMatrix m_variables;
152 VariableMatrix m_wrt;
153
154 gch::small_vector<detail::AdjointExpressionGraph> m_graphs;
155
156 Eigen::SparseMatrix<double> m_H{m_variables.rows(), m_wrt.rows()};
157
158 // Cached triplets for gradients of linear rows
159 gch::small_vector<Eigen::Triplet<double>> m_cached_triplets;
160
161 // List of row indices for nonlinear rows whose graients will be computed in
162 // Value()
163 gch::small_vector<int> m_nonlinear_rows;
164};
165
166} // namespace slp
Definition hessian.hpp:30
Hessian(Variable variable, SleipnirMatrixLike auto wrt)
Definition hessian.hpp:48
VariableMatrix get() const
Definition hessian.hpp:101
Hessian(Variable variable, Variable wrt)
Definition hessian.hpp:38
const Eigen::SparseMatrix< double > & value()
Definition hessian.hpp:124
Definition variable_matrix.hpp:29
int rows() const
Definition variable_matrix.hpp:949
static constexpr empty_t empty
Definition variable_matrix.hpp:39
Definition variable.hpp:40
Definition concepts.hpp:30