标签归档：算法

Levenberg-Marquardt算法

发表评论

匆匆一更。还没写完，对不起我太懒了T_T

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

Levenberg-Marquardt算法（下文简称LM算法）通常用于非线性最小二乘法的目标函数极小化。这是一个置信域方法（Trust-Region Method），为了防止步长太大而跳到非预期的局部极小值，这类算法自适应的调整步长。

假设f(x)为一个向量函数，每个分量都刻画了一个样本和模型预测之间的某种偏差，则非线性最小二乘目标函数为

(1)   \begin{eqnarray*} E&=&\transp{f}f \end{eqnarray*}

通常基于梯度的数值极小化算法，都是用一点的局域特征（如导数）建立该函数变化趋势的模型，以此来推测一个最可取的下降方向和步长。取一阶导数近似是线性模型，用切平面来近似函数在该点附近的变化趋势；取到二阶导数近似则是二阶模型，用一个抛物面来近似。当然，也可以取到更多阶，模型的几何直观也会越来越复杂，越来越精确的逼近真实的函数变化趋势，但是会付出更大的计算代价。当取到无穷阶时，模型在级数收敛区间内可以做到完全精确的刻画函数，此时这就是模型就是函数的泰勒展开了。

我们把f在任一点x附近展开f(x+\delta x)=f(x)+J(x)\delta x,这样E在x附近就是

(2)   \begin{eqnarray*} E(x+\delta x)&=&|f(x)+J(x)\delta x)|^2\\ &=&|f|^2+2\transp{f}J(x)\delta x + |J(x)\delta x|^2 \end{eqnarray*}

注意J是f而非E的导数，这里我们得到了E的二阶展开，即用一个抛物面来近似目标函数E在一点附近的形态。这个二阶模型显然是一个线性最小二乘问题——这正是我们所要的，把一个非线性最小二乘问题转换为一系列线性最小二乘问题来求解。

因此一个简单直接的想法就是，每一迭代步构造这样一个线性最小二乘问题，通过极小化它下降到下一个迭代点，循环直到收敛。但这个方案显然有问题。我们的二阶展开模型只能在一个小邻域上有效，当\delta x增大时，被我们忽略掉的高阶项就会越来越明显，最后让二阶模型彻底失效。这样，直接根据二阶模型计算出来的迭代点对目标函数而言可能并非我们想要的，它的函数值可能比上一点还要高（因而目标函数不降反增）；也可能尽管比上一点下降了但是却跳到了另一个我们不想要的凸区域里（因而算法会收敛到离初始猜测解更远的局部极小值点上）；更严重的情况，模型甚至可能是退化的，即一种极端形式的抛物面——就像一本被半卷起的书，此时我们在弯折的方向仍能找到极小值，但在其垂直方向已经无法找到极小值（或者说处处极小）。

为了解决以上问题，我们需要给模型加入一些约束，把迭代步长约束在一个合理的范围内，保证在这个范围内模型足够有效。同时，我们也希望在二阶模型退化的时候，约束能够拯救它，保证算法总能求出一个极小值。这个约束很简单，如下

(3)   \begin{eqnarray*} E(x+\delta x)&=&|f(x)+J(x)\delta x)|^2 + \lambda |C(x)\delta x|^2 \end{eqnarray*}

我们加入了一个二次约束项。单独看该项，它是一个中心在x点，开口向上的抛物面，或者说一个二次势阱。这个势阱像一个箍，把\delta x箍在x点附近的小邻域中。通过调节\lambda可以增强或减弱它对\delta x的约束：当\lambda充分大时，约束项成为主要贡献项，这时模型总有一个极小值，且步长较小；当\lambda较小时，约束项几乎可以忽略，算法以大步长快速往极小值接近。

再来看矩阵C(x). 这个矩阵通常被设置为对角的。此时，

(4)   \begin{eqnarray*} |C(x)\delta x|^2 &=& \transp{\delta x}\transp{C}C\delta x\\ &=&\sum_{i} \lambda_i \delta x_i^2 \end{eqnarray*}

\transp{C}C的本征值控制着每个方向上抛物面上升的快慢。可见，如果一个方向相应的本征值被设置的较小，那么该方向的约束项上升缓慢，约束较弱，步长跨度较大；否则反之。现实中使用的模型，目标函数对于各个参数的敏感度不同，有的参数略微调整就可以导致目标函数的强烈变化，有的则反之。C(x)这个参数矩阵使得我们可以针对各个不同的参数设置不同程度的步长约束，从而能够构造出一个稳定的算法。

（未完待补）

Notes on Conjugate Gradient Method

发表评论
版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

这是一个关于共轭梯度法的笔记。请大家注意的是，这是个笔记，并不是一个教程，因此少不了跳跃和欠解释的地方。对CG方法了解不多的同学请移步这里。

线性方程组和极小化问题

一个关于对称矩阵A的线性方程组Ax=b等价于求解如下极小值问题：

    \[ f(x) = \frac{1}{2} \transp{x} A x - \transp{b} x \]

这很容易说明，我们微分目标函数得

