当一个优化问题其中的Cost function 或是约束条件包含非线性方程式时，使得优化问题的解决方案变得相当棘手，这种问题称为非线性优化Nonlinear Optimization问题.

Nonlinear Optimization 问题简介

如果我们想要最小化函数 ，通常使用微积分中的方法：找到一组 ，使得对于所有 ，函数 对 的偏导数都等于零，即 。在优化问题中，一般假设函数 是平滑且可微的，因此可以通过泰勒展开式来解决优化问题：

这里的 表示梯度（gradient），而 则表示 Hessian 矩阵：

梯度下降法

梯度下降法的本质是沿着梯度最陡的方向进行移动，以达到函数的极小值。梯度下降的公式为：

这里 是梯度， 是每次移动的步长，也称为学习率。

牛顿法

如前所述，希望找到一组 ，使得 。牛顿法的思想是在当前点 和当前梯度 的基础上，找到 ，使得 。根据泰勒展开式，可以得到以下近似关系：

这里的 是 Hessian 矩阵，可以看作是梯度函数 的导数。根据上式，可以推导出：

因此，牛顿法利用以下式子不断更新 的值：

梯度下降法与牛顿法的比较

牛顿法收敛速度较快。梯度下降法将 乘以 ，而牛顿法将 乘以 。显然两者都更新 的方向不同。梯度下降法是沿着 更新 的值，而牛顿法是沿着 的方向更新。从直觉上来说，梯度下降法是一阶收敛，仅考虑当前位置，并寻找最陡的方向迈出一步；而牛顿法是二阶收敛，同时考虑迈出一步后梯度是否会变得更大。因此牛顿法的收敛速度通常比梯度下降法更快。

高斯牛顿法

如果牛顿法中要计算 Hessian matrix 的时间与空间的复杂度太大，而高斯牛顿法的实质就是去近似 Hessian matrix 进而降低梯度。高斯牛顿法的前提是这个最佳化问题必须为least square problem ,也就是以下式子：

以上的问题当然可以用梯度下降法或牛顿法来解，但是如果用高斯牛顿法的话会更有效率。

基本描述

给定 个函数 （通常称为残差）关于 个变量 满足 高斯-牛顿算法迭代地找到使得平方和最小化的 的取值

从初始猜测 开始，该方法通过迭代进行如下计算

符号 表示矩阵的转置。其中，如果 和 是列向量，则雅可比矩阵的元素为

在每次迭代中，可以通过重新排列前述方程来确定更新 的值：

通过将 , , 和 代入上述表达式，它可以转化为形式为 的常规矩阵方程。

特别的，如果 ，则迭代简化为

在数据拟合中，目标是找到参数 ，使得给定的模型函数 最佳拟合一些数据点 ，其中函数 是残差：

然后，高斯-牛顿算法可以用函数 的雅克比矩阵 来表示：

需要注意 是 的左伪逆。

具体解法

通常情况下，否则矩阵 不可逆，无法解出正规方程。

高斯-牛顿算法可以通过对向量函数 进行线性近似得到。使用泰勒定理，可以在每次迭代时写出：

其中 。找到使右边平方和最小的 ，即

问题转化为是一个线性最小二乘问题，可以显式求解，方法是给出正规方程。

正规方程是 个未知增量 的联立线性方程，可以通过乔列斯基分解（Cholesky decomposition）或更好地通过的QR分解来一步解决。对于大型系统，迭代方法，例如共轭梯度法可能更有效。如果的列之间存在线性依赖关系，迭代将失败，因为成为奇异的。

PS：当是复数，还可以使用共轭形式：。

具体案例

使用高斯-牛顿算法，通过最小化数据与模型预测之间的误差平方和来拟合模型。

在一项关于底物浓度和酶介导反应的反应速率之间关系的生物学实验中，通过测试得到了以下表中的数据。

希望找到最符合最小二乘意义下数据的曲线（模型函数）形式如下：

其中和是要确定的参数。

用和分别表示和rate的值，。令，，我们的目标是找到和，使得残差平方和

最小化。

残差向量的雅可比矩阵关于未知数是一个矩阵，其第行计算为

从和的初步估计开始，进行了高斯-牛顿算法的五次迭代后，得到了最优值和。残差平方和从初始值1.445减少到第五次迭代后的0.00784。右侧图中的绘图显示了由模型确定的最优参数的曲线与观测数据。

收敛性质

高斯-牛顿迭代在以下四个条件下保证收敛至局部最小点：函数在包含的开凸集上两次连续可微，雅可比矩阵具有满列秩，初始迭代点靠近，局部最小值很小。如果，则收敛是二次的。

可以证明增量Δ是的下降方向，并且如果算法收敛，则极限是的稳定点。然而，对于较大的最小值，收敛不被保证，甚至不能保证局部收敛，就像牛顿法那样，或者不满足通常的Wolfe条件。

高斯-牛顿算法的收敛速度可以达到二次。如果初始猜测远离最小值，或者矩阵的条件不良，则该算法可能收敛缓慢或根本不收敛。例如，考虑具有个方程和个变量的问题，如下：

