非线性最小二乘法
一.梯度下降法以及Jacobian矩阵计算
在2010年的关于L-K和AAM的博客里提到，模板匹配公式的一阶泰勒展开ΔT=J*Δp，J是用于梯度下降的Jacobian矩阵，是高维矢量函数值T=f(p)相对与参数矢量p变化时的增量(导数)。如果p是n维矢量，T是M维矢量，则J是一个[m*n]的矩阵。J在(i,j)处的元素值是 (əTi/əpj)。
对模板匹配问题，假设矢量ΔT是图像Patch和模板的差，Δp是模板匹配每次迭代的参数增量，模板匹配过程就是根据ΔT求Δp的过程：
Δp =(JTJ)-1JTΔT
其中JTJ是高斯牛顿法中的海塞矩阵H = JTJ。因为H并不总是可逆，所以
Δp =(H+λI)-1JTΔT
关于λ的说明：(无约束)梯度下降法求函数最小化问题时，λ的增加，减弱了海塞矩阵H的影响，λ从0到正方向的递增过程对应了梯度下降法从牛顿法到最速下降法的切换。优化过程中根据目标函数的结果实时调整λ值的方法即是所谓Levenberg-Marquardt法。 λ的另外一个作用是可以限制回归结果ΔP的幅值，称为Ridge regression。
所以求解模板匹配方法要在每步迭代时计算当前位置的梯度矩阵J，也就是估计输入图像（或模版图像）在参数变化时对应的变化量。J的计算可以使用解析法和数值法：
1.解析法：如果模板的几何变形可以用显式的公式来表达，比如Lucas-kanade（L-K）方法中图像Patch的变形可以用2个参数的位移，6个参数的仿射变换或8个参数透射变换Homography来表达。əTi/əpj可以用链式法则表示成图像的梯度和几何变化的内积的形式əTi/əpj=<əTi/əW,əW/əp>;再如AAM revisited文章介绍的face alignment方法，二维人脸特征点经过三角剖分形成多个三角面片，人脸这个非刚体变形模板可以用对应三角形的仿射变换来表达。
2.数值法：如果是几何参数和模板的几何对应关系是非线性的,比如机器视觉的成像投影函数是非线性的，再如图像模板可能不是简单地用图像的灰度矢量来表示的，而是用高级的图像冗余特征如梯度、Gabor小波或者SIFT/HOG等表达，这些高维的表达方式非常不适合计算梯度。数值法求J是在对当前每个参数pj增加一个小的扰动量(Purterbation)Δpj，计算因此产生的增量ΔTi和Δpj的比值：ΔTi/Δpj。因为P的每个分量是单独计算的，假设pj彼此之间是相互独立的，则矩阵J的列之间的相关性应该比较小。
另外，仅仅是为了提高J的鲁棒性，也可以用数值法计算，比如在上述解析法光流中求J需要计算图像或模板的梯度，其中就有很多变化，常规的Sobel算子只能考虑3×3邻域，必然受图像噪声的影响，此时可以考虑使用多尺寸模板度卷积的均值。再如Tim Cootes 2001年最初发表的AAM方法，计算人脸图像相对于人脸形状参数变化的J时，就是使用多次扰动结果的均值进行计算。
可以使用梯度下降法求解的非线性最小二乘问题：
(1)L-K方法image registration和AAM人脸对齐：根据输入图像和模板(或子空间模型)的差，推出形状参数的变化。为了减少计算量，Jacobian矩阵J和海塞矩阵H不是在图像I上而是在模板T(对于AAM是人脸训练样本的平均Appearence A0处)上进行计算的。第k次迭代时，根据参数Pk-1把输入图像warp到模板T附近（对于AAM是利用刚性参数和PCA形状参数到输入图像上采样，填充模板图像位置的像素值），然后利用warp后的图像和模板的差计算得到Δpk 并更新形状参数Pk = Pk-1+Δpk（有时是矩阵乘法）。(就是所谓Inverse Composition？）。文章Lucas-kanade20 years on: A unifying framework就是最初把L-K方法和AAM统一到模板匹配的框架下的。
(2)根据已知3D物体的空间结构和对应的一个2D图像投影估计3D物体的姿态（OpenCV的solvePnP()）：Pose矩阵[R|T]有6个自由度，该问题是有6个参数的梯度下降问题。误差函数优化的目标是使估计的投影点和实际的投影点之间的距离最小。一个典型的应用是根据Face Alignment的2D特征点估计3D人脸模型的姿态。
(3)机器视觉的光束平差(Bundle Adjustment)：根据3D空间物体（通常是根据图像特征点匹配恢复出来的稀疏点云）的多张2D图像上的投影点P2d，同时计算每个摄像机的6自由度的Pose [R|T]、所有3D空间点的坐标P3d=(x,y,z)T，以及每个图像的摄像机内参K（以及非线性的畸变参数）。2D图像i上投影的匹配误差是3D重建的投影点和2D观察点之间的距离：D isti=|P2d-K*[R|T]*P3d|,所有图像上的匹配误差最小化minE = min sum(Disti)。Bundle Adjustment在张正友的相机标定方法中用于估计相机内参和畸变参数，标定棋盘上的点的空间位置和相机姿态作为副产品，也可以同时得到。通过已知内参的单目相机实时估计相机姿态[R|T]和三维点云坐标的算法叫做SfM或SLAM。如果摄像机是变焦的，原理上也能实时计算出相机内参，称为相机自标定。
(4)两个RGBD图像（带颜色的深度图）配准的变换矩阵[R|T]（ICP算法）：用于计算RGBD Odometry。可参考[4]。刚体配准需要6个参数，把一个深度图像形成的点云投影到另外一个深度图像(可视为模板)上进行配准。匹配的误差是两个RGBD的RGB及深度Z的差（也可以是3D点云坐标XYZ的距离）。另外，也有利用深度图像的TSDF体数据进行等值面配准的方法。
Matlab里有求解非线性最小二乘问题的函数X=lsqnomlin(myfun,X0,LB,UB,options)，特别好用。只要定义好表达矢量差的函数myfun()，以及给定初值X0，设置options等参数（比如设置是否使用L-M算法等），程序会自动计算Jacobian求出最优解X。以CameraCalibration为例，myfun()中只要定义好所有3D-2D投影和检测到的2D点集的误差矢量，并且把所有待定参数，包括n个相机的K和pose、m个3D点坐标的初始值连接成参数矢量X0，送进优化函数，返回的优化结果X再解码成对应参数。实际上初始参数选择比较随意，甚至不需要太精确的线性解(张正友的文章中的3.1 Close-formsolution一节)步骤就能收敛到精确解。也有许多基于C++的开源库能很好地求解NLS问题，比如Eigen库中就有函数能够完成类似lsqnomlin()的功能。还有专门求解BA问题的开源库。
二. 通过学习的方法计算下降矩阵D
其实就是待定系数法，或者看成拟合超平面的回归问题。直接求解下面式子中的D
Δp = D *ΔT
在梯度下降过程的当前位置通过给p加一组扰动采样得到对应的ΔT样本。假设样本个数为m，参数P是np维，T是nt维的。则上式的矩阵维度是
[np*m]=[np*nt][nt*m]
Δp和Δt的采样方法和上面提到的估计J的数值法类似，不同的是：样本采样除了希望相互独立外(高维空间中的采样：蒙特卡洛法？)。m的个数也可以远远大于参数p的维度np。D的计算：
D = Δp*ΔTT*(ΔT*ΔTT)-1
其中(ΔT*ΔTT)是[nt*nt]维的矩阵，如果图像patch T的每个像素是由HOG或SIFT等特征表达的，比如图像Patch的size是32*32，每个pixel位置的HOG特征是128维的，T就是nt=32*32*128=131072维矢量，这样实时计算(ΔT*ΔTT)的逆矩阵时的计算量很大。解决的方法：
(1).SVD法或PCA降维，比如在CVPR2013 SDM文章中使用PCA对样本进行降维，最终的特征矢量T是乘以一个降维矩阵的结果。
(2).梯度下降法迭代求解，D(t) = D(t-1) + η*(-əE/əD)，E是拟合残差，η是学习率。这是神经网络用于更新连接权重的常用方法，它的好处在于可以支持对D矩阵的on-line学习，即每次只使用少量(mini-batch)甚至单个的训练样本对D进行更新。相当于一个使用线性激活函数且只有一层的神经网，适合于内存受限的场合。
(3).参考“Online Learning of LinearPredictors for Real-Time Tracking”，把对[nt*nt]维矩阵求逆的过程转化为对[np*np]矩阵求逆的过程(通常nt>>np)。具体见下图(其中模板Homography变形np=8)：
