非线性最优化基础 简单易懂理解

一 非线性优化中的基本概念

(1). 方向导数：函数在自变量的某个方向上的导数。以一个二元函数为例，自变量有x和y，该函数在某一点方向有(x,y)这个向量组成的无穷多个。

(2).梯度：梯度是一个矢量，在这个方向上的方向导数达到最大值。某一点的梯度就是某一点方向导数的最大值，也就是从该点出发函数值变化最剧烈的方向（梯度方向代表增加最快的方向，负梯度方向代表减少最快的方向）。梯度的求法即为改点所有方向的偏导数组成的向量。

几乎所有的优化算法，都是找某个梯度方向，然后进行迭代。所以梯度方向对于最优化有重要意义。

(3). Jacobian矩阵:

Rn→Rm是一个从欧式n维空间转换到欧式m维空间的函数。（注意理解方式）可以理解为这个函数由m个实函数组成: y1(x1,…,xn), …, ym(x1,…,xn)，或者理解为这个函数的输出形式为这样的m维的向量。这些函数的偏导数(如果存在)可以组成一个m行n列的矩阵，这就是所谓的Jacobian矩阵，这个矩阵的第i行是由第i个函数的梯度函数的转置yi(i=1,…,m)表示的。

如果把[y1,y2,…,ym]看成向量，那么完全也可以理解为这是Y向量对X向量的求导。

 

特殊的，当我们只有一个函数时，或者说函数是一个N*1的形式时(自变量N维，输出1维)，Jacobian矩阵只是一个向量，

Jacobian矩阵相当于求一个多元函数在某点的导数：

更常用的是用它做线性逼近求最优化，我们用p点的值F(p)去逼近F(X)：

(4). Hessian矩阵

对于这样的一个函数： ，Hessian矩阵是它的二阶偏导数组成的方块矩阵。Hessian矩阵只针对一个函数，而不像Jacobian矩阵可以针对一个函数组。

Hessian矩阵的形式：

如果是Jacobian矩阵是一个多元函数的一阶导数，那么Hessian矩阵就相当于这个函数的二阶导数。

 

(5). 对矩阵求导

之前的固有意识里，我们认为求导只能是单个变量之间进行。现在我们把导数扩展到矩阵和向量。即单个标量，行列向量，矩阵之间都可以互相求导。

 

矩阵求导的规则很简单，就是逐个求导。理解成矩阵或向量的每个元素关于另外矩阵的每个量的变化。具体规则在：

http://files.cnblogs.com/files/leoleo/matrix_rules.pdf

 

 

 

二 优化思想与方法

1、泰勒公式的线性逼近

之前的blog已经写过泰勒公式，这次主要谈谈理解和灵活运用。我们知道，泰勒公式是为了求解一个复杂的函数在某处附近（x0附近）的近似函数。最常用的是泰勒公式的前两阶展开式。

从函数逼近的角度理解

（t11&t12）

利用上面的两个式子，我们知道，我们可以用x0（x0理解为一个离x很近的数）处的值，去逼近离x0不远的x处的值。也就是说，泰勒公式可以用一个离目标点近处的点，去求大概的目标点的值。

 

 

从迭代的角度理解，x可以写成x_k，而 x0可以写成x_k-1，意思就是第k次迭代f(x_k)的值，可以由不远处第k-1次迭代时候f(x_k-1)的值逼近而来

（t21&t22）

从增量的角度理解，我们可以把(x-x0)这个理解成增量，也就是说x是由x0变化了(x-x0)这么大的而产生的。因为迭代时候要用增量的性质，所以我们可以按照增量的概念去理解。我们已知了x的值，我们现在去求x增加了一个Δx后的值。

（t31&t32）

 

注意下面的算法中，全都用到了迭代的思想，也就是每一步我们是想办法找一个最优的Δx，曲线救国一步步去逼近最优的x，而不是直接找一个最优的x，这一点要明确注意。所以每一步的目标：找最优的Δx！

 

2、最速下降法(一阶梯度法)——用梯度直接当Δx

顾名思义，我们是每次找一个函数值下降最快的方向，不断更新x的值，最终下降到最优值。

我们给一个函数f(x)我们要求f(x)的最小值，我们需要打破初等数学学习的函数的概念，认识到f(x)可能是一个1à1的函数，也可能是Nà1的函数。

但对于解优化的问题，f(x)不可能是NàM的函数，原因在于如果输出是一个M维的向量，我们没办法定义输出的最值，所以必须输出是一个标量数字，我们这个问题才可以解。

我们已经讨论了，给定一点，该点下降最快的方向，就是反梯度的方向，而梯度的求法我们之前已经给出。所以，下一个逼近的x的值应该是：

 

其中α是步长，决定了每次迭代需要往反梯度方向走多远。当f是多元函数时，梯度用Jacobian矩阵来表示。所以一个通用梯度下降的算法步骤是

（1）    给一个能量函数f(X)，计算该函数的Jacobian矩阵表达式J(X)

（2）    估计一个迭代初始值X0，给一个步长α，给一个迭代终止条件t，即Δx<t时，认为收敛，停止迭代。

