为了描述简便,用函数 $f(x)$ 来表示我们的目标函数,其中 $x$ 是包含了我们所有参数的向量,现在我们的目标就成了求解 $f(x)$ 在何处取得极小值。这个函数的参数这么多,很明显不能直接一步到位,但是从某个值出发(对于三维重建的问题,这个出发点已经有了,就是前面计算出来的相机姿态和三维点坐标)慢慢向正确的值靠近,一次靠近一丢丢($\delta x$),很多次后就可以很接近或者达到正确的值了。于是嘛,问题又变了,变成了我们一次靠近多少的问题了。。。问题变到这儿,貌似就有现成的方法可以解决啦。

梯度下降(最速下降)法

从名字就可以看出来,这个方法是用了梯度信息。梯度的方向是函数上升最快的方向,我们要求的是最小值,那就沿着梯度方向的反方向进行靠近就可以了嘛。于是梯度下降法就是每次沿着梯度反方向以一个给定的步长进行逼近,这样就可以一步一步的靠近最小值。

这个方法虽然可以保证每一次都是在下降,但是下降的速度却并不像它的名字那样很快,收敛的速度较慢。

牛顿法

这篇文章 很好的解释了牛顿法的原理,但是貌似只能求零点啊,这个求极值有啥关系,函数导数为0的地方不就是极值么。。假设在第$k$ 次迭代的时候,我们的参数向量为$x_k$ ,本次迭代的增量为 $\delta x_k$ (也就是本次迭代要求解的东西),那么函数就变成了$f(x_k+\delta x_k)$ ,首先将这个函数进行泰勒展开到二次项:
$$
f(x_k+\delta x_k) \approx f(x_k) + \delta x_k f’(x_k) +\frac{1}{2} (\delta x_k)^2 f’’(x_k) = \phi(\delta x_k)
$$
然后对$\phi(\delta x_k)$ 求导:
$$
\phi ‘(\delta x_k) = f’(x_k) + f’’(x_k)\delta x_k
$$
令上面的导数为0就可以求得$\delta x_k$ ,但是需要注意的是,$x_k、\delta x_k$ 都是向量,所以$f’(x_k)、f’’(x_k)$ 都是矩阵,它们分别叫做雅可比(Jacobi)矩阵(记为$J$)和海森(Hessian)矩阵(记为$H$),雅克比矩阵是函数对所有参数求一阶偏导得到,海森矩阵是函数对所有参数求二阶偏导得到。于是:
$$
J+H\delta x_k = 0 \
\delta x_k = -H^{-1}J
$$
这种方法下降速度比梯度下降要快很多(至于为啥快,我就不知道了,,),但是他不能保证每次迭代都是在下降的(但是总体是在下降的)。这样貌似问题解决了,但是,回头去看看我们的函数,要对这个函数求二阶偏导,简直太复杂了,大家选择放弃,,,,

高斯牛顿法

牛顿大神把问题解决了一半,高斯大神决定再帮他一把,解决海森矩阵的问题,至于解决方式嘛,,,,就是不求它。我们在泰勒展开的时候就只展开到一次项:
$$
f(x_k+\delta x_k) \approx f(x_k)+ J\delta x_k
$$
这样貌似等式右边就变成了类型$ax+b$ 这种线性的形式了啊,这懒偷得有点多了,必须找其他的约束补回来,,,比如我们的优化的是重投影误差,我们最希望它为零啊,,,于是嘛就直接令上式为零就得到了:
$$
f(x_k) + J\delta x_k = 0
$$
貌似$\delta x_k​$ 就可以求出来了:$\delta x_k = -J^{-1}f(x_k)​$ ,但是总感觉有啥不对,,,恩,$J​$ 的大小是$m\times n​$ 的,它可能不是一个方阵!!!所以它是不能直接求逆的,所以需要按照这样算:
$$
f(x_k) + J\delta x_k = 0 \
J^Tf(x_k) + J^TJ\delta x_k = 0 \
\delta x_k = -(J^TJ)^{-1}J^Tf(x_k)
$$

这里可以看出来高斯牛顿法是利用了我们目标函数的特性来使用雅克比矩阵近似海森矩阵。这里只是解决了不计算海森矩阵,但是它还是存在不能保证每次迭代都是在下降的问题。

Levenberg-Marquard算法

梯度下降可以保证每次迭代都是下降的,但是收敛慢,高斯牛顿收敛快但是不能保证每次迭代都是下降的。于是LM算法就出现啦,它把梯度下降和高斯牛顿法结合了起来,使得每次迭代都是下降的,而且收敛速度还比较快。至于怎么结合的,简直不能更粗暴,用一个参数$\lambda$ 来控制增量计算方式是靠近高斯牛顿迭代的方法还是靠近梯度下降的方法,它把上面计算增量的方程改成了下面的形式:
$$
\delta x_k = -[J^TJ+\lambda diag(J^TJ)]^{-1} J^Tf(x_k)
$$
可以看出来$\lambda$ 很小甚至为零时,增量计算方式和高斯牛顿迭代的方法是一致的;当$\lambda$ 很大的时候,大到$J^TJ$ 可以忽略的时候,就变成了:
$$
\delta x_k = -\lambda J^Tf(x_k)
$$
恩,长得和梯度下降的式子挺像的。目前的Bundle Adjustment貌似都是基于LM算法实现的。

