1 Ceres教程
官方文档 Tutorial——Ceres Solver
简明入门教程 一文助你Ceres 入门
1. 1 构建非线性最小二乘
Ceres 求解器包括两个独立的部分,一个建模的API提供丰富的工具一次性以一条语句构建一个优化问题以及一个求解器API控制求解过程。这个部分主要讲解如何利用Ceres构建优化问题。
Ceres求解如下形式的具有鲁棒的边界约束的非线性最小二乘问题。 在Ceres的语境中,表达式
被称作残差块(residual block), 其中
是代价函数(Cost Function) ,依赖于参数块(parameter blocks)
.
在大多数的优化问题中,经常包含一些标量组。例如平移向量的三个组件,和定义了相机姿态的四元数的四个标量。在Ceres中,这样的一组标量被成为参数块(parameter block).参数块也可以是一个标量。是一个损失函数,一个标量值的函数,用于减少在非线性最小二乘问题求解过程中离群值的影响。
分别是参数块
的上界和下界。
在特殊情况下, 例如当, 并且
时,就得到了最熟悉的无约束最小二乘问题。
CostFunction
对于目标函数中的每个项,一个CostFunction 负责计算一个残差向量和一个Jacobian矩阵。具体地,考虑一个参数块为的代价函数
。对于给定的
, CostFunction 用于计算代价向量
以及Jacobian矩阵。
class CostFunction
{
public:
virtual bool Evaluate(double const* const* parameters,
double* residuals,
double** jacobians) = 0;
const vector<int32>& parameter_block_sizes();
int num_residuals() const;
protected:
vector<int32>* mutable_parameter_block_sizes();
void set_num_residuals(int num_residuals);
};
CostFunction 的签名(输入参数块的大小和输出的大小)分别相应地被保存在CostFunction::parameter_block_size_
和CostFunction::num_residuals
。使用从这个类继承的代码时,需要使用对应的存取器来对两个成员变量进行设置。Problem::AddResidualBlock()
添加这些信息时会对其进行验证。
bool CostFunction::Evaluate(double const* const* parameters, double *residuals, double ** jacobians)
parameters
是一个大小为CostFunction::parameter_block_sizes_.size()
数组的数组。parameters[i]
是一个大小为parameter_block_sizes_[i]
的数组,包含代价函数依赖的第i个参数块。parameters
不能是nullptr.
residuals
是一个大小为num_residuals_
的数组, 也不能是nullptr
.jacobians
是一个大小为CostFunction::parameter_block_sizes_.size()
的数组的数组。
如果jacobians
是nullptr
, 则用户只计算残差项。jacobians[i]
是一个大小为num_residuals ✖ parameter_block_sizes_[i]
的行主数组。
如果jacobians[i]
不是nullptr
, 用户需要计算残差向量residuals
对应参数parameters[i]
的Jacobian矩阵,并保存在这个数组里。类似于
如果jacobians[i]
是nullptr
, 则这个计算会被跳过。这种情况只当对应的参数块被设置为常数时出现。
返回值 用于表征残差项或Jacobians的计算是否成功。这个可以用于传达计算雅可比矩阵过程中的数值故障。
SizedCostFunction
class SizedCostFunction
通常情况下,参数块的大小和残差项的大小在编译时都是已知的。这些值可以作为SizedCostFunction
的模板参数,再使用时只需要对CostFunction::Evaluate().
函数进行重载即可。
template<int kNumResiduals, int... Ns> # 在ceres的头文件中实际支持10个参数块
class SizedCostFunction : public CostFunction{
public:
virtual bool Evaluate(double const* const* parameters,
double *residuals,
double **jacobians) const = 0;
};
LocalParameterization
class LocalParameterization
在许多优化问题中,尤其在传感器融合问题中, 需要对存在于流形空间中的量进行建模,例如一个传感器的旋转和转动用四元数表示。
流形是局部空间类似欧几里德空间的空间。更精确地说,流形上的每个点都有一个和这个流形正切的线性空间,这个线性空间有着和流形内置维度一样的维度,小于或等于流形嵌入的环境空间维度。
例如三维空间中球体上一个点的正切空间是一个与这个球体在该点相切的二维平面。正切空间有两个特点:
一是正切空间均为欧几里德空间,可以应用基本的向量运算。
二是正切空间中的运动可以转化为沿流形的运动,垂直于正切空间的运动不会转化为流形中的运动。回到我们的球体示例,在与球体相切的二维平面上移动并投影回球体将使您远离起始点,但在同一点沿法线移动,而投影回球体将使您回到该点。
除了数学上的精确性外,正确建模流形值量并注意其几何结构也有实际的好处:
一是在优化过程中自然地将值约束在流形上,使得使用者不用进行四元数规范化。
二是减小优化问题的维数到自然的大小。例如限制为直线的量是一维对象,与直线所在的环境空间的尺寸无关。
在正切空间中工作不仅降低了优化算法的计算复杂度,而且改善了算法的数值性能。
流形上的一个基础操作是, 计算沿x处的正切空间中增量移动的结果,然后投影回x所属的流形。也被称为回撤(Retraction),
是欧几里德空间中向量加法的推广。正式地说,
是一个从流形
和 它的正切空间
到流形
的映射,并具有以下恒等式。
这样确保了正切空间以为中心,且零向量是恒等向量。
举两个例子:
欧几里德空间是最简单的流形。它和它的正切空间均具有n维度,
是熟悉的向量加法运算。
另一个更有趣的例子是, 三维特殊正交群,3✖3 旋转矩阵的空间。
是一个三维的流形,嵌入在
空间中。在
上的
通过指数映射实现,从正切空间
到流形。这个指数映射的定义如下所示:
其中 .也即
LocalParameterization
接口允许用户定义并关联参数块所在的流形。通过定义Plus
操作和在
时的导数。
class LocalParameterization{
public:
virtual ~LocalParameterization(){}
virtual bool Plus(const double* x, const double* delta, double* x_plus_delta) const = 0;
virtual bool ComputeJacobian(const double* x, double* jacobian) const = 0;
virtual bool MultiplyByJacobian(const double* x,
const int num_rows,
const double* global_matrix,
double* local_matrix) const;
virtual int GlobalSize() const = 0;
virtual int LocalSize() const = 0;
}
int LocalParameterization::GlobalSize()
参数块所在的环境空间的维度。
int LocalParameterization::LocalSize()
增量所在正切空间的维度。
bool LocalParameterization::Plus(const double* x, const double *delta, double *x_plus_delta) const
实现。
bool LocalParameterization::ComputeJacobian(const double*x, double *jacobian) const
计算Jacobian矩阵,以主行的形式。bool MultiplyByJacobian(const double *x, const int num_rows, const double *global_matrix, double *local_matrix) const
local_matrix = global_matrix * jacobian
global_matrix
是一个num_rows ✖ GlobalSize
的主行矩阵, local_matirx
是一个num_rows ✖ LocalSize
的主行矩阵。 jacobian
是一个由在x处计算的LocalParameterization::ComputeJacobian()
。
这个方法只用在GradientProblem
上, 对于大多数普通的使用场景,使用默认的应用即可。
Ceres 求解器集成了许多LocalParameterization
通常的使用例子。 另外包含一系列关于在流形上进行高质量操作的库是Hauke Strasdat 和他的同事开发的
Sophus
库。
IdentityParameterization
当 和
的维度一致时,
操作
和x在欧几里德流形上的运算是一样的。
QuaternionParameterization
另一个通常出现在从运动恢复结构(SFM)问题中的例子是相机的旋转通常使用四元数进行参数化。这是一个在四维空间中的3维流形。
右边的两个四维向量的乘法的产物是标准的四元数产物。
EigenQuaternionParameterization
Eigen
使用一个不同的内存结构对四元数的元素进行保存。例如, Eigen以的形式保存四元数,实数部分(w)放置在最后一个元素。在通过构造函数构建Eigen 四元数时,元素的参数时按照
的顺序传递的。
由于Ceres 在指向double
型数据的参数块上进行操作,需要进行一些不同的参数化操作。EigenQuaternionParameterization
使用和QuaternionParameterization
一样的Plus
操作,但需要考虑Eigen内置的内存元素顺序。
SubsetParameterization
假设x是一个两维的向量, 用户希望第一位保持常量,是一个标量,
被定义为
SubsetParameterization
将这种结构通过定义需要设置为常量的坐标集合推广至可以保持参数块任何一个部分为常量。
HomogenerousVectorParameterization
在计算机视觉中,齐次向量通常用于表示射影几何中的对象,例如投影空间中的点。在三角化是不佳时,使用这种过度参数化的表示方法是非常有用的。在这里使用其次向量而不是欧几里德向量是有利的,可以表示无穷远或者接近无穷远的点。HomogenerousVectorParameterization
为一个嵌入在n维空间的n-1维的流形定义了一个LocalParameterization
,在这个n维空间中不需要考虑流形中向量的尺度,例如投影空间的元素。假设了n维向量的最后一个坐标是这个齐次向量的尺度部分,例如有限距离的点在这种表示方法中是那些尺度信息不为零的点。更进一步的,
HomogenerousVectorParameterization::Plus
保留了参数的尺度信息。
LineParameterization
这个类为线的参数化提供了方法,使用原点和方向向量进行定义。所以参数向量的大小是环境空间维度的两倍,其中第一部分是原点坐标,第二部分是方向。这种局部参数化的方式是格拉斯曼尼仿射流形的一种特殊情况。
注意这是线条的参数化方法,而不是一个点约束在一条线上。这在对空间中的线条进行优化时很有作用。例如,找到一条线,使得个3D点到这条线的距离之和最短。
ProductParameterization
考虑一个在刚体旋转空间上的优化问题,该空间是一个
和
的笛卡尔积。 假设使用四元数来表示旋转,Ceres 为此提供了一个本地的参数化模块,并且
不需要参数化模块,或使用
IdentityParameterization
进行参数化。我们如何创建一个本地的参数化模块,来代表一个刚体旋转。
在具体案例下,当一个参数块是一系列流形的乘积时,并且具备每个流形的部分参数化模块,ProductParameterization
可以被用于构建一个笛卡尔积的局部参数。对于刚体变换的例子,也即具有一个大小为7的参数块,前四个参数是以四元数表示的旋转,一个局部的参数可以用下面的方式进行构建。
ProductParameterization se3_param(new QuaternionParameterization(),
new IdentityParameterization(3));
Problem
class Problem
保持有界约束非线性最小二乘问题的鲁棒性。使用Problem::AddResidualBlock()
和Problem::AddParameterBlock()
来创建一个最小二乘问题。
例如一个包含3个大小分别为3,4和5的参数块,以及两个残差块大小分别为2和6.
double x1[] = {1.0, 2.0, 3.0};
double x2[] = {1.0, 2.0, 3.0, 5.0};
double x3[] = {1.0, 2.0, 3.0, 5.0, 7.0};
Problem problem;
problem.AddResidualBlock(new MyUnaryCostFunction(...), x1);
problem.AddResidualBlock(new MyBinaryCostFunction(...), x2, x3);
Problem::AddResidualBlock()
如同名字所展示的那样,将一个残差块添加到问题中。这个增加了一个CostFunction
, 一个可选的LossFunction
并将CostFunction
与一系列的参数块连接起来。
代价函数包含的信息包括了期望的参数块大小,会对输入的参数块大小与保存在parameter_blocks
的大小进行比较。如果不相等则会导致程序中止。loss_function
可以是nullptr
,在这种情形下,代价项只是残差的二范数的平方。
用户可以通过Problem::AddParameterBlock()
显式地添加参数块,但这会导致额外的正确性检查,这是因为在调用Problem::AddResidualBlock()
时会隐式地添加参数块如果未添加,所以显式地调用Problem::AddParameterBlock()
并不是必须的。Problem::AddParameterBlock()
显式地添加参数块到Problem
中。 也允许用户将LocalParameterization
对象关联到参数块。 使用相同参数重复调用会被忽略,使用相同的浮点数指针但不同大小会导致未定义行为。
也可以通过使用Problem::SetParameterBlockConstant()
函数将参数设置为常值,通过调用函数SetParameterBlockVariable()
取消这一设置。
事实上,您可以将任意数量的参数块设置为常量,Ceres非常聪明,能够找出您构建的问题的哪一部分取决于可以自由更改且只需花费时间解决的参数块。例如,如果你构造了一个有一百万个参数块和两百万个残差块的问题,然后将除了一个参数块以外的所有参数块都设置为常量,并且说只有10个残差块依赖于这个非常量参数块。如果你定义了一个有一个参数块和10个残差块块的问题,那么Ceres在解决这个问题上花费的计算努力将是相同的。Problem
默认具有cost_function
,loss_function
和local_parameterization
的指针所有权,这些指针在Problem
的生存周期里一直保持着未释放的状态。如果希望对这些对象的析构进行控制,则需要对Problem::Options
中对应的选项进行控制。
即使Problem
具有cost_function
, loss_function
的所有权,但用户依然不具有在其他残差块使用它们的权力。每个cost_function 和loss_function的析构函数只会被调用一次,无论有多少残差块与之关联。**class Problem::Options**
选项结构体用于控制Problem
Ownership Problem::Options::cost_function_ownership
Default: TAKE_OWNERSHIP
。这个选项控制了Problem
对象是否拥有代价函数的所有权。如果设置为占有所有权,则Problem 会在析构函数时删除掉代价函数的指针。这个析构函数只会被删除一次因为允许共享代价函数。Ownership Problem::Options::loss_function_ownership
Default: TAKE_OWNERSHIP
.这个选项控制是否包含损失函数的所有权。Ownership Problem::Options::local_parameterization_ownership
Default: TAKE_OWNERSHIP
,这个选项控制是否包含本地参数的所有权。bool Problem::Options::enable_fast_removal
Default: false
。如果为真,则为更快的Problem::RemoveResidualBlock()
和Problem::RemoveParameterBlock()
操作提供内存交换。
在默认情况下,这两个操作的时间长短和问题的规模成正比。如果移除参数块或残差块的次数较少,这样的时间消耗是可以接受的。因此,如果电脑有足够的内存,可以使能这个选项使得移除参数块的时间仅依赖于该参数块的残差块的数目,而移除残差块的时间是常值。bool problem::Options::disable_all_safety_checks
Default:false
. 默认设置为假,Ceres 在进行问题构建时会进行一系列的安全检查。安全检查会导致一些微量但可观测的性能下降,一般在构建问题的5%左右。如果能够确保问题的构建是正确的,而且5%的时间也是不可容忍的,则可以把这个选项设置为真。Context *Problem::Options::context
Default: nullptr
. 被用于求解问题的Ceres 全局上下文,这有助于减少计算时间,因为Ceresk可以复用复杂的对象。这个上下文对象可以是空指针,在这种情况下,Ceres会创建一个。Ceres 不保有这个指针的所有权。EvaluationCallback *Problem::Options::evaluation_callback
Default:nullptr
.使用这个回调的接口,Ceres会在将要评价残差或者Jacobians矩阵时,通知你。如果存在一个evaluation_callback
, Ceres会根据在调用CostFunction::Evaluate()
使用的值,更新用户的参数块 ,在EvaluationCallback::PrepareForEvaluation()
. 可以利用这个回调函数分享在损失函数之间共享计算,通过在Ceres调用CostFunction::Evaluate()
前在EvaluationCallback::PrepareForEvaluation()
中做分享计算。
ResidualBlockId Problem::AddResidualBlock(CostFunction cost_function, LossFunction loss_function, const vector
template<typename Ts...> ResidualBlockId Problem::AddResidualBlock(CostFunction *cost_function, LossFunction *loss_function, double *x0, Ts...xs)
添加一个残差块到总代价函数。 代价函数中包含这个代价函数中所期望的参数块数量。这个函数会去检查参数块的大小是否与parameter_blocks
中是否保持一致。如果不匹配则会导致程序中止。loss_function
可以是空指针,在这种情况下代价项只是残差项范数的平方。
多个参数块可以通过vector<double*>
或则double*
指针进行传递。
使用范例:
double x1[] = {1.0, 2.0, 3.0};
double x2[] = {1.0, 2.0, 5.0, 6.0};
double x3[] = {3.0, 6.0, 2.0, 5.0, 1.0};
vector<double*> v1;
v1.push_back(x1);
vector<double*> v2;
v2.push_back(x2);
v2.push_back(x1);
Problem problem;
problem.AddResidualBlock(new MyUnaryCostFunction(...), nullptr, x1);
problem.AddResidualBlock(new MyBinaryCostFunction(...), nullptr, x2, x1);
problem.AddResidualBlock(new MyUnaryCostFunction(...), nullptr, v1);
problem.AddResidualBlock(new MyBinaryCostFunction(...), nullptr, v2);
void Problem::AddParameterBlock(double values, int size, LocalParameterization local_parameterization )
添加一个适当大小的参数块到问题中。参数相同的调用会被忽略,同样的doule型指针,但大小不一样,会导致未定义的行为。
void Problem::AddParameterBlock(double *values, int size)
同上
void Problem::RemoveResidualBlock(ResidualBlockId residual_block)
从问题中移除一个残差块。由于Ceres求解器使用引用记数的方法,残差块允许在不同代价函数或损失函数对象之间共享。所以当一个残差块被删除时,残差块中对应的代价函数和损失函数的引用记数会减少,如果减少到零,则会调用析构函数删除该残差块。
如果Problem::Options::enable_fast_removal
为真,则移除的速度很快,几乎是常值。否则是线性增长的,需要对问题进行一个完全的扫描。
移除一个残差块对问题依赖的参数块没有影响。但移除一个残差块或者参数块会更改内置的顺序,使得从求解器返回的jacobian或残差无法解释。如果依赖于评价的jacobian, 不要使用remove。
void Problem::RemoveParameterBlock(const double *values)
从问题中移除一个参数块。任何依赖于此参数块的残差块也会被移除。如果这个参数块的参数化存在,则会存在到问题被删除(析构)。
移除一个残差块或者参数块会更改内置的顺序,使得从求解器返回的jacobian或残差无法解释。如果依赖于评价的jacobian, 不要使用remove。
void Problem::SetParameterBlockConstant(const double *values)
使得指示的参数块在优化过程中保持为常值。
void Problem::SetParameterBlockVariable(double *values)
允许指定的参数块在优化过程中变化。
bool Problem::IsParameterBlockConstant(const double *values) const
查询参数块是否被设置为常量。
void Problem::SetParameterization(double values, LocalParameterization local_parameterization)
为其中的一个参数块设置参数化。local_parameterization
默认属于Problem
。可以为多个参数设置同一个参数化,析构函数会仅删除一次本地参数化。以nullptr
调用这个函数会清除之前设置的所有参数化。
LocalParameterization Problem::GetParameterization(const double values) const
返回与这个参数块关联的局部参数化对象指针,如果没有,则返回空指针。
void Problem::SetParameterLowerBound(double *values, int index, double lower_bound)
为参数块values
对应index
位置的参数设置下界。 默认的下界是-std::numeric_limits<double>::max()
, 在求解器中视为.
void Problem::SetParameterUpperBound(double *values, int index, double upper_bound)
同上。
double Problem::GetParameterLowerBound(const double *values, int index)
获取参数块中指定位置参数的下界。如果未对下界进行设置,则返回
double Problem::GetParameterUpperBound(const double *values, int index)
同上,返回上界。
int Problem::NumParameterBlocks() const
返回问题中参数块的数量。恒等于parameter_blocks.size()
和parameter_block_sizes().size()
int Problem::NumParameters() const
所有参数块的参数向量的大小。
int Problem::NumResidualBlocks() const
返回问题中残差块的数量。恒等于residual_blocks.size()
int Problem::ParameterBlockSize(const double *values) const
参数块的数量大小。
int Problem::ParameterBlockLocalSize(const double *values) const
返回参数块的局部参数化的数量。 如果没有与这个参数块相关联的局部参数化, 则有ParameterBlockLocalSize
= ParameterBlockSize
.
bool Problem::HasParameterBlock(const double *values) const
判断给定的参数是否出现在问题中。
void Problem::GetParameterBlock(vector
使用当前在问题中的参数块的指针填充传递的参数parameter_blocks
向量。 经过这个调用,parameter_block_size() == NumParameterBlocks
.
void Problem::GetResidualBlocks(vector
使用当前问题中的残差块填充传递的参数residual_blocks
向量。经过这个调用,residual_blocks.size()== NumResidualBlocks
.
void Problem::GetParametersBlocksForResidualBlock(const ResidualBlockId residual_block, vector
获取所有依赖给定残差块的所有参数块。
void Problem::GetResidualBlocksForParameterBlock(const double values, vector
获取所有依赖于给定参数块的所有残差块。
const CostFunction *Problem::GetCostFunctionForResidualBlock(const ResidualBlockId residual_block) const
为给定的残差块获取CostFunction
指针。
const LossFunction *Problem::GetLossFunctionForResidualBlock(const ResidualBlockId residual_block) const
为给定的残差块获取LossFunction
指针。
bool EvaluateResidualBlock(ResidualBlockId resdual_block_id, bool apply_loss_function, double cost, double residual, double jacobians) const*
对残差块residuals
中的残差项和残差项和参数的 jacobians [i]进行评价,把标量的代价保存到cost
中。
如果residuals
是nullptr
, 则不计算残差项。
如果jacobians
是nullptr
, 则不计算 Jacobians矩阵, 如果jacobians[i]
是nullptr
,则对应参数块的Jacobian不计算。
如果参数是常量,则无法计算Jacobians.
返回的代价和jacobians 是已经鲁棒化并且也应用了局部参数化,例如使用QuaternionParameterization
定义的4维四元数参数的jacobian 的维度是`num_residuals3,而不是
numresiduals*4。<br />
apply_loss_function选项允许用户决定是否使用损失函数。<br />如果
EvaluationCallBack` 与问题关联了,则在_new_point = true 时,EvaluationCallback::PrepareForEvaluation()
方法会在每次调用该函数时被调用。这是因为假设了用户在上次调用评价或者求解时,参数的值可能发生了改变。为了提高效率,如果已知参数值没有发生改变,可以通过Problem::EvaluateResidualBlockAssumingParametersUnChanged()
来进行评价。
bool EvaluateResidualBlockAssumingParametersUnChanged(ResidualBlockId residual_block_id, bool apply_loss_function, double cost, double residual, double jacobians) const**
与上面的评价函数一样,除非有一个EvaluationCallback
与问题关联。则在new_point = false 时,则对应的EvaluationCallback::PrepareForEvaluation()
方法会在每次调用该函数时被调用。
这意味着有一个EvaluationCallback
与问题关联时,在调用该函数前,需要用户自己去调用EvaluationCallback::PrepareForEvaluation()
函数。这是因为在本函数中,我们假设参数的值没有发生变化。
bool Problem::Evaluate(const Problem::EvaluateOptions &options, double cost, vector
对Problem
进行评价。任何输出的指针都有可能是空指针。使用到的残差块和参数块通过下面所展示的Problem::EvaluateOptions
结构体进行评价。
评价会使用存储在内存位置中的值,该内存位置由构造问题时使用的参数块指针指向,例如在以下代码中:
Problem problem;
double x = 1;
problem.Add(new MyCostFunction, nullptr, &x);
double cost = 0.0;
problem.Evaluate(Problem::EvaluateOptions(), &cost, nullptr, nullptr, nullptr);
如果未使用局部参数化,则梯度向量的大小是所有参数块大小的总和。如果参数块具有局部参数化,则它会向梯度向量提供“LocalSize”条目。
该函数不能在问题求解过程中被调用。
如果一个EvaluationCallback
与问题关联了,则会在调用该函数时调用对应的PrepareForEvaluation 方法当 new_point=true时.
class Problem::EvaluateOptions
选项结构体被用于控制Problem::Evaluate().
vector<double*> Problem::EvaluateOptions::parameter_blocks
应该进行评价的参数块的集合。这个向量决定了参数块在梯度向量和雅可比矩阵列中出现的顺序。如果parameter_blocks
为空,则假定它等于包含所有参数块的向量。一般来说,在这种情况下,参数块的顺序取决于它们添加到问题中的顺序,以及用户是否删除了任何参数块。vector<ResidualBlockId> Problem::EvaluateOptions::residual_blocks
应执行评估的残差块集合。该向量确定残差出现的顺序,以及雅可比矩阵行的排序方式。如果残差块为空,则假定它等于包含所有残差块的向量。bool Problem::EvaluateOptions::apply_loss_function
即使问题中的残差块可能包含损失函数,将apply_loss_function设置为false将关闭损失函数对成本函数输出的应用。例如,如果用户希望通过研究解算前后的残差分布来分析解算质量,则可使用该方法。int Problem::EvaluateOptions::num_threads
评价时使用的线程数量(需要OpenMP)。
EvaluationCallback
class EvauluationCallback
在Ceres 对残差或Jacobians 进行评价时,用于接收回调函数的接口。
class Evaluation{
public:
virtual ~EvaluationCallBack() {}
virtual void PrepareForEvaluation()(bool evaluate_jacobians
bool new_evaluation_point) = 0;
};
void EvaluationCallback::PrepareForEvaluation(bool evaluate_jacobians, bool new_evaluation_point)
Ceres 在计算残差或者jacobians 之前,每次都会先调用这个函数。
用户的参数(也即double 型的数据)在下次调用EvaluationCallback::PrepareForEvaluation().
前都是固定的。如果new_evaluation_point = true
, 则这是一个与上次评价的点不同的点。否则,是与之前评价的点(Jacobian 或者 残差)是一样的,因此用户可以使用上次评估缓存的结果。如果evaluate_jacobian
为真,则Ceres会在后续的代价评估中使用Jacobians。
通过使用这个接口,Ceres会在用户打算对残差或者雅可比矩阵进行评估前通知用户。也可以通过这个接口,在残差块之间共享计算(?),这是通过在Ceres调用CostFunction::Evaluate()
对所有残差块进行评估之前,在EvaluationCallback::PrepareForEvaluation()
中共享计算。也可以通过设置new_evaluation_point
参数实现在进行残差和Jacobian评价前缓存残差的评估结果。
这个回调函数的另一种使用方法是,当代价函数的计算被移到了GPU,这种情况下,这个预计算回调实际计算代价函数的评价,Ceres后续调用真实代价函数只是将GPU中的结果复制到对应的块,以供输入到求解器。
除来自回调函数的通知以外,Ceres没有提供任何其他数据共享机制。用户必须自己设计接口来访问之前已计算好的共享数据。一种方式是在每个代价函数中放置要共享的数据的指针,或者使用一个全局共享变量(*bug-prone)。在Ceres中,这个计算其他代价函数没有什么区别,只是在幕后,代价函数通过使用之前计算的数据,执行得更加快速。
rotation.h
Ceres求解器的许多应用场景都和一些涉及到旋转的变量有关。为了减轻各种旋转表示方式的负担,Ceres提供了一套便捷的模板函数。这些函数是模板,所以用户可以在Ceres求解器的自动求导框架下使用。
template<typename T> void AngleAxisToQuaternion(T const *angle_axis, T *quaternion)
将一个联合的轴角表示方式转化为四元数的形式。值angld_axis
是一个三维的数组,范数是一个以弧度为单位的角度,方向与旋转的轴对齐。quaternion
是一个4维的格子,包含对应的结果。
template<typename T> void QuaternionToAngleAxis(T const *quaternion, T *angle_axis)
将四元数转化为角度。quaternion
的值必须是一个单位四元数-不是规格化的,angle_axis
的值会被设置为范数是弧度形式的角度, 方向是旋转的轴。
template<typename T, int row_stride, int col_stride>
void RotationMatrixToAngleAxis(const MatrixAdapter<const T, row_stride, col_stride> &R,
T *angle_axis)
template<typename T, int row_stride, int col_stride>
void AngleAxisToRotationMatrix(T const* angle_axis,
const MatrisAdapter<T, row_stride, col_stride> &R)
template<typename T>
void RotationMatrixToAngleAxis(T const *R, T *angle_axis)
template<typename T>
void AngleAxisToRotationMatrix(T const *angle_axis, T *R)
角度表示方式和矩阵表示方式的互相转换的函数。将指针指向T而不是MatrixAdapter的函数,假定了输入的矩阵采用列主表示法,单位行跨距和列跨距为3。
template<typename T, int row_stride, int col_stride>
void QuaternionToScaledRotation(const T q[4], const MatrixAdapter<T, row_stride, col_stride> &R)
template<typename T>
void QuaternionToScaleRotation(const T q[4], T R[3*3])
将一个4维向量转化为一个3*3的尺度化的旋转矩阵。在旋转方向的选择方面,四元数会导致一个单位矩阵,而对于微小的
, 四元数
得到以下矩阵。
对应罗德里格近似,最后的矩阵是的叉乘矩阵。同时具备结合率的特性
,定义了从四元数到矩阵的单射。
不执行四元数的标准化,例如, 其中Q是一个正交矩阵,也即有
template<typename T>
void QuaternionToRotation(const T q[4], const MatrixAdapter<T, row_stride, col_stride> &R)
template<typename T>
void QuaternionToRotation(const T q[4], T R[3*3])
和上面的转换函数功能一致,除了旋转矩阵是通过Frobenius 范数规格化以外,也即有.
template<typename T>
void UnitQuaternionRotatePoint(const T q[4], const T pt[3], T result[3])
通过一个四元数q
旋转一个点pt
.
程序假定了输入的四元数是单位范数。如果传递了一个 的四元数,但不会获得使用单位四元数作为参数的结果的两倍。
template<typename T>
void QuaternionRotationPoint(const T q[4], const T pt[3], T result[3])
使用这个函数不用假定q具有单位范数,仅仅假设范数是非零的。
template<typename T>
void QuaternionProduct(const T z[4], const T w[4], T zw[4])
其中*是四元数之间的乘积。
template<typename T>
void CrossProduct(const T x[3], const T y[3], T x_cross_y[3])
template<typename T>
void AngleAxisRotatePoint(const T angle_axis[3], const T pt[3], T result[3])
Cubic Interpolation
优化问题通常与具有一些函数,以表格的形式给出离散的函数值。对这些函数和导数需要对这些值进行插值。对具有离散值的函数进行插值是一个广泛的研究领域,有许多库实现了各种插值方案。但在Ceres的自动求导框架内使用是较为麻烦的。为此Ceres提供了对于一维和二维离散函数进行插值的功能。
一维的插值是基于三维厄尔密样条函数,也被称为Catmull-Rom样条。这可以和获得一阶可微的插值函数。二维插值格式是一维格式的推广,其中假设插值函数在二维中是可分离的。class CubicInterpolator
给定一个无限一维网格作为输入,该网格提供以下接口。
struct Grid1D{
enum{ DATA_DIMENSION = 2; };
void GetValue(int n, double* f) const;
};
此处,GetValue
给出了对于任意整数n和表明函数插值的维度的枚举类型DATA_DIMENSION
.例如对一角度表示的旋转在时间上进行插值,则有DATA_DIMENSION=3
.CubicInterpolator
使用Cubic Hermite 样条曲线来对目标函数进行光滑的近似,可以用于评价在实数范围内的 。例如下面的代码实现了对一个四个数的数组进行插值。
const double x[] = {1.0, 2.0, 5.0, 6.0};
Grid1D<double, 1> array(x, 0, 4);
CubicInterpolator interpolator(array);
double f, dfdx;
interpolator.Evaluate(1.5, &f, &dfdx);
在上述示例中,我们使用Grid1D
作为一个模板类,使得C++的数组和CubicInterpolator
之间可以更好的对接。Grid1D
支持向量值函数,其中函数的各种坐标可以交错或堆叠。也允许使用任何可以转化为双精度的数字类型作为输入。class BiCubicInterpolator
以一个无限的二维网格作为输入,提供了以下的接口。
struct Grid2D{
enum {DATA_DIMENSION = 2;}
void GetValue(int row, int col, double* f) const;
};
其中,GetValue
提供一个插值后的函数在整数对
row
和col
以及表明被插值的函数维度的枚举类型DATA_DIMENSION
下的值。假设对一个具有三通道的彩色图像进行插值(Red, Green &Blue),则DATA_MENSION=3
.BiCubicInterpolator
基于R.Keys提出的一种三次卷积插值算法,可以建立一个光滑近似,计算实数域上任意一点的的值。
下面是对二维数组进行插值的示例代码。
const double data[] ={1.0, 3.0, -1.0, 4.0,
3.6, 2.1, 4.2, 2.0,
2.0, 1.0, 3.1, 5.2};
Grid2D<double, 1> array(data, 0, 3, 0, 4);
BiCubicInterpolator interpolator(array);
double f, dfdr, dfdc;
interpolator.Evaluate(1.5, 2.5, &f, &dfdr, &dfdc);
在以上代码中,模板类Grid2D
被用于将C++
数组为BiCubicInterpolator
构造为类似于一个二维表格。Grid2D
支持主行或主列的格式。支持向量值函数,其中函数的各种坐标可以交错或堆叠。也允许使用任何可以转化为双精度的数字类型作为输入。
1.2 求解非线性最小二乘问题
简介
高效使用Ceres需要熟悉一个非线性求解器的基本组成,在前一部分已经描述了如何使用和配置求解器,后续将简要介绍Ceres中的核心优化算法是如何工作的。
记为一个n维的向量,而
是一个关于
的m维函数,假设要求解如下所示的优化问题。
其中, L 和U分别是参数向量x的上界和下界。
求解的全局最小值是一个棘手的问题,因此我们退而求其次,求解一个局部的最小值。
在下文中,的雅可比矩阵
是一个
的矩阵,其中
,梯度向量
求解非线性优化问题的常用策略是求解原问题的近似队列。在每次迭代过程中,求解建立的近似问题来确定对于向量x的修正.对于非线性最小二乘,可以通过线性化
来构建一个近似,由此可以得到以下的线性最小二乘问题。
简单对问题队列求解,并更新
有可能导致算法不收敛。为了使算法收敛,需要对迭代步长
进行控制。根据
的步长控制方法,非线性优化算法可以分为两种。
- 信任区域。 信任区域方法使用一个模型函数在一个搜索空间的子集即信任区域内对目标函数进行近似。如果模型函数成功地最小化了真实目标函数,则信赖域被扩展;反之,则收缩,需要再次求解模型优化问题。
- 线搜索。 线搜索方式会首先找到一个下降方向,目标函数会沿该方向减少,然后计算步长,该步长决定了沿该方向移动的距离。下降方向可以通过许多方法计算,例如梯度下降法,牛顿法和拟牛顿法。步长可以精确或者不精确地确定。
信任区域方法某种程度上可以说是双线搜索方式,,信任区域方法首先选定一个步长,而后确定一个步长方向。 Ceres 在这两种类别的算法上都是实现了许多算法。
信任区域法
基本的信任区域搜索方法类似于以下的内容。
- 给定一个初始的点x和一个信任区域的半径
.
- 求解。
3.
4. 如果
5. 如果 则,
6. 如果 则,
7. 至第二步。
此处, 是信任区域的半径,
是用于在
的值域上定义的一个度量矩阵,而
用于衡量步长
的质量,例如线性模型预测非线性目标的减少的效果如何。去增加还是减少信任区域的半径取决于线性化反映非线性目标的效果,也即反映在
的值上面。
信任区域法中最关键的计算步骤是求解约束的优化问题。 目前已经有许多种方法用于求解该问题,每种都对应了一种不同的信任区域的算法。目前,Ceres 应用了两种信任区域的算法 Levenberg-Marquardt 和 Dogleg, 如果存在边界约束,则每种方法都对应一种线搜索, 用户可以通过设置
Solver::Options::trust_region_strategy_type
进行设置。
在非线性求解器的层面,块结构是不相关的,因此此处的讨论是一个定义在n维状态向量的优化问题。类似地,当问题在内部转化为纯非线性最小二乘问题时,忽略损失函数的存在。
Levenberg-Marquardt
龙文伯格-马尔库塔算法是求解非线性最小二乘方法的最常用算法,也是最早出现的信任区域方法。Ceres 支持该算法的准确步骤和模糊步骤。
对于式(3)的求解,可以通过以下形式的非约束问题进行求解。
此处,是拉格朗日乘子,是
的逆。在Ceres中,通过求解下式求解结果。
矩阵
是一个正定对角矩阵,通常是矩阵
对角化的平方根。
在进一步讨论问题时,在此做一些标记的简化。假设,矩阵已经连接到矩阵
的底部,同样的,一个零向量被添加到向量
的底部,后续的讨论则是以
和
展开,例如线性最小化问题。
在每次Levenberg-Marquardt算法迭代过程中对于(5) 式最小化问题求解是Ceres的主要计算消耗。Ceres为求解(5)式提供了一系列不同选项。 目前主要有两种方法因式分解和迭代。
因式分解方法是使用Cholesky 或QR分解求得(4) 式的精确解,由此得到一个准确步骤的L-M算法。每次LM的迭代求解(1)式的过程中,(4) 式的精确解的求解是否是必要的不是清晰的。事实上,我们可以发现不必要的案例,例如(4)式本身是(2)式的规则化版本。实际上,我们可以通过线性化问题来构造算法对非线性问题进行近似求解。这些算法被称作非精确牛顿法或者截断牛顿法。
模糊牛顿法通常包含两个内容。首先,一个简单的方法对具有线性方程的系统进行求解。通常采用共轭梯度方法。第二,一个迭代求解器的终止条件。终止条件通常是以下的形式。 此处,
表示L-M 方法的迭代次数,
被乘坐强制队列。Wright Holt 证明了一个基于(6)式使用模糊牛顿迭代方法的L-M算法,会在序列
中收敛且收敛的速率取决于强制序列
的选择。
Ceres同时支持精确和模糊递进的求解策略。当用户选择一种基于因式分解的线性求解器时,Ceres会采用精确递进的L-M;当用户选择一种迭代求解线性求解时,会采用模糊递进L-M算法。
Dogleg
另一种求解(3)式信任区域问题由M. J. D. Powell提出。核心观点是求解以下两个向量。 其中 向量
是(2)式的求解,
是最小化线性近似的向量,如果将我们约束在梯度方向移动。Dogleg 方法通过由求解信任区区域问题得到的
和
来找到向量
。Ceres提供了两种方式来求解信任区域问题,通过设置
Solver::Options::dogleg_type
来选择。TRADITIONAL_DOGLEG
如Powell 所描述的,分别使用Gauss-Newton 和 Cauchy 向量构建两种不同的线段,并在信任区域内查找沿此线最远的点(形似狗腿,这也是名字的来源)。对于推导和计算的更多细节,参加Madsen 等人的文章。SUBSPACE_DOGLEG
是一种更复杂的,包含由两个向量扩展而成的二维子空间,并在这个子空间中找到最小化信任区域问题的点。
在L-M方法上使用Dogleg方法的主要优点是如果以确切的递进计算,在目标函数值的计算上没有获得有效的减少, L-M 方法会以比
更小的值从头开始计算。 Dogleg方法则只需要计算Gauss-Newton 和 Cauchy向量的插值,因为这两者的计算都不依赖
。
Dogleg 方法只适用于基于精确因子分解的线性求解器。
Inner Iterations (内部插值法)
一些非线性最小二乘方法问题在参数块相互作用的方式中具有额外的结构,可以据此对信任区域法的步长计算方法进行改进。例如,考虑如下所示的回归问题。
对于给定的点对, 用户对参数
进行估计。需要注意的是,左边的表达式关于
是线性的,对于给定的
, 可以利用线性回归对
的优化值进行估计。可以从问题中解析地完全将变量
消除。这样的问题被称为可分离最小二乘问题,最著名的算法是Golub 和 Pereyra 发明的变量投影算法。
相同的结构也可以在数据丢失问题中的矩阵分解中找到,对应的算法被称为Wiberg算法。
Ruhe & Wedin 展示了关于求解分离非线性最小二乘法的各种算法的分析,并在变量投影法作为算法一。
实现变量投影是繁琐且耗费计算资源的。Ruhe & Wedin 展示了一种更为简单的算法,具有可比的收敛特性,也即算法II。 算法II 展示了在成功计算一次Newton迭代后,会增加一个额外的优化步骤对 进行估计。
这个思想可以推广到当残差关于 是非线性的。
在这次的案例中,我们为整个问题求解信任区间的步长,并将其作为起始点对 进行优化。对线性化的案例,这相当于执行一个单一的线性化最小二乘求解。对于非线性问题,任何方法都可以对
的优化问题求解。只是约束了
不会同时出现在统一参数块。
进一步推广,不仅仅是优化, 而是将对应于Hessian矩阵稀疏结构的图分解为一组不重叠的独立集,并对每个独立集进行优化。
设置Solver::Options::use_inner_iterations
为true
来使能Ruhe & Wedin 的算法的非线性推广。这个版本的Ceres 有一个更高的迭代复杂度,也在每次迭代上展示了更好的收敛特性。
可以设置Solver::Options::num_threads
为最大的数字来提升计算速度。
非单调步长
基础的信任区域法是一种梯度下降的算法,只在某点减少目标函数时才接受该点。
放宽这个条件使得这个算法可以使得算法在长期内更有效,但代价是目标函数的值局部增加。
这是因为以原则性的方式考虑非递减目标函数值允许算法跳过巨石,因为该方法不限于移动到狭窄的山谷中,同时保持其收敛性。? (可以避免出现局部最小值?)
设置Solver::Options::use_nonmonotonic_steps
为true
来使能Conn, Gould & Toint提出的信任区域算法。
尽管在目标函数的值可能大于优化过程中遇到的最小值,但返回给用户的最终参数是与所有迭代中的最小成本相对应的参数。
所有基于信任区域的方法都可以选择采用非单调步骤。
线搜索方法
Ceres 求解器目前的线搜索方法,无法求解具有边界约束的问题,只能求解无约束的问题。
线搜索的算法。
- 给定一个初始点
- Goto 2
此处是目标函数Hessian的近似,且
是点
处的梯度.根据
的选择,可以获得一系列搜索方向
. 步骤四是一个在
周围一维的优化或者线搜索.
线搜索方式的不同取决于对于搜索方向 和在其附近的一维优化方法的不同. 目前,Ceres求解器支持三种搜索方向,均旨在解决大尺度问题.
1.STEEPEST_DESCENT
这个选择对应的H(x) 为单位阵,这不是一个好的搜索方向,但是是这个问题最简单的形式.此处包含只是为了完整性.
2.NONLINEAR_CONJUGATE_GRADIENT
推广到非线性函数的共轭梯度方法.这个推广以多种形式展开,对应了不同的搜索方向.Ceres求解器目前支持FLETCHER_REEVES
, POLAK_RIBIERE
和HESTENSES_STIEFEL
方向。
BFGS
将割线法推广到多维,获取逆Hessian矩阵的稠密近似,并用于计算拟牛顿迭代。BFGS 是目前已知最通用的拟牛顿算法。LBFGS
一个有限记忆接近完全BFGS
的方法,其中最新的M个迭代被用于近似逆Hessian矩阵,并求解拟牛顿迭代。
目前 Ceres求解器支持一个基于回溯和一个基于插值的Armijo的线搜索方式,以分段/缩放插值 Wolfe条件线搜索算法。然而 为了保证满足BFGS
和LBFGS
的基本假设,应使用Wolfe 线搜索算法。
线性求解器
回到在上述问题中提及的信任区域方法,最关键的计算消耗是对下式最小二乘问题的求解。 使得
及
. 为了便于记号,丢弃关于
的依赖。 则求解(7) 式等价于求解以下的正态方程。
Ceres提供了一系列不同的选项对(8)式求解。
DENSE_QR
对于具有相对稠密的小型问题(数百个参数, 以及上千个残差项),DENSE_QR
可以较好地解决对应的问题。 令 为
的QR分解, 其中
是正交矩阵,
是上三角矩阵。则式(8)的求解通过下式获得。
Ceres使用Eigen
的稠密QR分解方法。
DENSE_NORMAL_CHOLESKY
& SPARSE_NORMAL_CHOLESKY
大型非线性最小二乘问题通常是稀疏的。在这样的例子中,使用稠密的QR分解是效率低下的。记 为正态方程的Cholesky 分解,其中R是一个上三角矩阵,则关于(8)式的解通过下式给出。
细心的读者可能发现在矩阵H的Cholesky 分解中的矩阵和
的QR分解中的矩阵R都是上三角矩阵。由于
是一个对角矩阵,
表明
.这里有两种Cholesky分解的形势——稀疏和稠密两种。
DENSE_NORMAL_CHOLESKY
对正态方程进行稠密Cholesky分解。Ceres 通过使用Eigen
的稠密LDLT分解方法。SPARSE_NORMAL_CHOLESKY
对正态方程进行稀疏Cholesky分解。这可以为大型稀疏问题节省大量时间和内存。Ceres 使用由Tim Davis 教授提出的SuiteSparse
或CXSparse
包提供的稀疏Cholesky分解方法,或者使用Eigen
内置的分解方法(实际上是一个CXSparse
的接口)。
CGNR
对于通常的稀疏问题,如果对于CHOLMOD
而言问题变量维数过大,或者一个稀疏线性代数库没有被链接到Ceres, 可以使用另一个选项是CGNR
求解器。这个求解器使用共轭梯度法对正态方式进行求解,而没有显式构建正态方程。 利用了以下关系。 共轭梯度的方法的收敛是基于条件数
.通常
是条件极差的,需要使用一个
Preconditioner
来获得可行的表现。目前只有JACOBI
预处理器可以在CGNR
中使用。 使用的对角方块线来对正态方程进行预调节。
当用户选择CGNR
作为线性求解器,Ceres会自动从精确步长算法转化为非精确步长算法。
DENSE_SCHUR
& SPARSE_SCHUR
可以使用SPARSE_NORMAL_CHOLESKY
对BA问题进行求解,BA问题具有特殊结构,可以构建更有效的用于求解(8)式的方法。
假设SfM问题包含个相机位姿和
个点, 而且变量向量
.其中
分别对应相机和点的参数。记相机的参数块大小为
,点的参数块的大小为
。(对于大多数问题, c = 6-9, s = 3).Ceres 没有对这些块的大小施加任何恒定要求,将这些值设为常值可以简化说明。
BA问题的关键特性是这里没有项 包含两个或者更多的点块。这也表明矩阵
具有以下的形式。
其中
是一个具有
个大小为
的矩阵块, 而
是一个具有
个大小为
的正对角矩阵。
是一个通常的稀疏矩阵块,对于每个观测值都有一个
%22%20aria-hidden%3D%22true%22%3E%0A%20%3Cuse%20xlink%3Ahref%3D%22%23E1-MJMATHI-63%22%20x%3D%220%22%20y%3D%220%22%3E%3C%2Fuse%3E%0A%20%3Cuse%20xlink%3Ahref%3D%22%23E1-MJMAIN-D7%22%20x%3D%22655%22%20y%3D%220%22%3E%3C%2Fuse%3E%0A%20%3Cuse%20xlink%3Ahref%3D%22%23E1-MJMATHI-73%22%20x%3D%221656%22%20y%3D%220%22%3E%3C%2Fuse%3E%0A%3C%2Fg%3E%0A%3C%2Fsvg%3E#card=math&code=c%20%5Ctimes%20s&id=nvcZG)的矩阵块。假设矩阵的组成。
且
对(8)式进行重新状态设置,并应用高斯消元法。
如上面所提到的,
是一个对角矩阵块,具备较小的对角块的大小
。因此,计算
的逆是相对简单的。可以观察到可以通过
消去
, 由此得到
矩阵
是矩阵
在矩阵
当中的舒尔补。也被成为简化相机矩阵,应为参与到(11) 式当中的变量是对应相机的。
是一个块结构对称正定矩阵,大小为
。 块
对应图像i和j当两幅图像的观察到至少一个共同点是非零的。
可以通过首先构造,用于求解
, 然后回代到
来获得
。从而维度为
的线性系统,被降为求取对角矩阵
的逆,一些矩阵-矩阵、矩阵-向量的乘法,以及求解维数为
的线性系统。对于大多数问题,相机位姿的维数是远小于点的数量的,
, 因此求解(11)是明显比式(10) 要明显更快。这就是舒尔补的技巧。
求解式(11) 依然是有待解决的问题。 求解对称正定矩阵可以通过Cholesky 分解并且依赖矩阵的结构,通常由两种解法。第一种是直接分解,在此方法中,我们保存因子为稠密矩阵。这种方法具有
空间复杂度和
的时间复杂度。只适用于只有数百个相机位姿的问题。 Ceres 将这种策略应用于
DENSE_SCHUR
。
但是一个相对稀疏的矩阵,因为大多数图像只能看到场景的一小部分。这使得我们可以采用第二种选项:稀疏直接法。这些方法将
保存为稀疏矩阵,使用行和列重新排序算法最大化Cholesky分解的稀疏性,并且聚焦于计算分解的非零部分。稀疏直接法,依赖于舒尔补的精确稀疏结构,允许光束平差法实现快速收敛。Ceres 在
SPARSE_SCHUR
求解器应用这种策略。
ITERATIVE_SCHUR
求解另一个BA问题的选项是在简化相机矩阵而不是
上运用预调节共轭梯度。其中一个原因是
的大小比
要小很多,另外也可以证明
。Ceres中在
上应用共轭梯度法的方法为
ITERATIVE_SCHUR
求解器。当用户选择ITERATIVE_SCHUR
作为线性求解器时,Ceres 会自动从精确步长算法切换到非精确步长算法。
使用共轭梯度方法的核心计算操作是对任一向量评价矩阵向量乘积
。可以使用两种方式对乘积进行评价,可以通过
Solver::Options::use_explicit_schur_complement
选项进行设置。根据问题的不同,两种方式的求解性能的差异会较为明显。
- 隐式是默认的选项。隐式评价适用于计算和储存舒尔补
开销较大的大型问题。因为PCG需要通过
和一个向量的乘积来访问
,一种对
的评价方法是观测下式。
我们可以在
上运用PCG,耗费的算力与在
上每次迭代运行PCG相当,同时可以获得一个效果更好的前置调节器。事实上,我们甚至不需要计算矩阵
,(12)可以仅仅通过
的列计算得到。
(12)式和在结构工程和偏微分方程中求解大型线性系统的区域分解法相关度较高。在区域分解的语境下,每个在光束平差问题的点是一个域,而相机组成了域之间的接口。迭代求解舒尔补则可以归属于迭代子结构的技术子类别中。
- 显式评价。矩阵向量乘积的隐式评价的计算复杂度和Jacobian矩阵的非零数的数量是一致的。对于小型或者中型问题,建立舒尔补的代价足够小,可以在内存中直接显式地构造,并用于评价乘积
.
当用户选择ITERATIVE_SCHUR
作为线性求解器时。Ceres
会自动从精确步长算法切换到非精确步长算法。
预调节器
求解(8)式的共轭梯度方法的收敛速度依赖于的特征值的分布。一个可供选择的上界是
,其中
是矩阵
的条件数。对于大多数光束平差问题,
是较大的,并且在(8)式上直接使用共轭梯度方法会导致极差的表现。
求解这个问题的方法是将(8)式替换为一个预调节系统。对于一个给定的线性系统和预调节器
,则得到对应的预调节系统
。对应的算法被称为预调节共轭梯度算法(Preconditioned Conjugate Gradients,PCG),其最差情况的复杂度取决于预调节矩阵
的条件数。
使用预调节器的计算消耗主要是计算
的消耗和评价对于任意向量
的乘积
.需要考虑两个存在冲突的因素:一是
占据
的结构多少使得条件数
尽可能低,二是构建和使用矩阵
的计算消耗.理想的预调节器的条件是
可以达到这样的效果,但不是一个具有可行性的选择,在这样的条件下,等价于求解一个未被预调节的问题.通常来说
包含
的信息越多,耗费的计算资源就越多.例如,基于不完全Cholesky分解的预调节器比基于Jacobi的预调节器具有更好的收敛特性,但通常也耗费更多的资源.
Ceres求解器提供了许多适用于一般稀疏性问题以及光束平差问题中遇到的特殊稀疏性结构的预调节器.
JACOBI
所有调节器中最简单的是对角或者Jacobi 预调节器,例如,任意结构类似矩阵
的块状结构矩阵可以推广到块Jacobi预调节器.Ceres中的
JACOBI
预调节器使用CGNR
选项进行求解时对应矩阵, 使用
ITERATIVE_SCHUR
选项进行求解时对应矩阵的对角块.
SCHUR_JACOBI
另一个对于ITERATIVE_SCHUR
的明显选项是舒尔补的矩阵, 例如
的块Jacobi预调节器.在Ceres当中对应的是
SCHUR_JACOBI
预调节器.
CLUSTER_JACOBI
& CLUSTER_TRIDIAGONAL
对于在从大量图像中重建场景时出现的光束平差问题,可以通过分析和利用场景中的相机点可见性结构来构造更有效的预调节器.
核心思想是根据场景的可视化结构对相机进行聚类.相机信息的相似度可以通过下式得到.
此处,是相机信息i下可视的场景点集合.
Ceres
使用这样的规则来创建CLUSTER_JACOBI
和CLUSTER_TRIDIAGONAL
预调节器.
这两种预调节器的性能取决于在构建预调节器时的聚类算法的速度和聚类质量.Ceres提供了两种可视化聚类算法CANONICAL_VIEWS
和SINGLE_LINKAGE
. 前者是简化的Simon 简化的Canonical Views 算法,后者是经典的Single Linkage 聚类算法.聚类算法的选项通过设置Solver::Options::visibility_clustering_type
进行设置.
SUBSET
这是一个用于求解具有一般稀疏性的问题的预调节器.给定一个问题的残差块的子集,使用Jacobian矩阵的行对应的子集来构建一个预调节器.
假设雅可比矩阵已经被水平分为如下所示.
其中,是
Solver::Options::residual_blocks_for_subset_preconditioner
中残差块中对应的行集合.预调节器是矩阵
预调节器的效率取决于矩阵接近
,或者选择的残差块如何近似于整个问题.
Ordering
在线性求解器中,变量消去的顺序会对方法的效率和精度有着显著影响.例如当进行稀疏Cholesky分解时,一个好的顺序的矩阵排序可以得到一个只占内存的Cholesky 因子,而一个差的顺序会导致一个完全稠密的因子.
Ceres允许用户对变量消除顺序设置一些提示.可以从毫无提示,也即由求解器根据用户的选择自行决定最佳顺序; 到一个确定的顺序,确定变量应该消除的顺序,以及中间许多的可能方法.ParameterBlockOrdering
类的对象被用于传递排序的信息.
形式上,排序是参数块的有序划分.每个参数块属于确切的一个组,每个组都有对应的唯一整数与之对应,这个整数决定了这个组在组的集合中的顺序.我们称这些组为消去组(群).
对于给定了排序的问题,Ceres确保以最小数字对应的消去群会被最先消去,后续按照数字的顺序消去对应的消去群.在每个消去群内部,Ceres可以根据自己的选择来消去参数块.例如,考虑以下的线性系统. 可以有两种方式进行求解,分别是先消去
,再求解
,后回代求解
.或者先消去
,再求解
,再回代求解
因此在使用时可以构建三种排序方式.
1.优先消除
.
2. 优先消除
.
3. 求解器决定消除的顺序.
可以通过把所有变量放置在同一个消去群,让Ceres 在求解时启发式地自动决定顺序.此时,群的编号无关紧要,这和没有排序是一样的.为了给每个变量排序,需要为每个变量创建一个消去群,并按期望的顺序进行排列.
如果用户选择了一种Schur求解器(DENSE_SCHUR
, SPARSE_SCHUR
, ITERATIVE_SCHUR
) 并选择一个特定的顺序,必须有一个重要的特性. 即优先级最高编号最小的消去群必须在Hessian对应的图中是独立的,也即在第一消去群中的任意两个参数块不能同时出现在残差块中.为了达到最佳性能,消去群应该尽可能的大.对于标准光束平差问题,这对应了在第一消去群里包含了所有的3D点,在第二消去群里包含了所有相机参数块.
如果用户将选择权留给Ceres, 则求解器使用近似最大独立集算法来识别第一个消除群.
Solver::Options
class Solver::Options
控制了求解器的所有行为.后续列出对应的设置选项和对应的默认值.
bool Solver::Options::IsValid(string *error) const
判断Options
结构体中的值是不是有效的,如果是则返回真.如果存在问题, 这个方法会返回false 和对应的错误error
,包含了对于原因的描述.
MinimizerType Solver::Options::minimizer_type
默认值:TRUST_REGION
.
在LINE_SEARCH
和TRUST_REGION
两种算法之间进行选择,
LineSearchDirectionType Solver::Options::line_search_direction_type
默认值: LBFGS
有STEEPEST_DESCENT
, NONLINEAR_CONJUGATE_GRADIENT
, BFGS
和LBFGS
.
LineSearchType Solver::Options::line_search_type
默认值:WOLFE
可以在ARMIJO
和WOLFE
(强Wolfe条件)中进行选择.为了保证满足线方向搜索方式BFGS
和LBFGS
的底层假设,需要使用WOLFE
方式进行搜索.
NonlinearConjugateGradientType Solver::Options::nonlinear_conjugate_gradient_type
默认值:FLETHER_REEVES
可供的选项有FLETCHER_REEVES
, POLAK_RIBIERE
以及HESTENES_STIEFEL
.
int Solver::Options::max_lbfgs_rank
默认值: 20
L-BFGS hessian 近似是对于Hessian 矩阵的逆的低阶近似.近似的阶次决定了使用近似的空间和时间复杂度.(泰勒级数展开的阶次?)阶次越高,近似的质量越好.而这种质量的上升受到了限制:
- 该方法仅使用割线信息,不使用实际导数.
- Hessian 近似受到正定的约束.
所以将近似的阶次增加到一个较大的数字,可能会消耗时间和空间复杂度,而没有增加求解质量.没有硬性规定选择最大的阶次.最好的选择通常需要一些特定的问题实验.
bool Solver::Options::use_approximate_eigenvalue_bfgs_scalling
默认值: false
作为BFGS
的更新步骤或LBFGS
右乘步骤,初始的近似Hessian的逆一般是单位阵.然而,Oren的研究表明,使用可以在更大的范围内提高收敛速度,其中
是选择的一个标量用来逼近逆Hessian矩阵的真实特征值.设置
use_approximate_eigenvalue_bfgs_scaling
为真,在BFGS
的第一次迭代前或者LBFGS
的每次迭代使能这个尺度化.
精确的近似特征值的尺度化等于
其中: 其中
是线搜索的目标函数,
是参数值向量.
需要注意的是使用这样的近似特征值尺度化并不会提升收敛速度,反而在一些特定类别的问题中会降低表现,这也是为什么默认设置为false.在问题对于参数变化敏感时,在这种情况下,单个标量因子无法捕捉到这种变化,且导致缩小低敏感部分参数对应的Jacobians近似.也有可能减少对Jacobians残差求解的鲁棒性.
LineSearchIterationType SOlver::Options::line_search_interpolation_type
默认值:CUBIC
用于逼近目标函数的多项式次数.有效值有BISECTION
,QUADRATIC
, CUBIC
.
double Solver::Options::min_line_search_step_size
线搜索会在满足以下条件时终止:
其中对应无穷范数,而
表现了参数值在第k次迭代时的变化.
double Solver::Options::line_search_sufficient_function_decrease
默认值: 1e-4
精确求解线搜索问题是不可行的.幸运的是,基于线搜索方式的算法在有效降低目标函数的值后会返回一个解,从而可以确保收敛而不是一个精确的解. 更为精确的表述是
这个条件被称为 Armijo 条件.
double Solver::Options::max_line_search_step_contraction
默认值:1e-3
在线搜索的每次迭代过程中,
对于收缩:double Solver::Options::min_line_search_step_contraction
默认值:0.6
在线搜索的每次迭代中:
对于收缩:
double Solver::Options::max_num_line_search_direction_restarts
默认值:5
终止优化算法前的最大重新开始线搜索次数.在当前算法无法找到一个新的下降方向时,会重启线搜索.通常意味着数值错误,或者近似的有效性崩溃.
double Solver::Options::line_search_sufficient_curvature_decrease
默认值: 0.9
double SOlver::Options::max_line_search_step_expansion
默认值: 10.0
在Wolfe线搜索的包围阶段,步长会一直增加直到找到一个符合Wolfe条件,或者找到了一个上界包含了满足条件的点.精确的,每次迭代的步长以以下的方式增长.TrustRegionStrategyType SOlver::Options::trust_region_strategy_type
默认值:LEVENBERG_MARQUARDT
可以设置的选项有LEVENBERG_MARQUARDT
和DOGLEG
两种方式.
DoglegType Solver::Options::dogleg_type
默认值:TRADITIONAL_DOGLEG
Ceres 支持TRADITIONAL_DOGLEG
和SUBSPACE_DOGLEG
两种方式.
bool Solver::Options::use_nonmonotonic_steps
默认值:false
放松信任区域法严格递减的步长的约束.
int Solver::Options::max_consective_nonmonotonic_steps
默认值:5
在非单调步长算法的步长选择算法的窗口大小.
int Solver::Options::max_num_iterations
默认值:50
求解器运行的最大迭代次数.
double Solver::Options::max_solver_time_in_seconds
默认值:1e6
求解器运行的最大时间长度
int Solver::Options::num_threads
默认:1
Ceres 评估Jacobians的线程数量.
double Solver::Options::initial_trust_region_radius
默认值:1e4
信任区域的初始大小.使用LEVENBERG_MARQUARDT
策略时,这个数的倒数是初始正则化参数。
double Solver::Options::max_trust_region_radius
默认值: 1e16
信任区域的半径不允许超过这个值。
double Solver::Options::min_trust_region_radius
默认值:1e-32
当信任区域小于这个值时,求解器会终止。
double Solver::Options::min_relative_decrease
默认值: 1e-3
接受信任区域步长之前,相对减少的较低阈值。
double Solver::Options::min_lm_diagonal
默认值: 1e-6
LEVENBERG_MARQUARDT
使用一个对角矩阵对信任区域的步长规则化。这是对角矩阵的值的下界。
double Solver::Options::max_lm_diagonal
默认值: le32
LEVENBERG_MARQUARDT
使用一个对角矩阵对信任区域的步长规则化。 这是对角矩阵的值的上界。
int Solver::Options::max_num_consecutive_invalid_steps
默认值:5
信任区域策略返回的步长有时候是数值无效的,通常是因为条件问题。求解器不停止或中断求解,而是以一个更小的信任区域或更好的调节问题来求解。这个参数被设置为连续的尝试次数在放弃最小化之前。
double Solver::Options::function_tolerance
默认值:1e-6
求解器会在满足以下条件时停止。
其中, 是目标函数值在当前Levenberg-Marquardt法迭代后的变化(上或者下) 。
double Solver::Options::gradient_tolerance
默认值:1e-10
求解器会在满足以下条件时停止。
其中对应无穷范数,
是边界约束下的投影,
是对与参数向量相关的所有局部参数化的加法。
double Solver::Options::parameter_tolerance
默认值:1e-8
求解器会在满足以下条件时停止。
其中是在当前迭代中线性求解器计算的步长。
LinearSolverType Solver::Options::linear_solver_type
默认值:SPARSE_NORMAL_CHOLESKY
/DENSE_QR
在Levenberg-Marquardt 算法的每次迭代中, 用于求解线性最小二乘问题的线性求解器类型。如果Ceres在编译时包含了对SuiteSparses
或CXSparse
或Eigen
的Cholesky 分解的支持,默认值是SPARSE_NORMAL_CHOLESKY
, 否则是DENSE_QR
。
PreconditionerType Solver::Options::preconditioner_type
默认值:JACOBI
迭代线性求解器的预调节器。默认的是块Jacobi预调节器。有效的值(以复杂度增加的顺序)包括IDENTITY
、JACOBI
、SCHUR_JACOBI
、CLUSTER_JACOBI
和CLUSTER_TRIDIAGONAL
。
VisibilityClusteringType Solver::Options::visibility_clustering_type
默认值: CANONICAL_VIEWS
在构建基于可视的预调节器时使用的聚类算法类型。最开始的基于可视的预调节论文和应用,只使用了规范视图算法。
对于大型稠密图优化,这个算法可以给出高质量但开销巨大的结果。这是由于其最差结果的复杂度是图的大小的立方。
另一个选项是使用SINGLE_LINKAGE
, 这是一种简单的阈值单链聚类算法,只关注Schur 补中的紧耦合块。这是一个快速的算法。
对于聚类算法的选择依赖于问题的稀疏性结构,通常来说可以先尝试CANONICAL_VIEWS
, 如果这个方法开销过大,则可以尝试SINGLE_LINKAGE
.
std::unordered_set<ResidualBlockId> resdiual_blocks_for_subset_preconditioner
SUBSET
预调节器适用于具有一般稀疏性的问题。对于给定问题的残差块的子集,使用Jacobian矩阵中的行对应的子集来构建一个预调节器。
假设Jacobian矩阵 已经被水平分为以下的形式。
其中是
Solver::Options::residual_blocks_for_subset_preconditioner
中残差块的对应行集。预调节器是矩阵.
预调节器的效率取决于矩阵近似
的程度,或者选择的残差块如何近似于整个问题。
如果设置了Solver::Options::preconditioner_type == SUBSET
, 则residual_blocks_for_subset_preconditioner
必须是非空的。
DenseLinearAlgebraLibrary Solver::Options::dense_linear_algebra_library_type
默认值:Eigen
Ceres支持多种用于稠密矩阵分解的稠密线性代数库。目前EIGEN
和LAPACK
是有效的选项。LAPACK
对应系统BLAS + LAPACK
库,可能有效也可能无效。
这个设置会影响DENSE_QR
, DENSE_NORMAL_CHOLESKY
和DENSE_SCHUR
求解器。对于中小型问题,Eigen
是一个较好的选择,但对于大型问题,一个优化的LAPACK+BLAS
应用可以在性能上产生巨大差异。
SparseLinearAlgebraLibrary Solver::Options::sparse_linear_algebra_library_type
默认值: 最可靠的线性代数库,根据 SUITE_SPARSE
> CX_SPARSE
> EIGEN_SPARSE
> NO_SPARSE
Ceres支持三种稀疏线性代数库,SuiteSparse
,通过将这个参数设置为SUITE_SPARSE
来使能;CXSparse
设置这个参数为CX_SPARSE
来使能。Eigen
通过设置这个参数为EIGEN_SPARSE
来使能。最后设置为NO_SPARSE
意味着不使用稀疏线性求解器;这与Ceres在编译时是否支持是无关的。SuiteSparse
是一个精密且复杂的稀疏线性代数运算库,是优选考虑的选择。
如果平台或者需求不允许你使用SuiteSparse
, 考虑使用CXSparse
,这是一个更小,更容易构建的库。可以预料的是它在大型问题上的表现是比不过SuiteSparse
的。
最后也可以考虑使用Eigen
库中的稀疏线性代数方法。目前为止,这个库的性能是这三者中最差的,但也许会在后续得到改善。
shared_ptr<ParameterBlockOrdering> Solver::Options::linear_solver_ordering
默认值:NULL
一个排序对象的实体,用于告知求解器应该以何种顺序去消去参数块。
如果为NULL
,则由求解器自行决定最佳的消去顺序。
bool Solver::Options::use_explicit_schur_complement
默认值:false
在ITERATIVE_SCHUR
中显式地计算Schur 补矩阵。
这个选项默认设置为去使能,ITERATIVE_SCHUR
利用Schur补的代数表达式,隐式地对Schur 补和向量的矩阵向量乘积进行评估。
这个评估的成本和Jacobians中非零的数量是成正比的。
对于中小型问题,在一定限度内,计算Schur补的代价足够小,显式地计算Schur补并用于评估矩阵向量乘积会更加有效。
使能这个选项,ITERATIVE_SCHUR
会显式地计算Schur补,这会显著提高ITERATIVE_SCHUR
求解器的性能。
bool Solver::Options::use_post_ordering
默认值: false
稀疏Cholesky 分解算法使用一个填充减少的顺序来排列Jacobian矩阵的列。一般有两种方式。
1.以某种顺序计算Jacobian矩阵,然后使用因子分解算法来对Jacobian矩阵的列。
2.因式分解已经改变列顺序的Jacobian矩阵。
第一选项会占用较大的内存。分解算法需要先复制已经改变顺序的Jacobians矩阵,因此,Ceres预先排列Jacobian矩阵的列,一般来说,这样不会造成性能损失。
极少数的一些例子下,值得使用一个更复杂的重排序算法来获得更好的运行性能,以多余的Jacobian矩阵为拷贝。设置use_post_ordering
为true
来使能这个切换。
bool Solver::Options::dynamic_sparsity
一些非线性最小二乘问题是符号稠密但数值稀疏的,例如在任意给定的状态,只有一小部分的Jacobian项是非零的,但非零项的数目和位置取决于状态。对于这类问题,在每次求解器迭代时对稀疏Jacobian进行分解,而不是对包含所有非零项的矩阵做一个整体的分解。
如果所求解的问题并不具有这样的特性,则最好将这个选项设置为假,否则可能会导致更差的性能。
这个设置只会影响SPARSE_NORMAL_CHOLESKY
求解器。
int Solver::Options::min_linear_solver_iterations
默认值: 0
线性求解器的最小迭代次数。这个选项只在线性求解器是迭代求解器时才有效。例如ITERATIVE_SCHUR
或CGNR
。
int Solver::Options::max_linear_solver_iterations
默认值:500
线性求解器的最大迭代次数。
double Solver::Options::eta
默认值: 1e-1
强制序列参数,截断牛顿求解器使用这个数来控制计算牛顿步长的相对精度。这个参数被传递到ConjugateGradientsSolver
,并在满足下式时终止迭代。
bool Solver::Options::jacobi_scaling
默认值: true
true
表示在将jacobian传递给线性求解器前,Jacobian矩阵通过列范数缩放了。这个提升了正态方程的数值条件。
bool Solver::Options::use_inner_iterations
默认值: true
使用一个简化的变量投影算法的非线性版本。本质上,这相当于使用坐标下降算法对每个牛顿/信赖域步骤进行进一步优化。
注意:内部迭代不能应用于关联了EvaluationCallback
的问题Problem
对象。
double Solver::Options::inner_iteration_tolerance
默认值: 1e-3
通常来说,内部迭代会在求解的早期过程具有显著的改善,但在后期会迅速下降,在这个时候,进行内部迭代的时间花费已经不值得了。
当由于内部迭代在目标函数中引起的相对减少小于inner_iteration_tolerance
时,在后续的信任区域最小化中使用内部迭代是不被允许的。
shared_ptr<ParameterBlockOrdering> Solver::Options::inner_iteration_ordering
默认值:NULL
如果Solver::Options::use_inner_iteration
为真,则用户具有两个选择。
1.让求解器启发式地决定在每次内部迭代时优化哪些参数块。为此,将Solver::Options::inner_iteration_ordering
为NULL
.
2.特殊化一个排列好的独立的集合。编号较低的群在内部优化过程中会在更高编号的群之前被优化。每个群都必须是独立的集合。不是所有参数块都需要在排序中被包含。
LoggingType Solver::Options::logging_type
默认值: PER_MINIMIZER_ITERATION
bool Solver::Options::minimizer_progress_to_stdout
默认值: false
默认情况下,Minimizer
过程的输出根据vlog
等级被打印到STDERR
。如果这个标志被设置为真,而且Solver::Options::logging_type
不是SILENT
, 输出的log会以STDOUT
的形式打印。
对于TRUST_REGION_MINIMIZER
过程,展示的结果如下。
iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time
0 4.185660e+06 0.00e+00 1.09e+08 0.00e+00 0.00e+00 1.00e+04 0 7.59e-02 3.37e-01
1 1.062590e+05 4.08e+06 8.99e+06 5.36e+02 9.82e-01 3.00e+04 1 1.65e-01 5.03e-01
2 4.992817e+04 5.63e+04 8.32e+06 3.19e+02 6.52e-01 3.09e+04 1 1.45e-01 6.48e-01
此处。cost
是目标函数的值。cost_change
在接受这次迭代后,目标函数值的变化。|gradient|
梯度的无穷范数。|step|
参数向量中的变化。tr_ratio
目标函数值的实际变化与信任区域模型的值的变化之比。tr_radius
信任区域半径的大小。ls_iter
用于计算信任区域步长的线性求解器迭代次数。对于基于直接/分解的求解器,这个值总是1,对于像ITERATIVE_SCHUR
的迭代求解器,是共轭梯度算法的迭代次数。iter_time
当前迭代耗费的时间。total_time
最小化耗费的总时长。
对于LINE_SEARCH_MINIMIZER
来说,过程展示如下。
0: f: 2.317806e+05 d: 0.00e+00 g: 3.19e-01 h: 0.00e+00 s: 0.00e+00 e: 0 it: 2.98e-02 tt: 8.50e-02
1: f: 2.312019e+05 d: 5.79e+02 g: 3.18e-01 h: 2.41e+01 s: 1.00e+00 e: 1 it: 4.54e-02 tt: 1.31e-01
2: f: 2.300462e+05 d: 1.16e+03 g: 3.17e-01 h: 4.90e+01 s: 2.54e-03 e: 1 it: 4.96e-02 tt: 1.81e-01
f
是目标函数的值。d
在接受这次迭代后,目标函数值的变化。g
梯度的无穷范数。h
参数向量中的变化。s
是通过线搜索计算的优化的步长。it
当前迭代耗费的时间。tt
最小化耗费的总时长。
vector<int> Solver::Options::trust_region_minimizer_iterations_to_dump
默认值:empty
信任区域问题求解时崩溃的迭代编号列表。用于测试和基准测试,如果为empty
, 则问题求解的迭代过程没有一次问题.
string Solver::Options::trust_region_problem_dump_directory
默认值:\tmp
问题保存的路径。如果Solver::Options::trust_region_minimizer_iterations_to_dump
是非空的,且Solver::Options::trust_region_problem_dump_format_type
没有被设置为CONSOLE
,则这个路径应该是非空的。
DumpFormatType Solver::Options::trust_region_problem_dump_format
默认值:TEXTFILE
当Solver::Options::trust_region_minimizer_iterations_to_dump
为真时,信任区域问题的输出格式。通常有三种方式。CONSOLE
以可读的格式输出线性最小二乘问题到stderr
. Jacobian矩阵以稠密矩阵的形式打印。向量.这应该只用于一些小型问题。
TEXTFILE
按照Solver::Options::trust_region_problem_dump_directory
指定的路径输出最小线性二乘问题的文本文件,可以通过MATLAB/Octave
来读数据。Jacobian矩阵是包含三元数组构成的文本文件,向量
被存储为包含一个值的列表的文本文件。
同时会输出一个MATLAB/Octave
脚本,一般名为ceres_solver_iteration_???.m
,用于将问题格式化并加载到内存中。
bool Solver::Options::check_gradients
默认值:false
用有限差分检测每个残差块计算的所有Jacobian矩阵。因为需要计算常用方法计算导数(例如用户指定的,自动求导的),然后计算有限差分,这通常会耗费大量计算。而后会对结果进行比较,如果差别很大,则优化会失败,细节会被保存在求解器的摘要里。
double Solver::Options::gradient_check_relative_precision
默认值: 1e-8
梯度检测器中的检测精度,如果Jacobian矩阵中的某个元素的相对差值超过了这个值,则该代价项的Jacobian将被转储。
double Solver::Options::gradient_check_numeric_derivative_relative_step_size
默认值:1e-6
当Solver::Options::check_gradients
被设置为真时。对于有限差分,每个维度都通过一些有轻微偏移的值进行评价,对于前向差分,数值导数如下。 有限差分是在每个维度进行的。使用一个相对步长而不是绝对步长的原因是,数值微分适用于参数较大的函数(10e9),当参数值较小时(10e-5),有可能构建出打破启发式有限差分的案例。
bool Solver::Options::update_state_every_iteration
默认值:false
如果update_state_every_iteration
设置为true
, Ceres 求解器会在每次迭代后,在用户调用IterationCallback
前,参数块会被更新到求解找到的最佳结果。这是因为IterationCallback
可以检查参数块的值,以便进行计算、可视化或终止。
如果update_state_every_iteration
为false
, 则这里不会有对应的保证,用户也不应该提供IterationCallback
对参数块进行查看或者解释这些变量的值。
vector<IterationCallback> Solver::Options::callbacks
Callbacks 会在Minimizer
每次迭代的末尾调用。以在这个向量中的顺序被调用。
参数块默认只在优化的末尾被调用,例如当Minimizer
终止时。这意味着在默认情况下,如果一个IterationCallback
检查参数块,无法观察到在优化过程中的变化。
设置Solver::Options::update_state_every_iteration
为真,来告知Ceres 在每次迭代末尾并在调用用户的回调函数之前更新这些参数块。
ParameterBlockOrdering
class ParameterBlockOrdering
是一个用于存储和操作具有以下语义的有序组/集合的类:
群的ID是非负整型值。元素是可以用作映射的键或集合元素的任何类型。
一个元素只能归属于一个群。一个群可以包含任意数量的元素。
群通过群id进行排序。
bool ParameterBlockOrdering::AddElementToGroup(const double *element, const int group)
添加一个元素到群。如果不存在具有这个id的群,则会创建一个。可以对同样的元素多次调用该方法。群id应该是非负数。返回值表明添加这个值是否成功。
void ParameterBlockOrdering::Clear()
清除排序。
bool ParameterBlockOrdering::Remove(const double *element)
移除对应的元素,无论该元素属于哪个群。如果这个元素不属于任何群,则调用该函数会导致崩溃,返回值表明这个元素是否真的被移除了。
void ParameterBlockOrdering::Reverse()
反转群的顺序。
int ParameterBlockOrdering::GroupId(const double *element) const
返回元素对应的群的id. 如果这个元素不是任意一个群的成员,返回-1.
bool ParameterBlockOrdering::IsMember(const double *element) const
返回真如果存在群包含对应的参数块。
int ParameterBlockOrdering::GroupSize(const int group) const
这个函数返回总是存在的,例如隐式地存在一个包含所有整数的群。
int ParameterBlockOrdering::NumElements() const
进行排序的元素数量。
int ParameterBlockOrdering::NumGroup() const
非空的群的数量。
IterationSummary
class IterationSummary
这个类描述了在每次迭代后minimizer的状态。
int IterationSummary::iteration
当前迭代的次数。
bool IterationSummary::step_is_valid
步长是数值有效,例如,所有的值都是有限的,并且这个步长减少了线性化模型的值。 在迭代序号为0时,这个变量的值是false
.
bool IterationSummary::step_isnonmonotonic
步长没有有效减少目标函数的值,但由于使用了非单调信任区域算法而被接受了。在迭代序号为零时,IterationSummary::step_is_nonmonotonic
为假的。
bool IterationSummary::step_is_successful
最小化器是否接受了这个步长。如果使用普通信任区域算法,这意味着目标函数值的相对减少比Solver::Options::min_relative_decrease
. 如果使用了非单调信任区域算法,也即Solver::Options::use_nonmonotonic_steps = true
,则即使相对减少是无效的,这个算法也会接受这个步长,这个步长也被声明是成功的。当IterationSummary::iteration=0
时,IterationSummary::step_is_successful
是false
。
double IterationSummary::cost
目标函数的值。
double IterationSummary::cost_change
目标函数的值在这次迭代中的改变。可以是正也可以是负。
double IterationSummary::gradient_max_norm
梯度向量的无穷范数。
double IterationSummary::gradient_norm
梯度向量的二范数
double IterationSummary::step_norm
这次迭代的步长的二范数。
double IterationSummary::relative_decrease
对于信任区域算法, 实际代价函数变化值与线性化近似代价函数变化值的比率。 在线搜索中不使用这个域。
double IterationSummary::trust_region_radius
在当前迭代最后的信任区域的大小。对于L-M算法,这个规则化参数是1.0 / member::IterationSummary::trust_region_radius.
double IterationSummary::eta
对于非精确步长的L-M算法,这是步长求解的相对经典。这个数字仅适用于可以不精确求解系统的迭代解算器。基于因式分解的精确的求解器的eta的值一般为0.0。
double ItrationSummary::step_size
通过线搜索方法计算得到的步长。当使用信任区域求解时,不具有这个项目。
int IterationSummary::line_search_function_evaluations
行搜索算法使用到的函数求值数量。
int IterationSummary::linear_solver_iterations
线性求解器求解信任区域步长的迭代次数。
线搜索方法不具有这个选项。
double IterationSummary::iteration_time_in_seconds
在当前迭代中的最小化过程中耗费的时间。
double IterationSummary::step_solver_time_in_seconds
信任区域法步长求解花费的时间。
double IterationSummary::cumulative_time_in_seconds
自用户调用Solve()
以来的时间。
IterationCallback
class IterationCallback
用于在最小化的每次迭代末尾指定回调函数的结构。
class IterationCallback
{
public:
virtual ~IterationCallback() {}
virtual CallbackReturnType operator() (const IterationSummary& summary) = 0;
};
求解器使用operator()
的返回值来决定是继续求解还是终止。用户可以返回以下三个值。
1.SOLVER_ABORT
表明回调函数检测到了一个异常情况。求解器返回前没有更新参数块(Solver::Options::update_state_every_iteration
被设置为真)。 此时求解器返回的Solver::Summary::termination_type
被设置为USER_FAILURE
。
2.SOLVER_TERMINATE_SUCCESSFULLY
表明没有必要继续优化(符合了一些用户设置的终止条件)。求解器返回地Solver::Summary::termination_type
为USER_SUCCESS
。
3.SOLVER_CONTINUE
表明求解器需要继续迭代优化。
以下是Ceres内置的用于输出优化过程记录的回调函数IterationCallback
class LoggingCallback : public IterationCallback {
public:
explicit LoggingCallback(bool log_to_stdout)
: log_to_stdout_(log_to_stdout) {}
~LoggingCallback() {}
CallbackReturnType operator()(const IterationSummary& summary) {
const char* kReportRowFormat =
"% 4d: f:% 8e d:% 3.2e g:% 3.2e h:% 3.2e "
"rho:% 3.2e mu:% 3.2e eta:% 3.2e li:% 3d";
string output = StringPrintf(kReportRowFormat,
summary.iteration,
summary.cost,
summary.cost_change,
summary.gradient_max_norm,
summary.step_norm,
summary.relative_decrease,
summary.trust_region_radius,
summary.eta,
summary.linear_solver_iterations);
if (log_to_stdout_) {
cout << output << endl;
} else {
VLOG(1) << output;
}
return SOLVER_CONTINUE;
}
private:
const bool log_to_stdout_;
};
CRSMatrix
class CRSMatrix
一种压缩行稀疏矩阵,主要用于向用户传递雅可比矩阵。int CRSMatrix::num_rows
行数int CRSMatrix::num_cols
列数vector<int> CRSMatrix::rows
CRSMatrix::rows
是一个CRSMatrix::num_rows
+1 大小的数组,指向CRSMatrix::cols
和CRSMatrix::values
数组。vector<int> CRSMatrix::cols
包含了和矩阵中的非零项一样多的项。对于每一行i
, cols[rows[i]]
….cols[rows[i+1] -1]
是行i
中非零项的列。vector<int> CRSMatrix::values
包含了矩阵中的所有非零项,对于每一行i
,values[rows[i]]
… values[rows[i+1] -1 ]
是行i
的非零索引。
对于 稀疏矩阵。
0 10 0 4
0 2 -3 2
1 2 0 0
则对应的三个数组是
-row0- ---row1--- -row2-
rows = [ 0, 2, 5, 7]
cols = [ 1, 3, 1, 2, 3, 0, 1]
values = [10, 4, 2, -3, 2, 1, 2]
Solver::Summary
class Slover::Summary
在终止之后,求解器的各个过程的总结。
string Solver::Summary::BriefReport() const
求解器终止后,关于求解器的一句话的简要描述。
string Solver::Summary::FullReport() const
终止后解算器状态的完整多行描述。
bool Solver::Summary::IsSolutionUsable() const
优化算法返回的值是否数值可靠。如果求解器的Solver::Summary::termination_type
被设置为CONVERGENCE, 则这个值可能是USER_SUCCESS
或NO_CONVERGENCE
.
MinimizerType Solver::Summary::minimizer_type
使用的最小化算法类型。
TerminationType Solver::Summary::termination_type
minimizer 终止的原因。
string Solver::Summary::message
求解器的终止原因。
double Solver::Summary::initial_cost
在优化前,问题的代价函数值。
double Solver::Summary::final_cost
优化后,问题的代价函数值。
double Solver::Summary::fixed_cost
由预处理器固定的残差块产生的代价函数,残差块固定的原因是因为其依赖的参数块是固定的。
vector<IterationSummary> Solver::Summary::iterations
以迭代次数排列的IterationSummary
。
int Solver::Summary::num_successful_steps
步长接受的迭代次数。 除非设置了Solver::Options::use_non_monotonic_steps
为真,否则这也是目标函数值下降的次数。
int Solver::Summary::num_inner_iteration_steps
内置迭代运行的次数。
double Solver::Summary::preprocessor_time_in_seconds
预处理器花费的时间。
double Solver::Summary::minimizer_time_in_steps
在Minimizer 中花费的时间。
double Sover::Summary::posetprocessor_time_in_seconds
在后处理器中花费的时间。
double Solver::Summary::total_time_in_seconds
在求解器中花费的时间。
double Solver::Summary::linear_solver_time_in_seconds
线性求解器中计算信任区域步长的时间。
int Solver::Summary::num_linear_solves
通过求解线性系统计算Newton 步长的次数。这个不包含内部迭代使用的线性求解。
double Solver::Summary::residual_evaluation_time_in_seconds
评估残差向量的时间消耗。
int Solver::Summary::num_residual_evaluations
残差项评估的次数。
double Solver::Summary::jacobian_evaluation_time_in_seconds
评估Jacobian矩阵花费的时间。
int Solver::Summary::num_jacobians_evaluations
Jacobian 和残差的估计次数。
double Solver::Summary::inner_iteration_time_in_seconds
内部迭代耗费的时间。
int Solver::Summary::num_parameter_blocks
问题中参数块的数量。
int Solver::Summary::num_parameters
问题中参数的数量。
int Solver::Summary::num_effective_parameters
问题的正切空间的维数(或这个问题的Jacobian矩阵的列数)。如果参数块和LocalParameterization
有关联,那这与Solver::Summary::num_parameters
是有所不同的。
int Solver::Summary::num_residual_block
问题中残差块的数量。
int Solver::Summary::num_residuals
问题中残差的数量。
int Solver::Summary::num_parameter_bloc_reduced
移除非活跃和常量参数块后,问题中的参数块数量。如果一个参数块没有任何残差块与之相关,则是非活跃的。
int Solver::Summary::num_residuals_reduced
精简后的残差数量。
int Solver::Summary::num_threads_given
给定的线程数量用于计算Jacobian和残差评估。
int Solver::Summary::num_thread_used
Ceres 进行Jacobian 和残差评估时使用的实际线程数量。这个数字与Solver::Summary::num_threads_given
并不一致,如果没有可用的OpenMP 或CXX_THREAD。
LinearSolverType Solver::Summary::linear_solver_type_given
用户给定的线性求解
LinearSolverType Solver::Summary::linear_solver_type_used
实际使用的线性求解器类型。可能与Solver::Summary::linear_solver_type_given
不同,如果Ceres确定了问题结构和线性求解器请求的不匹配或用户请求的线性求解器不存在。
vector<int> Solver::Summary::linear_solver_ordering_given
给定的消去群的顺序。
vector<int> Solver::Summary::linear_solver_ordering_used
实际使用的消去群的顺序。
std::string Solver::Summary::schur_structure_given
对于Schur类型的线性求解器,这个字符串是在问题中实例化的模板名称。
std::string Solver::Summary::schur_structure_used
实际使用的Schur 结构。
bool Solver::Summary::inner_iterations_given
如果用户在子优化过程中请求使用内部迭代,则为true
bool Solver::Summary::inner_iteration_used
如果问题实际求解过程中确实使用了内部迭代。
vector<int> inner_iteration_ordering_given
用户给定的内部迭代的参数群顺序。
vector<int> inner_iteration_ordering_used
实际使用的内部迭代的参数群顺序。
PreconditionerType Solver::Summary::preconditioner_type_given
用户指定的预调节器的类型。
PreconditionerType Solver::Summary::preconditioner_type_used
实际使用的预调节器的类型。
VisibilityClusteringType Solver::Summary::visibility_clustering_type
基于可视性的预调节器使用的聚类算法类型。只有当Solver::Summary::preconditioner_type
设置为CLUSTER_JACOBI
或者CLUSTER_TRIDIAGONAL
才有意义。
TrustRegionStrategyType Solver::Summary::trust_region_strategy_type
使用的信任区域类型,
DoglegType Solver::Summary::dogleg_type
求解信任区域问题的dogleg 策略的类型。
DenseLinearAlgebraLibraryType Solver::Summary::dense_linear_algebra_library_type
使用的稠密线性代数库的类型。
SparseLinearAlgebraLibraryType Solver::Summary::sparse_linear_algebra_library_type
使用的稀疏线性代数库的类型。
LineSearchDirectionType Solver::Summary::line_search_direction_type
线搜索方向的类型。
LinearSearchType Solver::Summary::line_search_type
线搜索方法的类型。
LineSearchInterpolationType Solver::Summary::line_search_interpolation_type
当使用线搜索时,用于近似目标函数的多项式阶次。
NonLinearConjugateGradientType Solver::Summary::nonlinear_conjugate_gradient_type
如果线搜索的方式为NONLINEAR_CONJUGATE_GRADIENT
, 则这个表示了具体的非线性共轭梯度的种类。
int Solver::Summary::max_lbfgs_rank
如果线搜索类型的方向是LBFGS, 则这个表明了Hessian 近似的阶次。
N Bugs & Problems
- 求解时出现报错,暂无解决办法