Covariance Estimation¶
Introduction¶
One way to assess the quality of the solution returned by a nonlinear least squares solver is to analyze the covariance of the solution.
Let us consider the nonlinear regression problem
i.e., the observation \(y\) is a random nonlinear function of the independent variable \(x\) with mean \(f(x)\) and identity covariance. Then the maximum likelihood estimate of \(x\) given observations \(y\) is the solution to the nonlinear least squares problem:
And the covariance of \(x^*\) is given by
Here \(J(x^*)\) is the Jacobian of \(f\) at \(x^*\). The above formula assumes that \(J(x^*)\) has full column rank.
If \(J(x^*)\) is rank deficient, then the covariance matrix \(C(x^*)\) is also rank deficient and is given by the MoorePenrose pseudo inverse.
Note that in the above, we assumed that the covariance matrix for \(y\) was identity. This is an important assumption. If this is not the case and we have
Where \(S\) is a positive semidefinite matrix denoting the covariance of \(y\), then the maximum likelihood problem to be solved is
and the corresponding covariance estimate of \(x^*\) is given by
So, if it is the case that the observations being fitted to have a covariance matrix not equal to identity, then it is the user’s responsibility that the corresponding cost functions are correctly scaled, e.g. in the above case the cost function for this problem should evaluate \(S^{1/2} f(x)\) instead of just \(f(x)\), where \(S^{1/2}\) is the inverse square root of the covariance matrix \(S\).
Gauge Invariance¶
In structure from motion (3D reconstruction) problems, the reconstruction is ambiguous up to a similarity transform. This is known as a Gauge Ambiguity. Handling Gauges correctly requires the use of SVD or custom inversion algorithms. For small problems the user can use the dense algorithm. For more details see the work of Kanatani & Morris [KanataniMorris].
Covariance
¶
Covariance
allows the user to evaluate the covariance for a
nonlinear least squares problem and provides random access to its
blocks. The computation assumes that the cost functions compute
residuals such that their covariance is identity.
Since the computation of the covariance matrix requires computing the
inverse of a potentially large matrix, this can involve a rather large
amount of time and memory. However, it is usually the case that the
user is only interested in a small part of the covariance
matrix. Quite often just the block diagonal. Covariance
allows the user to specify the parts of the covariance matrix that she
is interested in and then uses this information to only compute and
store those parts of the covariance matrix.
Rank of the Jacobian¶
As we noted above, if the Jacobian is rank deficient, then the inverse of \(J'J\) is not defined and instead a pseudo inverse needs to be computed.
The rank deficiency in \(J\) can be structural – columns which are always known to be zero or numerical – depending on the exact values in the Jacobian.
Structural rank deficiency occurs when the problem contains parameter blocks that are constant. This class correctly handles structural rank deficiency like that.
Numerical rank deficiency, where the rank of the matrix cannot be predicted by its sparsity structure and requires looking at its numerical values is more complicated. Here again there are two cases.
The rank deficiency arises from overparameterization. e.g., a four dimensional quaternion used to parameterize \(SO(3)\), which is a three dimensional manifold. In cases like this, the user should use an appropriate
LocalParameterization
. Not only will this lead to better numerical behaviour of the Solver, it will also expose the rank deficiency to theCovariance
object so that it can handle it correctly.More general numerical rank deficiency in the Jacobian requires the computation of the so called Singular Value Decomposition (SVD) of \(J'J\). We do not know how to do this for large sparse matrices efficiently. For small and moderate sized problems this is done using dense linear algebra.

class Covariance::Options¶

int Covariance::Options::num_threads¶
Default:
1
Number of threads to be used for evaluating the Jacobian and estimation of covariance.

SparseLinearAlgebraLibraryType Covariance::Options::sparse_linear_algebra_library_type¶
Default:
SUITE_SPARSE
Ceres Solver is built with support for SuiteSparse andEIGEN_SPARSE
otherwise. Note thatEIGEN_SPARSE
is always available.

CovarianceAlgorithmType Covariance::Options::algorithm_type¶
Default:
SPARSE_QR
Ceres supports two different algorithms for covariance estimation, which represent different tradeoffs in speed, accuracy and reliability.
SPARSE_QR
uses the sparse QR factorization algorithm to compute the decomposition\[\begin{split}QR &= J\\ \left(J^\top J\right)^{1} &= \left(R^\top R\right)^{1}\end{split}\]The speed of this algorithm depends on the sparse linear algebra library being used.
Eigen
’s sparse QR factorization is a moderately fast algorithm suitable for small to medium sized matrices. For best performance we recommend usingSuiteSparseQR
which is enabled by settingCovariance::Options::sparse_linear_algebra_library_type
toSUITE_SPARSE
.SPARSE_QR
cannot compute the covariance if the Jacobian is rank deficient.DENSE_SVD
usesEigen
’sJacobiSVD
to perform the computations. It computes the singular value decomposition\[U D V^\top = J\]and then uses it to compute the pseudo inverse of J’J as
\[(J'J)^{\dagger} = V D^{2\dagger} V^\top\]It is an accurate but slow method and should only be used for small to moderate sized problems. It can handle fullrank as well as rank deficient Jacobians.

double Covariance::Options::column_pivot_threshold¶
Default: \(1\)
During QR factorization, if a column with Euclidean norm less than
column_pivot_threshold
is encountered it is treated as zero.If
column_pivot_threshold < 0
, then an automatic default value of 20*(m+n)*eps*sqrt(max(diag(J’*J))) is used. Here m and n are the number of rows and columns of the Jacobian (J) respectively.This is an advanced option meant for users who know enough about their Jacobian matrices that they can determine a value better than the default.

int Covariance::Options::min_reciprocal_condition_number¶
Default: \(10^{14}\)
If the Jacobian matrix is near singular, then inverting \(J'J\) will result in unreliable results, e.g, if
\[\begin{split}J = \begin{bmatrix} 1.0& 1.0 \\ 1.0& 1.0000001 \end{bmatrix}\end{split}\]which is essentially a rank deficient matrix, we have
\[\begin{split}(J'J)^{1} = \begin{bmatrix} 2.0471e+14& 2.0471e+14 \\ 2.0471e+14& 2.0471e+14 \end{bmatrix}\end{split}\]This is not a useful result. Therefore, by default
Covariance::Compute()
will returnfalse
if a rank deficient Jacobian is encountered. How rank deficiency is detected depends on the algorithm being used.DENSE_SVD
\[\frac{\sigma_{\text{min}}}{\sigma_{\text{max}}} < \sqrt{\text{min_reciprocal_condition_number}}\]where \(\sigma_{\text{min}}\) and \(\sigma_{\text{max}}\) are the minimum and maximum singular values of \(J\) respectively.
SPARSE_QR
\[\operatorname{rank}(J) < \operatorname{num\_col}(J)\]Here \(\operatorname{rank}(J)\) is the estimate of the rank of \(J\) returned by the sparse QR factorization algorithm. It is a fairly reliable indication of rank deficiency.

int Covariance::Options::null_space_rank¶
When using
DENSE_SVD
, the user has more control in dealing with singular and near singular covariance matrices.As mentioned above, when the covariance matrix is near singular, instead of computing the inverse of \(J'J\), the MoorePenrose pseudoinverse of \(J'J\) should be computed.
If \(J'J\) has the eigen decomposition \((\lambda_i, e_i)\), where \(\lambda_i\) is the \(i^\textrm{th}\) eigenvalue and \(e_i\) is the corresponding eigenvector, then the inverse of \(J'J\) is
\[(J'J)^{1} = \sum_i \frac{1}{\lambda_i} e_i e_i'\]and computing the pseudo inverse involves dropping terms from this sum that correspond to small eigenvalues.
How terms are dropped is controlled by min_reciprocal_condition_number and null_space_rank.
If null_space_rank is nonnegative, then the smallest null_space_rank eigenvalue/eigenvectors are dropped irrespective of the magnitude of \(\lambda_i\). If the ratio of the smallest nonzero eigenvalue to the largest eigenvalue in the truncated matrix is still below min_reciprocal_condition_number, then the Covariance::Compute() will fail and return false.
Setting null_space_rank = 1 drops all terms for which
\[\frac{\lambda_i}{\lambda_{\textrm{max}}} < \textrm{min_reciprocal_condition_number}\]This option has no effect on
SPARSE_QR
.

bool Covariance::Options::apply_loss_function¶
Default: true
Even though the residual blocks in the problem may contain loss functions, setting
apply_loss_function
to false will turn off the application of the loss function to the output of the cost function and in turn its effect on the covariance.

class Covariance¶
Covariance::Options
as the name implies is used to control the covariance estimation algorithm. Covariance estimation is a complicated and numerically sensitive procedure. Please read the entire documentation forCovariance::Options
before usingCovariance
.

bool Covariance::Compute(const vector<pair<const double*, const double*>> &covariance_blocks, Problem *problem)¶
Compute a part of the covariance matrix.
The vector
covariance_blocks
, indexes into the covariance matrix blockwise using pairs of parameter blocks. This allows the covariance estimation algorithm to only compute and store these blocks.Since the covariance matrix is symmetric, if the user passes
<block1, block2>
, thenGetCovarianceBlock
can be called withblock1
,block2
as well asblock2
,block1
.covariance_blocks
cannot contain duplicates. Bad things will happen if they do.Note that the list of
covariance_blocks
is only used to determine what parts of the covariance matrix are computed. The full Jacobian is used to do the computation, i.e. they do not have an impact on what part of the Jacobian is used for computation.The return value indicates the success or failure of the covariance computation. Please see the documentation for
Covariance::Options
for more on the conditions under which this function returnsfalse
.

bool GetCovarianceBlock(const double *parameter_block1, const double *parameter_block2, double *covariance_block) const¶
Return the block of the crosscovariance matrix corresponding to
parameter_block1
andparameter_block2
.Compute must be called before the first call to
GetCovarianceBlock
and the pair<parameter_block1, parameter_block2>
OR the pair<parameter_block2, parameter_block1>
must have been present in the vector covariance_blocks whenCompute
was called. OtherwiseGetCovarianceBlock
will return false.covariance_block
must point to a memory location that can store aparameter_block1_size x parameter_block2_size
matrix. The returned covariance will be a rowmajor matrix.

bool GetCovarianceBlockInTangentSpace(const double *parameter_block1, const double *parameter_block2, double *covariance_block) const¶
Return the block of the crosscovariance matrix corresponding to
parameter_block1
andparameter_block2
. Returns crosscovariance in the tangent space if a local parameterization is associated with either parameter block; else returns crosscovariance in the ambient space.Compute must be called before the first call to
GetCovarianceBlock
and the pair<parameter_block1, parameter_block2>
OR the pair<parameter_block2, parameter_block1>
must have been present in the vector covariance_blocks whenCompute
was called. OtherwiseGetCovarianceBlock
will return false.covariance_block
must point to a memory location that can store aparameter_block1_local_size x parameter_block2_local_size
matrix. The returned covariance will be a rowmajor matrix.
Example Usage¶
double x[3];
double y[2];
Problem problem;
problem.AddParameterBlock(x, 3);
problem.AddParameterBlock(y, 2);
<Build Problem>
<Solve Problem>
Covariance::Options options;
Covariance covariance(options);
vector<pair<const double*, const double*> > covariance_blocks;
covariance_blocks.push_back(make_pair(x, x));
covariance_blocks.push_back(make_pair(y, y));
covariance_blocks.push_back(make_pair(x, y));
CHECK(covariance.Compute(covariance_blocks, &problem));
double covariance_xx[3 * 3];
double covariance_yy[2 * 2];
double covariance_xy[3 * 2];
covariance.GetCovarianceBlock(x, x, covariance_xx)
covariance.GetCovarianceBlock(y, y, covariance_yy)
covariance.GetCovarianceBlock(x, y, covariance_xy)