11d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Ceres Solver - A fast non-linear least squares minimizer 21d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Copyright 2013 Google Inc. All rights reserved. 31d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// http://code.google.com/p/ceres-solver/ 41d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 51d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Redistribution and use in source and binary forms, with or without 61d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// modification, are permitted provided that the following conditions are met: 71d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 81d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// * Redistributions of source code must retain the above copyright notice, 91d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// this list of conditions and the following disclaimer. 101d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// * Redistributions in binary form must reproduce the above copyright notice, 111d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// this list of conditions and the following disclaimer in the documentation 121d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// and/or other materials provided with the distribution. 131d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// * Neither the name of Google Inc. nor the names of its contributors may be 141d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// used to endorse or promote products derived from this software without 151d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// specific prior written permission. 161d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 171d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 181d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 191d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 201d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 211d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 221d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 231d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 241d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 251d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 261d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 271d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// POSSIBILITY OF SUCH DAMAGE. 281d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 291d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Author: sameeragarwal@google.com (Sameer Agarwal) 301d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 311d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#ifndef CERES_PUBLIC_COVARIANCE_H_ 321d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#define CERES_PUBLIC_COVARIANCE_H_ 331d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 341d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include <utility> 351d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include <vector> 361d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "ceres/internal/port.h" 371d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "ceres/internal/scoped_ptr.h" 381d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "ceres/types.h" 3979397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez#include "ceres/internal/disable_warnings.h" 401d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 411d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingnamespace ceres { 421d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 431d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingclass Problem; 441d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 451d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingnamespace internal { 461d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingclass CovarianceImpl; 471d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling} // namespace internal 481d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 491d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// WARNING 501d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// ======= 511d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// It is very easy to use this class incorrectly without understanding 521d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// the underlying mathematics. Please read and understand the 531d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// documentation completely before attempting to use this class. 541d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 551d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 561d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// This class allows the user to evaluate the covariance for a 571d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// non-linear least squares problem and provides random access to its 581d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// blocks 591d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 601d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Background 611d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// ========== 621d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// One way to assess the quality of the solution returned by a 631d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// non-linear least squares solve is to analyze the covariance of the 641d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// solution. 651d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 661d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Let us consider the non-linear regression problem 671d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 681d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// y = f(x) + N(0, I) 691d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 701d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// i.e., the observation y is a random non-linear function of the 711d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// independent variable x with mean f(x) and identity covariance. Then 721d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// the maximum likelihood estimate of x given observations y is the 731d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// solution to the non-linear least squares problem: 741d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 751d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// x* = arg min_x |f(x)|^2 761d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 771d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// And the covariance of x* is given by 781d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 791d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// C(x*) = inverse[J'(x*)J(x*)] 801d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 811d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Here J(x*) is the Jacobian of f at x*. The above formula assumes 821d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// that J(x*) has full column rank. 831d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 841d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// If J(x*) is rank deficient, then the covariance matrix C(x*) is 851d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// also rank deficient and is given by 861d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 871d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// C(x*) = pseudoinverse[J'(x*)J(x*)] 881d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 891d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Note that in the above, we assumed that the covariance 901d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// matrix for y was identity. This is an important assumption. If this 911d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// is not the case and we have 921d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 931d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// y = f(x) + N(0, S) 941d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 951d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Where S is a positive semi-definite matrix denoting the covariance 961d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// of y, then the maximum likelihood problem to be solved is 971d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 981d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// x* = arg min_x f'(x) inverse[S] f(x) 991d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1001d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// and the corresponding covariance estimate of x* is given by 1011d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1021d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// C(x*) = inverse[J'(x*) inverse[S] J(x*)] 1031d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1041d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// So, if it is the case that the observations being fitted to have a 1051d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// covariance matrix not equal to identity, then it is the user's 1061d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// responsibility that the corresponding cost functions are correctly 1071d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// scaled, e.g. in the above case the cost function for this problem 1081d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// should evaluate S^{-1/2} f(x) instead of just f(x), where S^{-1/2} 1091d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// is the inverse square root of the covariance matrix S. 1101d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1111d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// This class allows the user to evaluate the covariance for a 1121d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// non-linear least squares problem and provides random access to its 1131d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// blocks. The computation assumes that the CostFunctions compute 1141d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// residuals such that their covariance is identity. 1151d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1161d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Since the computation of the covariance matrix requires computing 1171d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// the inverse of a potentially large matrix, this can involve a 1181d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// rather large amount of time and memory. However, it is usually the 1191d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// case that the user is only interested in a small part of the 1201d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// covariance matrix. Quite often just the block diagonal. This class 1211d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// allows the user to specify the parts of the covariance matrix that 1221d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// she is interested in and then uses this information to only compute 1231d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// and store those parts of the covariance matrix. 1241d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1251d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Rank of the Jacobian 1261d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// -------------------- 1271d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// As we noted above, if the jacobian is rank deficient, then the 1281d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// inverse of J'J is not defined and instead a pseudo inverse needs to 1291d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// be computed. 1301d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1311d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// The rank deficiency in J can be structural -- columns which are 1321d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// always known to be zero or numerical -- depending on the exact 1331d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// values in the Jacobian. 1341d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1351d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Structural rank deficiency occurs when the problem contains 1361d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// parameter blocks that are constant. This class correctly handles 1371d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// structural rank deficiency like that. 1381d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1391d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Numerical rank deficiency, where the rank of the matrix cannot be 1401d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// predicted by its sparsity structure and requires looking at its 1411d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// numerical values is more complicated. Here again there are two 1421d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// cases. 1431d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1441d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// a. The rank deficiency arises from overparameterization. e.g., a 1451d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// four dimensional quaternion used to parameterize SO(3), which is 1461d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// a three dimensional manifold. In cases like this, the user should 1471d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// use an appropriate LocalParameterization. Not only will this lead 1481d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// to better numerical behaviour of the Solver, it will also expose 1491d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// the rank deficiency to the Covariance object so that it can 1501d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// handle it correctly. 1511d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1521d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// b. More general numerical rank deficiency in the Jacobian 1531d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// requires the computation of the so called Singular Value 1541d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Decomposition (SVD) of J'J. We do not know how to do this for 1551d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// large sparse matrices efficiently. For small and moderate sized 1561d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// problems this is done using dense linear algebra. 1571d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1581d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Gauge Invariance 1591d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// ---------------- 1601d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// In structure from motion (3D reconstruction) problems, the 1611d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// reconstruction is ambiguous upto a similarity transform. This is 1621d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// known as a Gauge Ambiguity. Handling Gauges correctly requires the 1631d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// use of SVD or custom inversion algorithms. For small problems the 1641d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// user can use the dense algorithm. For more details see 1651d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1661d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Ken-ichi Kanatani, Daniel D. Morris: Gauges and gauge 1671d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// transformations for uncertainty description of geometric structure 1681d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// with indeterminacy. IEEE Transactions on Information Theory 47(5): 1691d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 2017-2028 (2001) 1701d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1711d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Example Usage 1721d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// ============= 1731d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1741d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// double x[3]; 1751d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// double y[2]; 1761d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1771d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Problem problem; 1781d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// problem.AddParameterBlock(x, 3); 1791d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// problem.AddParameterBlock(y, 2); 1801d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// <Build Problem> 1811d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// <Solve Problem> 1821d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1831d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Covariance::Options options; 1841d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Covariance covariance(options); 1851d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1861d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// vector<pair<const double*, const double*> > covariance_blocks; 1871d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// covariance_blocks.push_back(make_pair(x, x)); 1881d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// covariance_blocks.push_back(make_pair(y, y)); 1891d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// covariance_blocks.push_back(make_pair(x, y)); 1901d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1911d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// CHECK(covariance.Compute(covariance_blocks, &problem)); 1921d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 1931d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// double covariance_xx[3 * 3]; 1941d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// double covariance_yy[2 * 2]; 1951d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// double covariance_xy[3 * 2]; 1961d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// covariance.GetCovarianceBlock(x, x, covariance_xx) 1971d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// covariance.GetCovarianceBlock(y, y, covariance_yy) 1981d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// covariance.GetCovarianceBlock(x, y, covariance_xy) 1991d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 20079397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandezclass CERES_EXPORT Covariance { 2011d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling public: 20279397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez struct CERES_EXPORT Options { 2031d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling Options() 2041d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#ifndef CERES_NO_SUITESPARSE 20579397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez : algorithm_type(SUITE_SPARSE_QR), 2061d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#else 20779397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez : algorithm_type(EIGEN_SPARSE_QR), 2081d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#endif 2091d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling min_reciprocal_condition_number(1e-14), 2101d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling null_space_rank(0), 2111d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling num_threads(1), 2121d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling apply_loss_function(true) { 2131d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling } 2141d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 2151d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Ceres supports three different algorithms for covariance 2161d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // estimation, which represent different tradeoffs in speed, 2171d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // accuracy and reliability. 2181d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2191d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 1. DENSE_SVD uses Eigen's JacobiSVD to perform the 2201d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // computations. It computes the singular value decomposition 2211d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2221d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // U * S * V' = J 2231d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2241d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // and then uses it to compute the pseudo inverse of J'J as 2251d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2261d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // pseudoinverse[J'J]^ = V * pseudoinverse[S] * V' 2271d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2281d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // It is an accurate but slow method and should only be used 2291d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // for small to moderate sized problems. It can handle 2301d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // full-rank as well as rank deficient Jacobians. 2311d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 23279397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // 2. EIGEN_SPARSE_QR uses the sparse QR factorization algorithm 23379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // in Eigen to compute the decomposition 2341d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2351d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Q * R = J 2361d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2371d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // [J'J]^-1 = [R*R']^-1 2381d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 23979397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // It is a moderately fast algorithm for sparse matrices. 2401d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 24179397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // 3. SUITE_SPARSE_QR uses the SuiteSparseQR sparse QR 24279397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // factorization algorithm. It uses dense linear algebra and is 24379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // multi threaded, so for large sparse sparse matrices it is 24479397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // significantly faster than EIGEN_SPARSE_QR. 24579397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // 24679397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // Neither EIGEN_SPARSE_QR not SUITE_SPARSE_QR are capable of 24779397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // computing the covariance if the Jacobian is rank deficient. 2481d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling CovarianceAlgorithmType algorithm_type; 2491d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 2501d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // If the Jacobian matrix is near singular, then inverting J'J 2511d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // will result in unreliable results, e.g, if 2521d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2531d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // J = [1.0 1.0 ] 2541d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // [1.0 1.0000001 ] 2551d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2561d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // which is essentially a rank deficient matrix, we have 2571d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2581d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // inv(J'J) = [ 2.0471e+14 -2.0471e+14] 2591d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // [-2.0471e+14 2.0471e+14] 2601d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2611d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // This is not a useful result. Therefore, by default 2621d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Covariance::Compute will return false if a rank deficient 2631d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Jacobian is encountered. How rank deficiency is detected 2641d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // depends on the algorithm being used. 2651d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2661d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 1. DENSE_SVD 2671d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2681d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // min_sigma / max_sigma < sqrt(min_reciprocal_condition_number) 2691d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2701d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // where min_sigma and max_sigma are the minimum and maxiumum 2711d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // singular values of J respectively. 2721d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 27379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // 2. SUITE_SPARSE_QR and EIGEN_SPARSE_QR 2741d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2751d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // rank(J) < num_col(J) 2761d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2771d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Here rank(J) is the estimate of the rank of J returned by the 27879397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // sparse QR factorization algorithm. It is a fairly reliable 27979397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // indication of rank deficiency. 2801d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2811d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling double min_reciprocal_condition_number; 2821d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 2831d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // When using DENSE_SVD, the user has more control in dealing with 2841d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // singular and near singular covariance matrices. 2851d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2861d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // As mentioned above, when the covariance matrix is near 2871d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // singular, instead of computing the inverse of J'J, the 2881d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Moore-Penrose pseudoinverse of J'J should be computed. 2891d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2901d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // If J'J has the eigen decomposition (lambda_i, e_i), where 2911d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // lambda_i is the i^th eigenvalue and e_i is the corresponding 2921d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // eigenvector, then the inverse of J'J is 2931d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2941d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // inverse[J'J] = sum_i e_i e_i' / lambda_i 2951d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2961d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // and computing the pseudo inverse involves dropping terms from 2971d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // this sum that correspond to small eigenvalues. 2981d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2991d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // How terms are dropped is controlled by 3001d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // min_reciprocal_condition_number and null_space_rank. 3011d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 3021d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // If null_space_rank is non-negative, then the smallest 3031d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // null_space_rank eigenvalue/eigenvectors are dropped 3041d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // irrespective of the magnitude of lambda_i. If the ratio of the 3051d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // smallest non-zero eigenvalue to the largest eigenvalue in the 3061d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // truncated matrix is still below 3071d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // min_reciprocal_condition_number, then the Covariance::Compute() 3081d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // will fail and return false. 3091d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 3101d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Setting null_space_rank = -1 drops all terms for which 3111d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 3121d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // lambda_i / lambda_max < min_reciprocal_condition_number. 3131d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 31479397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // This option has no effect on the SUITE_SPARSE_QR and 31579397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // EIGEN_SPARSE_QR algorithms. 3161d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling int null_space_rank; 3171d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 3181d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling int num_threads; 3191d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 3201d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Even though the residual blocks in the problem may contain loss 3211d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // functions, setting apply_loss_function to false will turn off 3221d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // the application of the loss function to the output of the cost 3231d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // function and in turn its effect on the covariance. 3241d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 3251d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // TODO(sameergaarwal): Expand this based on Jim's experiments. 3261d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling bool apply_loss_function; 3271d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling }; 3281d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 3291d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling explicit Covariance(const Options& options); 3301d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling ~Covariance(); 3311d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 3321d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Compute a part of the covariance matrix. 3331d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 3341d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // The vector covariance_blocks, indexes into the covariance matrix 3351d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // block-wise using pairs of parameter blocks. This allows the 3361d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // covariance estimation algorithm to only compute and store these 3371d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // blocks. 3381d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 3391d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Since the covariance matrix is symmetric, if the user passes 3401d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // (block1, block2), then GetCovarianceBlock can be called with 3411d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // block1, block2 as well as block2, block1. 3421d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 3431d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // covariance_blocks cannot contain duplicates. Bad things will 3441d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // happen if they do. 3451d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 3461d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Note that the list of covariance_blocks is only used to determine 3471d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // what parts of the covariance matrix are computed. The full 3481d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Jacobian is used to do the computation, i.e. they do not have an 3491d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // impact on what part of the Jacobian is used for computation. 3501d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 3511d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // The return value indicates the success or failure of the 3521d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // covariance computation. Please see the documentation for 3531d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Covariance::Options for more on the conditions under which this 3541d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // function returns false. 3551d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling bool Compute( 3561d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling const vector<pair<const double*, const double*> >& covariance_blocks, 3571d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling Problem* problem); 3581d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 3591d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Return the block of the covariance matrix corresponding to 3601d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // parameter_block1 and parameter_block2. 3611d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 3621d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Compute must be called before the first call to 3631d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // GetCovarianceBlock and the pair <parameter_block1, 3641d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // parameter_block2> OR the pair <parameter_block2, 3651d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // parameter_block1> must have been present in the vector 3661d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // covariance_blocks when Compute was called. Otherwise 3671d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // GetCovarianceBlock will return false. 3681d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 3691d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // covariance_block must point to a memory location that can store a 3701d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // parameter_block1_size x parameter_block2_size matrix. The 3711d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // returned covariance will be a row-major matrix. 3721d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling bool GetCovarianceBlock(const double* parameter_block1, 3731d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling const double* parameter_block2, 3741d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling double* covariance_block) const; 3751d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 3761d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling private: 3771d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling internal::scoped_ptr<internal::CovarianceImpl> impl_; 3781d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling}; 3791d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 3801d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling} // namespace ceres 3811d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 38279397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez#include "ceres/internal/reenable_warnings.h" 38379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez 3841d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#endif // CERES_PUBLIC_COVARIANCE_H_ 385