（3）    按照下面的式子更新 得到新的变量。

（4）    判断Δx<t，是否收敛。不收敛则回到(3)

（5）    输出取得最优解的自变量Xk 和最优值f(Xk)

 

（2）牛顿法（二阶梯度法）——解一个二阶泰勒展开的增量方程去找Δx

最速下降法中，我们只是选了一个当前下降最快的梯度方向作为迭代方向，因为我们知道梯度就是最快的方向。

现在我们用泰勒公式t32去推，有没有更好的梯度下降最快的方向。

通过t32我们知道，f(x)每次以这个值进行更新 ，现在我们想找 一个Δx，让f(x)更新的这个值最大。怎么做？

对这个更新量，求关于Δx的导数，令导数为0即可。

可以得

若为多元函数，用Hessian矩阵和Jacobian矩阵表示二阶和一阶导数，则有

所以更新公式变成了

所以通用牛顿法的步骤：

（1）    给一个能量函数f(X)，计算该函数的Jacobian矩阵表达式J(X)和Hessian矩阵表达式H(X)

（2）    估计一个迭代初始值X0，给一个步长α，给一个迭代终止条件t，即Δx<t时，认为收敛，停止迭代。

（3）    按照下面的式子更新 得到新的变量。

（4）    判断Δx<t，是否收敛。不收敛则回到(3)

（5）    输出取得最优解的自变量Xk 和最优值f(Xk)

 

注意：

（1）    由于考虑了二阶导数，牛顿法比梯度下降法有更快的收敛

（2）    牛顿法要求hessian矩阵必须要可逆而且正定，否则没法解。而且Hessian矩阵对于大规模问题计算起来比较复杂。

 

 

3、高斯牛顿法——解高斯牛顿方程去找Δx

高斯牛顿法是牛顿法的改进和替代，只用于解决最小二乘问题。它用其他的方式代替了牛顿法中增量方程的形式。

我们现在用一种简单的方式来定义最小二乘法拟合问题，首先设待拟合函数为h(θ), θ为待拟合参数，Xi是输入的数据，一组观测数据为Xi，Yi ，误差项fi(θ)代表第i组数据的残差，误差函数为误差项组成的向量f(θ)。  

(在我们的简单formulation里面，Xi为N维度，Yi为1维度，h(θ)是一个N*1的函数，如果不是这样N*1的函数，残差项就是一个无法定义。但实际上，Yi可以为任意M维度，但如果这时候的话，fi(θ)就必须是一个M*1的函数，这样才能保证每个误差项通过一个具体的标量定义出大小)。

简单定义一组(第i组)观测数据的残差(误差)项fi(X)，（实际上可以有别的定义方式，但必须保证误差项的输出是1维标量）：

 

所以所有的误差函数应该写成

由于向量二范数代表向量各项的平方和，所以f(θ)的二范数所组成的函数就是最终的优化目标。

最小二乘法优化目标可以写成

我们把θ换成x，这样更能从函数的角度解决这个问题，即f(θ)=f(X)。首先我们把误差函数（不是误差函数的平方和）一阶泰勒展开：

从迭代的角度，我们不是找最优的x，而是找最优的Δx。找的方法就是找一个Δx，让这个被一阶泰勒展开逼近的最小二乘法误差项最小。

我们对需要求最小值的函数进行展开，求它关于Δx的导数，并另导数为0求解Δx

展开得到：

（注意求导：如果A是标量,  A’=A；  d(X’MX)/dX=(d(X’MX)/dX)’=2*M’X，注意求导后转置，注意X’和X相乘代表二次的X）

求导，另导数为0，得到：

这个方程就是著名的高斯牛顿方程，又叫正规方程，它是一个线性方程组，我们现在来解线性方程组。J(x)是误差函数对待求参数 求导的Jacobian矩阵，f(x)是误差函数，注意误差函数目前是一个向量函数哦。

我们得到Δx的更新量为

然后按照下面的式子更新

只需要将牛顿法算法描述中的更新式换成上面的式子即可，并将学习速率或者步长α设为1即为Gassian-Newton算法。

 

算法评价：

该算法无需求复杂的Hessian矩阵，而把 作为其近似。缺点是这种近似矩阵也可能是奇异矩阵，而且步长也不好控制。但这个是解最小二乘最重要的方法，也是其他方法的基础核心。

 

4、 LM方法

 

 

对角矩阵D的含义：

拉格朗日数乘的含义：

在数学中的最优化问题中，拉格朗日乘数法（以数学家约瑟夫·拉格朗日命名）是一种寻找多元函数在其变量受到一个或多个条件的约束时的极值的方法。这种方法可以将一个有n个变量与k个约束条件的最优化问题转换为一个解有n+ k个变量的方程组的解的问题。这种方法中引入了一个或一组新的未知数，即拉格朗日乘数，又称拉格朗日乘子，或拉氏乘子，它们是在转换后的方程，即约束方程中作为梯度（gradient）的线性组合中各个向量的系数。

比如，要求 f(x,y)在g(x,y)=c时的最大值时，我们可以引入新变量拉格朗日乘数λ，这时我们只需要下列拉格朗日函数的极值：