(1)   \begin{eqnarray*} \td f &=& \frac{1}{2}(\td\transp{x}Ax+\transp{x}A\td x)-\transp{b} \td x\\ &=& \transp{(Ax - b)} \td x \end{eqnarray*}

所以\td f=0意味着Ax=b.

设x^*为问题的解，e为偏离极小值点的位移，即

(2)   \begin{eqnarray*} Ax^*&=&b\\ e&=&\Delta x=x-x^* \end{eqnarray*}

我们定义残差或负梯度r=b-Ax，容易看出，r正比于偏离梯度零点x^*的位移e:

    \[ r=-Ae \]

最速下降法和共轭梯度法简述

最速下降法：搜索方向为本轮迭代初始点的梯度方向，搜索到梯度与搜索方向正交的位置开始下一轮迭代。

共轭梯度法：搜索方向为本轮迭代初始点的梯度方向（残差方向）用A-内积做正交化，扣除掉之前所有搜索方向的分量给出，搜索到梯度与搜索方向正交的位置开始下一轮迭代。

共轭梯度法的理解

首先考虑一个简单情形，即A=I时。此时f退化，且x的不同分量解耦，f的极小值问题分解为各分量上的极小值问题。我们只需要在各分量方向上找到极小值，问题就得解了。具体的说，我们从任意一个初始点出发，在方向d_0上寻找极小值，然后从这个极小值出发继续在d_1寻找d_1方向的极小，以此类推最后到达全局极小值。

A为正定对称矩阵的一般情形，其实只是上述情形中的x表述在一个新坐标系下的结果（下文详细解释）。为了求解这种情形，我们考察上述方法的求解轨迹在新坐标系里的对应。

首先，A=I时的各个正交搜索方向在新坐标系里对应着一组A-正交的方向，因此我们需要设法找出这样一组正交方向，然后各自求解这些方向上的极小值。

求解特定方向上极小值的方法不变，极小点就是这个方向上梯度正交于搜索方向时的点（原因很简单，这说明搜索方向上梯度分量为零）。

最后一点，如何找到这样一组A-正交方向？

一般方法是找一组现成的完备基，进行A-正交化。问题是计算量比较大。

记迭代到第i步时，所搜索过的这样一组A-正交的搜索方向单位矢量为\{d_0,d_1,...,d_{i-1}\},它们撑起的子空间

(3)   \begin{eqnarray*} D_i=\mathrm{span}\{d_0,d_1,...,d_{i-1}\} \end{eqnarray*}

即第i步前已经搜索过的子空间。

共轭梯度法的一个关键之处在于它可以从每次迭代的r_i中快速构造出新的搜索方向d_i，而不需要对r_i进行完整的正交化手续（即不需依次扣除r_i中d_0,d_1…d_{i-1}方向的分量）。原因在于r_i自然的和D_{i-1}A-正交，因此构造d_i只需要从r_{i}中扣除掉d_{i-1}分量即可。

因此，共轭梯度算法有两处关键，一处为通过利用A-正交关系来把搜索问题转换到一个特殊坐标系，使得各搜索方向上的极小值问题解耦；一是利用r_i和D_{i-1}是A-正交的这一规律来充分简化A-正交搜索方向的构造。

接下来我们详细论证这个算法。

第一个关键点，坐标变换和子问题解耦

设

(4)   \begin{eqnarray*} A&=&\transp{T}T\\ y&=&T x\\ b^{\prime}&=&T^{-T} b\\ \end{eqnarray*}

其中T^{-T}表示矩阵T求逆后做转置，很拗口，不过这个不重要。

于是我们有

(5)   \begin{eqnarray*} f(x) &=& \frac{1}{2} \transp{y} y - b^{\prime T} y\\ &=& \frac{1}{2}\sum_{n}y_n^2 - b^{\prime}_n y_n \equiv \frac{1}{2}\sum_{n} P_n(y_n) \end{eqnarray*}

可见f已经被分解为一系列独立的子问题P_n(y_n)，T正是将问题解耦所用的坐标变换。

在新的\{y_n\}坐标系中，整个极小值问题可以通过在一组正交方向上依次寻找极小点来解决。但求解T需要对A做分解，并非一个计算上经济的解决方法。这一个变换的真正意义在于，我们可以通过它来找到求解轨迹在\{x_n\}坐标系中的对应。

T是线性变换，因此\{y_n\}中的一组直线仍被变换为\{x_n\}中的一组直线。\{y_n\}坐标系中，我们依次沿着一组正交方向求解极小值，因此只需要找到这组方向在\{x_n\}坐标系中的对应（一组未必正交的方向），就可以等价的直接在这组方向之上寻找极小值，而不需要经过坐标变换。

考虑这组正交方向单位矢\{y_n\}

(6)   \begin{eqnarray*} &&\transp{y_i} y_j = \delta_{ij}\\ &\Leftrightarrow& \transp{(T x_i)} T x_j=\delta_{ij}\\ &\Leftrightarrow& \transp{x_i} \transp{T}T x_j = \delta_{ij}\\ &\Leftrightarrow& \transp{x_i} A x_j = \delta_{ij} \end{eqnarray*}

