gtsam 4.1.1
gtsam
RegularImplicitSchurFactor.h
Go to the documentation of this file.
1
8#pragma once
9
13
14#include <iosfwd>
15#include <map>
16#include <string>
17#include <vector>
18
19namespace gtsam {
20
38template<class CAMERA>
40
41public:
43 typedef boost::shared_ptr<This> shared_ptr;
44
45protected:
46
47 // This factor is closely related to a CameraSet
48 typedef CameraSet<CAMERA> Set;
49
50 typedef typename CAMERA::Measurement Z;
51 static const int D = traits<CAMERA>::dimension;
52 static const int ZDim = traits<Z>::dimension;
53
54 typedef Eigen::Matrix<double, ZDim, D> MatrixZD;
55 typedef Eigen::Matrix<double, D, D> MatrixDD;
56 typedef std::vector<MatrixZD, Eigen::aligned_allocator<MatrixZD> > FBlocks;
57
58 FBlocks FBlocks_;
59 const Matrix PointCovariance_;
60 const Matrix E_;
61 const Vector b_;
62
63public:
64
67 }
68
70
80 RegularImplicitSchurFactor(const KeyVector& keys, const FBlocks& Fs,
81 const Matrix& E, const Matrix& P, const Vector& b)
82 : GaussianFactor(keys), FBlocks_(Fs), PointCovariance_(P), E_(E), b_(b) {}
83
86 }
87
88 const FBlocks& Fs() const {
89 return FBlocks_;
90 }
91
92 const Matrix& E() const {
93 return E_;
94 }
95
96 const Vector& b() const {
97 return b_;
98 }
99
100 const Matrix& getPointCovariance() const {
101 return PointCovariance_;
102 }
103
105 void print(const std::string& s = "", const KeyFormatter& keyFormatter =
106 DefaultKeyFormatter) const override {
107 std::cout << " RegularImplicitSchurFactor " << std::endl;
108 Factor::print(s);
109 for (size_t pos = 0; pos < size(); ++pos) {
110 std::cout << "Fblock:\n" << FBlocks_[pos] << std::endl;
111 }
112 std::cout << "PointCovariance:\n" << PointCovariance_ << std::endl;
113 std::cout << "E:\n" << E_ << std::endl;
114 std::cout << "b:\n" << b_.transpose() << std::endl;
115 }
116
118 bool equals(const GaussianFactor& lf, double tol) const override {
119 const This* f = dynamic_cast<const This*>(&lf);
120 if (!f)
121 return false;
122 for (size_t k = 0; k < FBlocks_.size(); ++k) {
123 if (keys_[k] != f->keys_[k])
124 return false;
125 if (!equal_with_abs_tol(FBlocks_[k], f->FBlocks_[k], tol))
126 return false;
127 }
128 return equal_with_abs_tol(PointCovariance_, f->PointCovariance_, tol)
129 && equal_with_abs_tol(E_, f->E_, tol)
130 && equal_with_abs_tol(b_, f->b_, tol);
131 }
132
134 DenseIndex getDim(const_iterator variable) const override {
135 return D;
136 }
137
139 SymmetricBlockMatrix* info) const override {
140 throw std::runtime_error(
141 "RegularImplicitSchurFactor::updateHessian non implemented");
142 }
143 Matrix augmentedJacobian() const override {
144 throw std::runtime_error(
145 "RegularImplicitSchurFactor::augmentedJacobian non implemented");
146 return Matrix();
147 }
148 std::pair<Matrix, Vector> jacobian() const override {
149 throw std::runtime_error(
150 "RegularImplicitSchurFactor::jacobian non implemented");
151 return std::make_pair(Matrix(), Vector());
152 }
153
155 Matrix augmentedInformation() const override {
156 // Do the Schur complement
157 SymmetricBlockMatrix augmentedHessian =
159 return augmentedHessian.selfadjointView();
160 }
161
163 Matrix information() const override {
164 Matrix augmented = augmentedInformation();
165 int m = this->keys_.size();
166 size_t M = D * m;
167 return augmented.block(0, 0, M, M);
168 }
169
172
174 void hessianDiagonalAdd(VectorValues &d) const override {
175 // diag(Hessian) = diag(F' * (I - E * PointCov * E') * F);
176 for (size_t k = 0; k < size(); ++k) { // for each camera
177 Key j = keys_[k];
178
179 // Calculate Fj'*Ej for the current camera (observing a single point)
180 // D x 3 = (D x ZDim) * (ZDim x 3)
181 const MatrixZD& Fj = FBlocks_[k];
182 Eigen::Matrix<double, D, 3> FtE = Fj.transpose()
183 * E_.block<ZDim, 3>(ZDim * k, 0);
184
185 Eigen::Matrix<double, D, 1> dj;
186 for (int k = 0; k < D; ++k) { // for each diagonal element of the camera hessian
187 // Vector column_k_Fj = Fj.col(k);
188 dj(k) = Fj.col(k).squaredNorm(); // dot(column_k_Fj, column_k_Fj);
189 // Vector column_k_FtE = FtE.row(k);
190 // (1 x 1) = (1 x 3) * (3 * 3) * (3 x 1)
191 dj(k) -= FtE.row(k) * PointCovariance_ * FtE.row(k).transpose();
192 }
193
194 auto result = d.emplace(j, dj);
195 if(!result.second) {
196 result.first->second += dj;
197 }
198 }
199 }
200
205 void hessianDiagonal(double* d) const override {
206 // diag(Hessian) = diag(F' * (I - E * PointCov * E') * F);
207 // Use eigen magic to access raw memory
208 typedef Eigen::Matrix<double, D, 1> DVector;
209 typedef Eigen::Map<DVector> DMap;
210
211 for (size_t pos = 0; pos < size(); ++pos) { // for each camera in the factor
212 Key j = keys_[pos];
213
214 // Calculate Fj'*Ej for the current camera (observing a single point)
215 // D x 3 = (D x ZDim) * (ZDim x 3)
216 const MatrixZD& Fj = FBlocks_[pos];
217 Eigen::Matrix<double, D, 3> FtE = Fj.transpose()
218 * E_.block<ZDim, 3>(ZDim * pos, 0);
219
220 DVector dj;
221 for (int k = 0; k < D; ++k) { // for each diagonal element of the camera hessian
222 dj(k) = Fj.col(k).squaredNorm();
223 // (1 x 1) = (1 x 3) * (3 * 3) * (3 x 1)
224 dj(k) -= FtE.row(k) * PointCovariance_ * FtE.row(k).transpose();
225 }
226 DMap(d + D * j) += dj;
227 }
228 }
229
231 std::map<Key, Matrix> hessianBlockDiagonal() const override {
232 std::map<Key, Matrix> blocks;
233 // F'*(I - E*P*E')*F
234 for (size_t pos = 0; pos < size(); ++pos) {
235 Key j = keys_[pos];
236 // F'*F - F'*E*P*E'*F e.g. (9*2)*(2*9) - (9*2)*(2*3)*(3*3)*(3*2)*(2*9)
237 const MatrixZD& Fj = FBlocks_[pos];
238 // Eigen::Matrix<double, D, 3> FtE = Fj.transpose()
239 // * E_.block<ZDim, 3>(ZDim * pos, 0);
240 // blocks[j] = Fj.transpose() * Fj
241 // - FtE * PointCovariance_ * FtE.transpose();
242
243 const Matrix23& Ej = E_.block<ZDim, 3>(ZDim * pos, 0);
244 blocks[j] = Fj.transpose()
245 * (Fj - Ej * PointCovariance_ * Ej.transpose() * Fj);
246
247 // F'*(I - E*P*E')*F, TODO: this should work, but it does not :-(
248 // static const Eigen::Matrix<double, ZDim, ZDim> I2 = eye(ZDim);
249 // Matrix2 Q = //
250 // I2 - E_.block<ZDim, 3>(ZDim * pos, 0) * PointCovariance_ * E_.block<ZDim, 3>(ZDim * pos, 0).transpose();
251 // blocks[j] = Fj.transpose() * Q * Fj;
252 }
253 return blocks;
254 }
255
257 return boost::make_shared<RegularImplicitSchurFactor<CAMERA> >(keys_,
259 throw std::runtime_error(
260 "RegularImplicitSchurFactor::clone non implemented");
261 }
262
263 bool empty() const override {
264 return false;
265 }
266
268 return boost::make_shared<RegularImplicitSchurFactor<CAMERA> >(keys_,
270 throw std::runtime_error(
271 "RegularImplicitSchurFactor::negate non implemented");
272 }
273
274 // Raw Vector version of y += F'*alpha*(I - E*P*E')*F*x, for testing
275 static
276 void multiplyHessianAdd(const Matrix& F, const Matrix& E,
277 const Matrix& PointCovariance, double alpha, const Vector& x, Vector& y) {
278 Vector e1 = F * x;
279 Vector d1 = E.transpose() * e1;
280 Vector d2 = PointCovariance * d1;
281 Vector e2 = E * d2;
282 Vector e3 = alpha * (e1 - e2);
283 y += F.transpose() * e3;
284 }
285
286 typedef std::vector<Vector2, Eigen::aligned_allocator<Vector2>> Error2s;
287
291 void projectError2(const Error2s& e1, Error2s& e2) const {
292
293 // d1 = E.transpose() * (e1-ZDim*b) = (3*2m)*2m
294 Vector3 d1;
295 d1.setZero();
296 for (size_t k = 0; k < size(); k++)
297 d1 += E_.block<ZDim, 3>(ZDim * k, 0).transpose()
298 * (e1[k] - ZDim * b_.segment<ZDim>(k * ZDim));
299
300 // d2 = E.transpose() * e1 = (3*2m)*2m
301 Vector3 d2 = PointCovariance_ * d1;
302
303 // e3 = alpha*(e1 - E*d2) = 1*[2m-(2m*3)*3]
304 for (size_t k = 0; k < size(); k++)
305 e2[k] = e1[k] - ZDim * b_.segment<ZDim>(k * ZDim)
306 - E_.block<ZDim, 3>(ZDim * k, 0) * d2;
307 }
308
309 /*
310 * This definition matches the linearized error in the Hessian Factor:
311 * LinError(x) = x'*H*x - 2*x'*eta + f
312 * with:
313 * H = F' * (I-E'*P*E) * F = F' * Q * F
314 * eta = F' * (I-E'*P*E) * b = F' * Q * b
315 * f = nonlinear error
316 * (x'*H*x - 2*x'*eta + f) = x'*F'*Q*F*x - 2*x'*F'*Q *b + f = x'*F'*Q*(F*x - 2*b) + f
317 */
318 double error(const VectorValues& x) const override {
319
320 // resize does not do malloc if correct size
321 e1.resize(size());
322 e2.resize(size());
323
324 // e1 = F * x - b = (2m*dm)*dm
325 for (size_t k = 0; k < size(); ++k)
326 e1[k] = FBlocks_[k] * x.at(keys_[k]);
327 projectError2(e1, e2);
328
329 double result = 0;
330 for (size_t k = 0; k < size(); ++k)
331 result += dot(e1[k], e2[k]);
332
333 double f = b_.squaredNorm();
334 return 0.5 * (result + f);
335 }
336
337 // needed to be GaussianFactor - (I - E*P*E')*(F*x - b)
338 // This is wrong and does not match the definition in Hessian,
339 // but it matches the definition of the Jacobian factor (JF)
340 double errorJF(const VectorValues& x) const {
341
342 // resize does not do malloc if correct size
343 e1.resize(size());
344 e2.resize(size());
345
346 // e1 = F * x - b = (2m*dm)*dm
347 for (size_t k = 0; k < size(); ++k)
348 e1[k] = FBlocks_[k] * x.at(keys_[k]) - b_.segment<ZDim>(k * ZDim);
349 projectError(e1, e2);
350
351 double result = 0;
352 for (size_t k = 0; k < size(); ++k)
353 result += dot(e2[k], e2[k]);
354
355 // std::cout << "implicitFactor::error result " << result << std::endl;
356 return 0.5 * result;
357 }
361 void projectError(const Error2s& e1, Error2s& e2) const {
362
363 // d1 = E.transpose() * e1 = (3*2m)*2m
364 Vector3 d1;
365 d1.setZero();
366 for (size_t k = 0; k < size(); k++)
367 d1 += E_.block<ZDim, 3>(ZDim * k, 0).transpose() * e1[k];
368
369 // d2 = E.transpose() * e1 = (3*2m)*2m
370 Vector3 d2 = PointCovariance_ * d1;
371
372 // e3 = alpha*(e1 - E*d2) = 1*[2m-(2m*3)*3]
373 for (size_t k = 0; k < size(); k++)
374 e2[k] = e1[k] - E_.block<ZDim, 3>(ZDim * k, 0) * d2;
375 }
376
378 mutable Error2s e1, e2;
379
384 void multiplyHessianAdd(double alpha, const double* x, double* y) const {
385
386 // Use eigen magic to access raw memory
387 typedef Eigen::Matrix<double, D, 1> DVector;
388 typedef Eigen::Map<DVector> DMap;
389 typedef Eigen::Map<const DVector> ConstDMap;
390
391 // resize does not do malloc if correct size
392 e1.resize(size());
393 e2.resize(size());
394
395 // e1 = F * x = (2m*dm)*dm
396 for (size_t k = 0; k < size(); ++k) {
397 Key key = keys_[k];
398 e1[k] = FBlocks_[k] * ConstDMap(x + D * key);
399 }
400
401 projectError(e1, e2);
402
403 // y += F.transpose()*e2 = (2d*2m)*2m
404 for (size_t k = 0; k < size(); ++k) {
405 Key key = keys_[k];
406 DMap(y + D * key) += FBlocks_[k].transpose() * alpha * e2[k];
407 }
408 }
409
410 void multiplyHessianAdd(double alpha, const double* x, double* y,
411 std::vector<size_t> keys) const {
412 }
413
417 void multiplyHessianAdd(double alpha, const VectorValues& x,
418 VectorValues& y) const override {
419
420 // resize does not do malloc if correct size
421 e1.resize(size());
422 e2.resize(size());
423
424 // e1 = F * x = (2m*dm)*dm
425 for (size_t k = 0; k < size(); ++k)
426 e1[k] = FBlocks_[k] * x.at(keys_[k]);
427
428 projectError(e1, e2);
429
430 // y += F.transpose()*e2 = (2d*2m)*2m
431 for (size_t k = 0; k < size(); ++k) {
432 Key key = keys_[k];
433 static const Vector empty;
434 std::pair<VectorValues::iterator, bool> it = y.tryInsert(key, empty);
435 Vector& yi = it.first->second;
436 // Create the value as a zero vector if it does not exist.
437 if (it.second)
438 yi = Vector::Zero(FBlocks_[k].cols());
439 yi += FBlocks_[k].transpose() * alpha * e2[k];
440 }
441 }
442
446 void multiplyHessianDummy(double alpha, const VectorValues& x,
447 VectorValues& y) const {
448
449 for (size_t k = 0; k < size(); ++k) {
450 static const Vector empty;
451 Key key = keys_[k];
452 std::pair<VectorValues::iterator, bool> it = y.tryInsert(key, empty);
453 Vector& yi = it.first->second;
454 yi = x.at(key);
455 }
456 }
457
461 VectorValues gradientAtZero() const override {
462 // calculate Q*b
463 e1.resize(size());
464 e2.resize(size());
465 for (size_t k = 0; k < size(); k++)
466 e1[k] = b_.segment<ZDim>(ZDim * k);
467 projectError(e1, e2);
468
469 // g = F.transpose()*e2
470 VectorValues g;
471 for (size_t k = 0; k < size(); ++k) {
472 Key key = keys_[k];
473 g.insert(key, -FBlocks_[k].transpose() * e2[k]);
474 }
475
476 // return it
477 return g;
478 }
479
483 void gradientAtZero(double* d) const override {
484
485 // Use eigen magic to access raw memory
486 typedef Eigen::Matrix<double, D, 1> DVector;
487 typedef Eigen::Map<DVector> DMap;
488
489 // calculate Q*b
490 e1.resize(size());
491 e2.resize(size());
492 for (size_t k = 0; k < size(); k++)
493 e1[k] = b_.segment<ZDim>(ZDim * k);
494 projectError(e1, e2);
495
496 for (size_t k = 0; k < size(); ++k) { // for each camera in the factor
497 Key j = keys_[k];
498 DMap(d + D * j) += -FBlocks_[k].transpose() * e2[k];
499 }
500 }
501
503 Vector gradient(Key key, const VectorValues& x) const override {
504 throw std::runtime_error(
505 "gradient for RegularImplicitSchurFactor is not implemented yet");
506 }
507
508};
509// end class RegularImplicitSchurFactor
510
511template<class CAMERA>
513
514template<class CAMERA>
516
517// traits
518template<class CAMERA> struct traits<RegularImplicitSchurFactor<CAMERA> > : public Testable<
519 RegularImplicitSchurFactor<CAMERA> > {
520};
521
522}
523
Base class to create smart factors on poses or cameras.
Factor Graph Values.
Global functions in a separate testing namespace.
Definition: chartTesting.h:28
FastVector< Key > KeyVector
Define collection type once and for all - also used in wrappers.
Definition: Key.h:86
ptrdiff_t DenseIndex
The index type for Eigen objects.
Definition: types.h:75
double dot(const V1 &a, const V2 &b)
Dot product.
Definition: Vector.h:194
std::uint64_t Key
Integer nonlinear key type.
Definition: types.h:69
std::function< std::string(Key)> KeyFormatter
Typedef for a function to format a key, i.e. to convert it to a string.
Definition: Key.h:35
bool equal_with_abs_tol(const Eigen::DenseBase< MATRIX > &A, const Eigen::DenseBase< MATRIX > &B, double tol=1e-9)
equals with a tolerance
Definition: Matrix.h:84
A manifold defines a space in which there is a notion of a linear tangent space that can be centered ...
Definition: concepts.h:30
Definition: SymmetricBlockMatrix.h:52
Eigen::SelfAdjointView< constBlock, Eigen::Upper > selfadjointView(DenseIndex I, DenseIndex J) const
Return the square sub-matrix that contains blocks(i:j, i:j).
Definition: SymmetricBlockMatrix.h:161
A helper that implements the traits interface for GTSAM types.
Definition: Testable.h:151
A set of cameras, all with their own calibration.
Definition: CameraSet.h:35
static SymmetricBlockMatrix SchurComplement(const std::vector< Eigen::Matrix< double, ZDim, ND >, Eigen::aligned_allocator< Eigen::Matrix< double, ZDim, ND > > > &Fs, const Matrix &E, const Eigen::Matrix< double, N, N > &P, const Vector &b)
Do Schur complement, given Jacobian as Fs,E,P, return SymmetricBlockMatrix G = F' * F - F' * E * P * ...
Definition: CameraSet.h:153
This is the base class for all factor types.
Definition: Factor.h:56
const KeyVector & keys() const
Access the factor's involved variable keys.
Definition: Factor.h:125
KeyVector keys_
The keys involved in this factor.
Definition: Factor.h:73
virtual void print(const std::string &s="Factor", const KeyFormatter &formatter=DefaultKeyFormatter) const
print
Definition: Factor.cpp:29
KeyVector::const_iterator const_iterator
Const iterator over keys.
Definition: Factor.h:68
size_t size() const
Definition: Factor.h:136
An abstract virtual base class for JacobianFactor and HessianFactor.
Definition: GaussianFactor.h:39
boost::shared_ptr< This > shared_ptr
shared_ptr to this class
Definition: GaussianFactor.h:42
VectorValues hessianDiagonal() const
Return the diagonal of the Hessian for this factor.
Definition: GaussianFactor.cpp:27
This class represents a collection of vector-valued variables associated each with a unique integer i...
Definition: VectorValues.h:74
std::pair< VectorValues::iterator, bool > emplace(Key j, Args &&... args)
Emplace a vector value with key j.
Definition: VectorValues.h:183
iterator insert(const std::pair< Key, Vector > &key_value)
Insert a vector value with key j.
Definition: VectorValues.cpp:91
Vector & at(Key j)
Read/write access to the vector value with key j, throws std::out_of_range if j does not exist,...
Definition: VectorValues.h:137
std::pair< iterator, bool > tryInsert(Key j, const Vector &value)
insert that mimics the STL map insert - if the value already exists, the map is not modified and an i...
Definition: VectorValues.h:207
RegularImplicitSchurFactor.
Definition: RegularImplicitSchurFactor.h:39
const Matrix E_
The 2m*3 E Jacobian with respect to the point.
Definition: RegularImplicitSchurFactor.h:60
void projectError(const Error2s &e1, Error2s &e2) const
Calculate corrected error Q*e = (I - E*P*E')*e.
Definition: RegularImplicitSchurFactor.h:361
bool empty() const override
Test whether the factor is empty.
Definition: RegularImplicitSchurFactor.h:263
GaussianFactor::shared_ptr clone() const override
Clone a factor (make a deep copy)
Definition: RegularImplicitSchurFactor.h:256
void hessianDiagonalAdd(VectorValues &d) const override
Add the diagonal of the Hessian for this factor to existing VectorValues.
Definition: RegularImplicitSchurFactor.h:174
RegularImplicitSchurFactor(const KeyVector &keys, const FBlocks &Fs, const Matrix &E, const Matrix &P, const Vector &b)
Construct from blocks of F, E, inv(E'*E), and RHS vector b.
Definition: RegularImplicitSchurFactor.h:80
RegularImplicitSchurFactor()
Constructor.
Definition: RegularImplicitSchurFactor.h:66
const Vector b_
2m-dimensional RHS vector
Definition: RegularImplicitSchurFactor.h:61
std::pair< Matrix, Vector > jacobian() const override
Return the dense Jacobian and right-hand-side , with the noise models baked into A and b.
Definition: RegularImplicitSchurFactor.h:148
Eigen::Matrix< double, D, D > MatrixDD
camera Hessian
Definition: RegularImplicitSchurFactor.h:55
void gradientAtZero(double *d) const override
Calculate gradient, which is -F'Q*b, see paper - RAW MEMORY ACCESS.
Definition: RegularImplicitSchurFactor.h:483
static const int ZDim
Measurement dimension.
Definition: RegularImplicitSchurFactor.h:52
Matrix information() const override
Compute full information matrix
Definition: RegularImplicitSchurFactor.h:163
~RegularImplicitSchurFactor() override
Destructor.
Definition: RegularImplicitSchurFactor.h:85
void updateHessian(const KeyVector &keys, SymmetricBlockMatrix *info) const override
Update an information matrix by adding the information corresponding to this factor (used internally ...
Definition: RegularImplicitSchurFactor.h:138
FBlocks FBlocks_
All ZDim*D F blocks (one for each camera)
Definition: RegularImplicitSchurFactor.h:58
VectorValues gradientAtZero() const override
Calculate gradient, which is -F'Q*b, see paper.
Definition: RegularImplicitSchurFactor.h:461
Eigen::Matrix< double, ZDim, D > MatrixZD
type of an F block
Definition: RegularImplicitSchurFactor.h:54
const Matrix PointCovariance_
the 3*3 matrix P = inv(E'E) (2*2 if degenerate)
Definition: RegularImplicitSchurFactor.h:59
void multiplyHessianAdd(double alpha, const VectorValues &x, VectorValues &y) const override
Hessian-vector multiply, i.e.
Definition: RegularImplicitSchurFactor.h:417
std::map< Key, Matrix > hessianBlockDiagonal() const override
Return the block diagonal of the Hessian for this factor.
Definition: RegularImplicitSchurFactor.h:231
Matrix augmentedInformation() const override
Compute full augmented information matrix
Definition: RegularImplicitSchurFactor.h:155
double error(const VectorValues &x) const override
Print for testable.
Definition: RegularImplicitSchurFactor.h:318
GaussianFactor::shared_ptr negate() const override
Construct the corresponding anti-factor to negate information stored stored in this factor.
Definition: RegularImplicitSchurFactor.h:267
boost::shared_ptr< This > shared_ptr
shared_ptr to this class
Definition: RegularImplicitSchurFactor.h:43
Error2s e1
Scratch space for multiplyHessianAdd.
Definition: RegularImplicitSchurFactor.h:378
void print(const std::string &s="", const KeyFormatter &keyFormatter=DefaultKeyFormatter) const override
print
Definition: RegularImplicitSchurFactor.h:105
DenseIndex getDim(const_iterator variable) const override
Degrees of freedom of camera.
Definition: RegularImplicitSchurFactor.h:134
bool equals(const GaussianFactor &lf, double tol) const override
equals
Definition: RegularImplicitSchurFactor.h:118
Matrix augmentedJacobian() const override
Return a dense Jacobian matrix, augmented with b with the noise models baked into A and b.
Definition: RegularImplicitSchurFactor.h:143
RegularImplicitSchurFactor This
Typedef to this class.
Definition: RegularImplicitSchurFactor.h:42
static const int D
Camera dimension.
Definition: RegularImplicitSchurFactor.h:51
void projectError2(const Error2s &e1, Error2s &e2) const
Calculate corrected error Q*(e-ZDim*b) = (I - E*P*E')*(e-ZDim*b)
Definition: RegularImplicitSchurFactor.h:291
void multiplyHessianAdd(double alpha, const double *x, double *y) const
double* Hessian-vector multiply, i.e.
Definition: RegularImplicitSchurFactor.h:384
void hessianDiagonal(double *d) const override
add the contribution of this factor to the diagonal of the hessian d(output) = d(input) + deltaHessia...
Definition: RegularImplicitSchurFactor.h:205
void multiplyHessianDummy(double alpha, const VectorValues &x, VectorValues &y) const
Dummy version to measure overhead of key access.
Definition: RegularImplicitSchurFactor.h:446
Vector gradient(Key key, const VectorValues &x) const override
Gradient wrt a key at any values.
Definition: RegularImplicitSchurFactor.h:503