协方差估计
介绍
鉴别最小二乘问题返回值好坏的一个方法是分析解的协方差。我们考虑非线性回归问题。
换句话上,观察值 y是一个对变量x的具有单位协方差的随机非线性函数,然后,对于观测y,x的最大似然估计是非线性最小二乘问题:
如果
注意上面我们假设y的协方差矩阵是单位矩阵。这是一个很重要的假设。如果不符合这个条件,那么
相应的x的协方差估计将如:
所以,如果观察的结果有一个不等于单位矩阵的协方差矩阵,用户需要自己缩放成本函数。例如上面的问题应该求
规范不变性
在SfM问题(3D重建)中,重建无法区分相似变换,这被称作 Gauge Ambiguity。正确地处理Gauges需要使用SVD或自定义的反演算法。对于小问题用户可以使用稠密算法。更多细节请看Kanatani & Morris的论文Gauges and gauge transformations for uncertainty description of geometric structure with indeterminacy。
协方差
Covariance 允许用户求解非线性最小二乘问题的协方差,并提供对其块的随机访问。计算假设成本函数计算残差的协方差是单位矩阵。
由于计算协方差矩阵一个潜在的大矩阵的逆,这可能需要大量的时间和内存。然而,通常用户只对协方差的很小一部分感兴趣,通常是对角块。Covariance 允许用户指定协方差矩阵中感兴趣的部分,并仅计算和存储这一部分。
Jacobian 的秩
正如我们上面所提到的,如果Jacobian矩阵的秩不足,那么
当问题包含常量参数块时,就会出现结构秩的不足。这个类正确地处理结构秩的不足。
数值秩的不足,由于它的稀疏结构不能预测矩阵的秩,并且需要看它的数值,这是比较复杂的。这里有两种情况。
- 秩不足由于过参数化,一个四维的四元数来参数化三维的SO(3),像这样的情况,用户需要适当的本地参数化。不仅仅因为这将导致更好解的数字表现,它还会将秩缺失暴露给协方差对象,使其能够正确处理。
- 在Jacobian矩阵中,更一般的数值秩缺乏需要计算
J′J 的奇异值分解(SVD)的计算。我们不知道如何有效地利用稀疏矩阵。对于小的和中等大小的问题,这是通过稠密的线性代数来完成的。
Covariance::Options
class Covariance::Options
Covariance::Options::num_threads
Default: 1
计算Jacobian和协方差所使用的线程数。
Covariance::Options::sparse_linear_algebra_library_type
Default: SUITE_SPARSE
Ceres Solver支持 SuiteSparse 和 EIGEN_SPARSE,EIGEN_SPARSE总是可用的。
Covariance::Options::algorithm_type
Default: SPARSE_QR
Ceres Solver支持两种不同的算法计算协方差,代表速度,精度和可靠性的权衡。
- SPARSE_QR 用稀疏QR分解来计算分解
QR(J⊤J)−1=J=(R⊤R)−1
算法的速度依赖于使用的稀疏线性代数库。Eigen的稀疏QR分解是一种中等速度的算法,适合小到中等大小的矩阵。为了更好的性能我们推荐使用 SuiteSparseQR,通过将Covaraince::Options::sparse_linear_algebra_library_type设置为 SUITE_SPARSE。
SPARSE_QR不能计算秩不足的情况。 - DENSE_SVD 使用 Eigen的JacobiSVD去计算。他计算SVD
USV⊤=J
然后用它去求解J’J 的伪逆矩阵
(J′J)†=VS†V⊤
这个方法很精确也很慢,只应该用来计算中等大小的问题。计算秩不足和满秩情况都可以用。
Covariance::Options::min_reciprocal_condition_number
Default:
如果Jacobian矩阵接近奇异,那么
这本质上是一个秩不足的矩阵,我们有
这不是一个有用的结果。因此,默认 Covariance::Compute() 将会返回 false,如果遇到一个秩不足的矩阵.鉴别秩不足依赖于算法
- DENSE_SVD
σminσmax<min_reciprocal_condition_number‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√
σmin 和σmax 分别表示J 的最小和最大奇异值。 - SPARSE_QR
rank(J)<num_col(J)
rank(J) 是稀疏QR分解计算的J的秩。这是一个相当可靠的秩不足的迹象。
Covariance::Options::null_space_rank
当使用 DENSE_SVD,用户在处理正定和近似正定协方差矩阵有更多的控制。我们计算
如果
计算伪逆涉及到从这个和中去掉对应的特征值小的项。
怎样去掉项,由min_reciprocal_condition_number和null_space_rank控制。
如果null_space_rank是非负的,最小的 null_space_rank 特征值/特征向量被丢弃在相应的
这个选项对SPARSE_QR没有影响。
Covariance::Options::apply_loss_function
Default: true
即使问题的残差块包含损失函数,将apply_loss_function设置成false将关闭应用损失函数。
Covariance::Options::apply_loss_function
Default: true
即使问题的残差块包含损失函数,将apply_loss_function设置成false将关闭应用损失函数。
class Covariance
Covariance::Options 正如名字一样,是用来控制协方差估计算法的。Covariance求解是一种复杂的数值敏感的计算。请参阅Covariance::Options 的所有文档,在使用Covariance之前。
Covariance::Compute
计算一部分协方差矩阵。向量 covariance_blocks,利用参数块对协方差矩阵块进行索引。这使得协方差估计算法只能计算和存储这些块。
由于协方差矩阵是对称的,如果用户传递block1, block2,那么GetCovarianceBlock可以通过 block1, block2和block2, block1调用。
covariance_blocks包含重复,如果重复会发生不好的事情。
注意,covariance_blocks的列表只用于确定计算协方差矩阵的哪些部分。完整的Jacobian矩阵用来做计算,也就是说,它们对Jacobian矩阵的计算部分没有影响。
返回值表示协方差计算的成功或失败。请参阅 Covariance::Options的文档关于更多该函数返回false的条件的选项。
GetCovarianceBlock
返回与parameter_block1和parameter_block2相对应的互协方差矩阵的块。
Compute 必须在第一次调用GetCovarianceBlock之前,并且在调用 Compute时,必须在矢量covariance_blocks中出现arameter_block1, parameter_block2 或者 parameter_block2, parameter_block1这两个参数否则,GetCovarianceBlock将返回false。
covariance_block必须指向一个parameter_block1_size × parameter_block2_size大小的内存位置。返回的协方差将是一个行优先矩阵。
GetCovarianceBlockInTangentSpace
返回与parameter_block1和parameter_block2相对应的互协方差矩阵的块。在切空间返回互协方差如果本地参数化与参数块有关,否则就会在环境空间中返回互协方差。
Compute 必须在第一次调用GetCovarianceBlockInTangentSpace之前,并且在调用 Compute时,必须在矢量covariance_blocks中出现arameter_block1, parameter_block2 或者 parameter_block2, parameter_block1这两个参数否则,GetCovarianceBlockInTangentSpace将返回false。
covariance_block必须指向一个parameter_block1_size × parameter_block2_size大小的内存位置。返回的协方差将是一个行优先矩阵。
例子
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)