可见，\{y_n\}这一组正交的方向在变换到\{x_n\}坐标系中后，虽然不是\{x_n\}中通常意义的正交矢量组，却可以在A-内积\langle \bullet,\bullet \rangle_A意义下成为一组A-正交矢量。这里的A-内积，指的是以对称矩阵A定义的内积运算

(7)   \begin{eqnarray*} \langle a,b \rangle_A \equiv \transp{a} A b \end{eqnarray*}

需要注意的是，上述结论是可逆的，如果\langle x_i,x_j \rangle_A = 0,同样也可以得出，x_i,x_j对应着\{y_n\}坐标系中两个正交的方向。因此，可以直接在\{x_n\}坐标系下，寻找一组A-正交的方向，然后在这一组方向上依次极小化目标函数，就可以找到整个问题的极小值。

第二个关键点，D_i子空间结构和正交关系

算法的更新规则为

(8)   \begin{eqnarray*} e_{i+1}&=&e_i + \alpha_i d_i\\ r_{i+1}&=&-A(e_i+\alpha_i d_i)=r_i-\alpha_i A d_i\\ d_{i}&=&r_i + \sum_{j<i} \beta_{ij} d_{j}\\ \beta_{ij}&=& - \frac{\transp{r_i}d_j}{\transp{d_j}Ad_j}\\ &=& - \frac{\langle r_i,d_j \rangle }{\langle d_j,d_j \rangle _A} \end{eqnarray*}

其中\alpha_i为第i步搜索方向上的位移长度，它的大小可以通过r_{i+1}应该和d_i正交这一要求定出

(9)   \begin{eqnarray*} &&\transp{r_{i+1}}d_i=0\\ &\Rightarrow&\transp{(e_i+\alpha_i d_i)} A d_i =0\\ &\Rightarrow&\alpha_i = - \frac{\langle d_i,e_i \rangle_A }{\langle d_i,d_i \rangle_A}= \frac{\langle d_i,r_i \rangle}{\langle d_i,d_i \rangle_A} \end{eqnarray*}

算法将参数空间划分成了一个层级结构，下面我们来看看这个结构的性质。

根据上述更新规则，可以发现D_i可以由另外几个矢量组撑起。首先，由d_i的更新规则可以发现，d_i由r_i和d_{j<i}线性组合而成，另一方面，我们有初始条件d_0=r_0，归纳而知，d_i可以由\{r_{j\leq i}\}线性组合出来。于是我们有

(10)   \begin{eqnarray*} D_i=\mathrm{span}\{r_0,r_1,...,r_{i-1}\} \end{eqnarray*}

接着我们考虑r_i的更新规则。可以看到，D_{i+1}=\mathrm{span}\{D_{i},Ad_i\}.于是我们可以从D_1=\{d_0\}或D_1=\{r_0\}构造出D_i

(11)   \begin{eqnarray*} D_i&=&\mathrm{span}\{d_0,Ad_0,A^2d_0,...,A^{i-1}d_0\} = \mathrm{span}\{d_0,A D_{i-1}\}\\ &=&\mathrm{span}\{r_0,Ar_0,A^2r_0,...,A^{i-1}r_0\} = \mathrm{span}\{r_0,A D_{i-1}\} \end{eqnarray*}

接下来我们再考虑各子空间和向量之间的正交关系。首先，e_i在D_i的补空间里。这点是明显的，因为每一个优化步都会优化掉e_j在相应d_j方向上的分量，因而剩余的e_i就只位于D_i的补空间里。于是我们有

(12)   \begin{eqnarray*} e_i = \sum_{j \geq i} \delta_j d_j \end{eqnarray*}

这样e和d之间的A-正交关系就很明显了。由于\langle d_i,d_j \rangle_A = \delta_{ij}，而e_i只包含d_j>i分量，因而

(13)   \begin{eqnarray*} \langle e_i,d_{j<i} \rangle_A =0 \end{eqnarray*}

即e_i和子空间D_i也是A-正交的。这意味着r_i正交于D_i

(14)   \begin{eqnarray*} &&\transp{e_i} A d_{j<i}=0\\ &\Rightarrow&\transp{(A e_i)} d_{j<i}=0\\ &\Rightarrow&\langle r_i,d_{j<i}\rangle =0 \end{eqnarray*}

到这里就可以解释为什么我们有\langle r_i,d_{j<i-1} \rangle_A =0了：

(15)   \begin{eqnarray*} &&\langle r_i,D_i \rangle =0\\ &\Rightarrow&\langle r_i,\mathrm{span}\{d_0,A D_{i-1}\} \rangle =0\\ &\Rightarrow&\langle r_i,A D_{i-1} \rangle =0\\ &\Rightarrow&\langle r_i,D_{i-1} \rangle_A =0 \end{eqnarray*}

至此我们完成了第二个关键点的论证，且更清楚的了解了参数空间，尤其是D_i子空间的结构。