最优解是。（实际上，对于，最优解是，因为，但）。如果，那么问题实际上是线性的，并且该方法可以在一次迭代中找到最优解。如果，那么该方法会线性收敛，误差在每次迭代中渐近地减小。然而，如果，那么该方法甚至不会局部收敛。

求解超定方程组

高斯-牛顿迭代

是解形式为的超定方程系统的有效方法，其中

且，这里是雅可比矩阵的Moore-Penrose逆（也称为伪逆）。它可以被视为牛顿法的扩展，在孤立的正则解处也有相同的局部二次收敛。

如果解决方案不存在但初始迭代接近于点，在该点上平方和达到一个小的局部最小值，那么高斯-牛顿迭代会线性收敛至。点通常被称为过度确定系统的最小二乘解。

从牛顿法推导得到高斯牛顿法

接下来，高斯-牛顿算法将从函数优化的牛顿法中通过一种逼近被推导。因此，在某些正则条件下，高斯-牛顿算法的收敛速度可以是二次的。一般来说（在较弱的条件下），收敛速度是线性的。

用于最小化参数的牛顿法的递推关系是

其中g表示S的梯度向量，H表示S的Hessian矩阵。

由于，梯度由以下给出

Hessian的元素是通过对梯度元素与的微分来计算的：

高斯-牛顿法是通过忽略二阶导数项（这个表达式中的第二项）得到的。也就是说，Hessian被近似为

其中是Jacobian矩阵Jr的元素。请注意，当在精确拟合附近评估精确Hessian时，我们有接近零的，因此第二项也变得接近零，从而证明了这个近似。梯度和近似Hessian可以通过矩阵表示为

以下是英文学术论文中的一节，将其翻译成中文:

这些表达式被代入上述递推关系中以获得操作方程

高斯-牛顿法的收敛并非在所有情况下都能得到保证。近似式

需要满足两种情况才能忽略二阶导数项，从而预期收敛：

• 函数值在最小值周围的幅度较小。
• 函数仅“轻微”非线性，因此幅度相对较小。

改进高斯-牛顿法

使用高斯-牛顿法时，残差平方和S可能不会在每次迭代中减小。然而，由于Δ是一种下降方向，除非是一个稳定点，否则成立 对于所有足够小的。因此，如果发散，解决方案之一是在更新公式中采用增量向量Δ的分数：

换句话说，增量向量太长，但仍指向“下坡”，所以沿部分路程前进将减小目标函数S。通过使用线性搜索算法可以找到的最佳值，即通过在区间中找到最小化S的值来确定的大小，通常使用直接搜索方法或Armijo线搜索等回溯线搜索。通常，应被选择为满足Wolfe条件或Goldstein条件的值。

在增量向量的方向使得最佳分数α接近于零的情况下，处理发散的另一种方法是使用Levenberg-Marquardt算法，即信赖域方法。正规方程式以一种使得增量向量朝向最陡降方向旋转的方式进行修改，

其中D是正对角矩阵。

请注意，当D是单位矩阵I且时，，因此Δ的方向逼近于负梯度的方向。

所谓的Marquardt参数 也可以通过线性搜索来优化，但这种方法效率低，因为每次改变 时都必须重新计算移位向量。更高效的策略是：当发散发生时，增加Marquardt参数，直到S减小为止。然后保留从一个迭代到下一个的值，如果可能的话，逐渐减小，直到达到截断值为止，然后可以将Marquardt参数设为零；这时S的最小化就变成了标准的Gauss–Newton最小化。

应用于大规模优化

对于大规模优化，Gauss–Newton方法有其应用前景，因为通常（尽管并非总是如此），矩阵 比近似Hessian 更为稀疏。在这种情况下，步长计算本身通常需要用适用于大规模稀疏问题的近似迭代方法来完成，例如共轭梯度法。

为了使这种方法有效，至少需要一种高效的方法来计算产品

对于某个向量，通过稀疏矩阵存储，通常可以以压缩格式（例如，无零条目）存储 的行，使得由于转置而直接计算上述产品复杂。然而，如果将定义为矩阵 的第i行，则有以下简单关系：

因此，每一行都会相加并独立地对产品做出贡献。除了保持实际的稀疏存储结构外，这个表达式也很适合并行计算。注意，每一行都是相应残差 的梯度；考虑到这一点，上面的公式强调了残差独立地对问题产生影响的事实。

MATLAB实现

function beta = gaussnewton(r, J, beta0, maxiter, tol)
    beta = beta0;
    for i = 1:maxiter
        Jbeta = J(beta);
        delta = -pinv(Jbeta' * Jbeta) * Jbeta' * r(beta);
        beta = beta + delta;
        if norm(delta) < tol
            break;
        end
    end
end
function beta = gaussnewton(r, beta0, maxiter, tol)
    beta = beta0;
    for i = 1:maxiter
        [rbeta, Jbeta] = jacobian(r(beta));
        delta = -pinv(Jbeta' * Jbeta) * Jbeta' * rbeta;
        beta = beta + delta;
        if norm(delta) < tol
            break;
        end
    end
end