拉格朗日乘数法所得的极点会包含原问题的所有极值点，但并不保证每个极值点都是原问题的极值点。拉格朗日乘数法的正确性的证明牵涉到偏微分，全微分或链法。

 

最终我们要解

对Δx进行求导，我们得到

 

  上式即为LM算法的增量方程，求得Δx作为更新量，其他的与guassian-newton步骤一样。

 

 

五、李代数变换与求导

1、基础知识

据我所知，我能用到李代数的地方也许只有计算机视觉中的相机变换了。所以我们只在旋转矩阵和仿射变换矩阵的角度讨论李群。

群是一种集合加上一种运算的代数结构，记作G=(A,·) 其中A是集合，·是运算。

李群是具有连续光环性质的群，离散群就不属于李群例如整数群。旋转矩阵这样的群是李群，因为他是刚体连续变化的运动。

 

引入一个概念：任意一个向量a，都可以找一个反对称矩阵a^来表示它，从而可以把叉积变成内积的计算方法。

a×b=a^·b

叉积是是表示两向量之间旋转的好方法，a和b是两个向量，那么从a转到b，可以以a×b的这个向量为对称轴，而|a×b|=|a||b|sinθ，所以旋转角是arcsin(|a×b|/|a||b|)。所以知道叉积后就可以描述一个绕轴的旋转变化，不过此变换仅在三维空间中成立。

同理，任意一个反对称矩阵，也可以找到一个向量来表示它。从而把点积变叉积。

 

2、李群SO(3)和SE(3)

简明给出直观的理解！！别扯没用的！

SO(3)就是旋转矩阵R构成的集合；

SE(3)就是变换矩阵 构成的集合；

李代数so(3)就是旋转矩阵R的向量表示，它是三维向量，向量长度代表旋转角，方向代表旋转轴；

李代数se(3)就是变换矩阵T的向量表示，它是六维向量，前三维是平移向量，后三维为旋转向量和so(3)一模一样；

 

 

 

 

我们略过推导，直接上关键结论

SO3

(1). 旋转矩阵是一个对称矩阵，旋转矩阵的导数经过推导，可以表示成一个由三维向量生成的反对称矩阵φ^乘以旋转矩阵本身的形式。（下式是旋转矩阵SO(3)关于时间的导数）

(2). Φ是SO(3)群对应的李代数，是三维向量（定义在R^3上），同样可以表示旋转。它与SO(3)李群的映射关系为：

这称为李代数的指数映射，可以从李代数映射到李群。

(3). 既然φ是一个旋转向量，那么我们可以表示成θa的形式，θ为模长，a是一个单位向量只表示方向。这样表示，代表一个以a为方向，以θ为角度的旋转。

我们有Rodriguez公式：

 

这给出了指数映射的显式计算方法，给出了李代数à李群的映射，给出了 给定 旋转轴和旋转角à计算旋转矩阵的方法。

 

(4). 从旋转矩阵到李代数（从旋转矩阵到旋转角和旋转轴），我们给出下面的变换来解得转轴a和转角θ：tr为矩阵的迹，代表主对角线元素之和

旋转轴上的向量经过旋转矩阵不变，解得a

 

这个叫对数映射，给出了李群à李代数的映射，给出了由旋转矩阵计算转角和转轴的方法。

 

 

 

 

SE3

(1).ξ是SE3对应的李代数，它是一个六维向量，ρ为平移的3维，φ为旋转的3维。

ξ也可以表示成一个矩阵，变换关系如下，但这里不再是反对称矩阵了（^符号只是扩展了向量到矩阵这个含义）

 

（2）指数映射（李代数à李群SE(3)，旋转+平移向量à变换矩阵）

需要将李代数写成ξ=[ρ,φ]=[ρ,θa]，这样拆解后，里面各元素即可映射到下面的式子里进行计算。

 

（3）对数映射（李群SE(3)à李代数，变换矩阵à平移+旋转向量）

 

 

总体来看，变换关系：

 

3. 李代数求导

(1) SO(3)求导

李代数最方便的求导方式就是利用导数的定义，进行扰动求导。

一般李代数李群（变换矩阵）作为自变量，其他因素作为因变量

求导形式：求空间位置对于旋转矩阵的导数

扰动方式：给空间点乘以一个旋转矩阵作为扰动

目标函数：函数是关于R旋转矩à空间位置的函数

扰动方式： 根据前面的李代数变化，意义是扰动矩阵，乘以空间点位置就相当于给了目标函数一个小变化。

无需计算Jacobian矩阵，非常方便。以后空间位置对旋转矩阵求导时，直接把R·p运算出的向量作反对称变换即可求出，无需多余操作。

 

 

(2) SE(3)求导

给变换矩阵T左乘一个扰动 ，设扰动项李代数 ，其中T的形式为 有：

 

最后一个运算符 代表把T·p一个其次坐标的空间点变换，变成一个4*6的矩阵。所以当空间位置对变换矩阵T求导时，直接按照这个变换写就行了。

 

 

 

参考高翔《视觉SLAM十四讲》 电子工业出版社