雅克比矩阵

上面提到了很多次雅克比矩阵,甚至于每次迭代都需要计算一次雅克比矩阵,雅克比又是个啥?对于我们三维重建的问题,我们期望优化的参数有相机姿态$[R\ \ \ t]$ ,这里一共有6个未知数(为啥是6个?因为旋转矩阵 $R$ 的自由度为3,事实上也可以用罗德里格斯公式转化为一个三维的向量$(r_1,r_2,r_3)$,位移向量$t$ 是长度为3的向量$(t_1,t_2,t_3)$),三维点的坐标$(x,y,z)$ 也就是一共就,9个参数。那么雅克比矩阵就一共有9列,每一列就是函数$f(x)$ 对不同的参数进行求偏导,这一列的元素值就是这个偏导数在一组参数下的取值:
$$
J = [\frac{\partial f}{\partial r_1} \ \ \frac{\partial f}{\partial r_2} \ \ \frac{\partial f}{\partial r_3} \ \ \frac{\partial f}{\partial t_1} \ \ \frac{\partial f}{\partial t_2} \ \ \frac{\partial f}{\partial t_3} \ \ \frac{\partial f}{\partial x} \ \ \frac{\partial f}{\partial y} \ \ \frac{\partial f}{\partial z}]
$$
假设一共有m个三维点,n幅图像,那么雅克比矩阵就应该有$i\times j$行,总的雅克比矩阵可以:
$$
J = \left[\begin{matrix}
\frac{\partial f}{\partial r_1^1} &\frac{\partial f}{\partial r_2^1}&\frac{\partial f}{\partial r_3^1} &\frac{\partial f}{\partial t_1^1} &\frac{\partial f}{\partial t_2^1} &\frac{\partial f}{\partial t_3^1} &\frac{\partial f}{\partial x^1} &\frac{\partial f}{\partial y^1} &\frac{\partial f}{\partial z^1} \
\frac{\partial f}{\partial r_1^2} &\frac{\partial f}{\partial r_2^2}&\frac{\partial f}{\partial r_3^2} &\frac{\partial f}{\partial t_1^2} &\frac{\partial f}{\partial t_2^2} &\frac{\partial f}{\partial t_3^2} &\frac{\partial f}{\partial x^1} &\frac{\partial f}{\partial y^1} &\frac{\partial f}{\partial z^1} \
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \
\frac{\partial f}{\partial r_1^n} &\frac{\partial f}{\partial r_2^n}&\frac{\partial f}{\partial r_3^n} &\frac{\partial f}{\partial t_1^n} &\frac{\partial f}{\partial t_2^n} &\frac{\partial f}{\partial t_3^n} &\frac{\partial f}{\partial x^1} &\frac{\partial f}{\partial y^1} &\frac{\partial f}{\partial z^1} \
\frac{\partial f}{\partial r_1^1} &\frac{\partial f}{\partial r_2^1}&\frac{\partial f}{\partial r_3^1} &\frac{\partial f}{\partial t_1^1} &\frac{\partial f}{\partial t_2^1} &\frac{\partial f}{\partial t_3^1} &\frac{\partial f}{\partial x^2} &\frac{\partial f}{\partial y^2} &\frac{\partial f}{\partial z^2} \
\frac{\partial f}{\partial r_1^2} &\frac{\partial f}{\partial r_2^2}&\frac{\partial f}{\partial r_3^2} &\frac{\partial f}{\partial t_1^2} &\frac{\partial f}{\partial t_2^2} &\frac{\partial f}{\partial t_3^2} &\frac{\partial f}{\partial x^2} &\frac{\partial f}{\partial y^2} &\frac{\partial f}{\partial z^2} \
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \
\frac{\partial f}{\partial r_1^n} &\frac{\partial f}{\partial r_2^n}&\frac{\partial f}{\partial r_3^n} &\frac{\partial f}{\partial t_1^n} &\frac{\partial f}{\partial t_2^n} &\frac{\partial f}{\partial t_3^n} &\frac{\partial f}{\partial x^2} &\frac{\partial f}{\partial y^2} &\frac{\partial f}{\partial z^2} \
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \
\frac{\partial f}{\partial r_1^n} &\frac{\partial f}{\partial r_2^n}&\frac{\partial f}{\partial r_3^n} &\frac{\partial f}{\partial t_1^n} &\frac{\partial f}{\partial t_2^n} &\frac{\partial f}{\partial t_3^n} &\frac{\partial f}{\partial x^m} &\frac{\partial f}{\partial y^m} &\frac{\partial f}{\partial z^m} \
\end{matrix}\right]
$$
三维重建一般都有成千上万个三维点和几十上百幅图像,所以雅克比矩阵非常庞大,每次迭代都需要计算一个如此大的矩阵,计算量也是非常巨大的,好在$J$ 是一个稀疏矩阵,也就是说有很多很多项都是0,比如当一个三维点在某一幅图像上没有投影点的时候,对应的那一行就为全为0了。