Non-Linear Programming

Posted on Edited on In Waline:

本文结合网络资源对非线性规划算法进行整理。

非线性规划

1 最速下降法

1.1 解决的问题

最速梯度下降法解决的问题是无约束优化问题，而所谓的无约束优化问题就是对目标函数的求解，没有任何的约束限制的优化问题，比如求 min⁡f(x),f:Rn→R\min f(x), f:R^n\rightarrow R。

求解这类的问题可以分为两大类：最优条件法和迭代法。

  • 最优条件法是是指当函数存在解析形式，能够通过最优性条件求解出显式最优解。对于无约束最优化问题，如果f(x)在最优点x附近可微，那么x是局部极小点的必要条件为：df(x∗)df(x^*) 我们常常就是通过这个必要条件去求取可能的极小值点，再验证这些点是否真的是极小值点。当上式方程可以求解的时候，无约束最优化问题基本就解决了。
  • 实际中，这个方程往往难以求解。这就引出了第二大类方法：迭代法。

1.2 最速梯度下降法

1.2.1 算法流程

  1. 选取初始点 x0x^0，给定终止误差 ε\varepsilon，令 k=0k=0；
  2. 计算 ∇f(xk)\nabla f(x^k)，若 ∣∣∇f(xk)∣∣<ε||\nabla f(x^k)||<\varepsilon，终止迭代，输出 x∗≈xkx^*\approx x^k，否则，进入 3；
  3. 取 pk=−∇f(xk)p^k=-\nabla f(x^k)；
  4. 进行一维搜索，求 tkt_k 使得 f(xk+tkpk)=min⁡t≥0f(xk+tpk)f(x^k+t_kp^k)=\underset{t\ge 0}{\min}f(x^k+tp^k)，令 xk+1=xk+tkpk,k=k+1x^{k+1}=x^k+t_kp^k,k=k+1，转 2。

1.2.2 步长 tkt_k 的确定

1.2.2.1 方法一：一维寻优法

上述第 4 步中，f(xk+tpk)f(x^k+tp^k) 已成为步长 tt 的一元函数，故可以用任何的一维寻优法求出 tkt_k，即

f(xk+1)=f(xk−tk∇f(xk))=min⁡t∈Rf(xk−t∇f(xk))f(x^{k+1})=f(x^k-t_k\nabla f(x^k))=\underset{t\in \mathbb{R}}{\min} f(x^k-t\nabla f(x^k))

1.2.2.2 方法二：微分法

令

φ(t)=f(xk+tpk)\varphi(t)=f(x^k+tp^k)

则在最优点 tkt_k 有

φ′(tk)=∇f(xk+tkpk)Tpk=−(pk+1)Tpk=0\varphi'(t_k)=\nabla f(x^k+t_kp^k)^Tp^k=-(p^{k+1})^Tp^k=0

这也证明了最速下降法在两个相邻点之间的搜索方向是正交的。

1.3 最速下降法与梯度下降法的区别

  • 梯度下降法的步长 tt 是给定的

xk+1=xk−t∇f(xk)x^{k+1}=x^k-t\nabla f(x^k)

  • 最速下降法的步长通过优化函数搜索得到

tk=argmin⁡⁡t∈Rf(xk−tk∇f(xk))t_k=\underset{t\in\mathbb{R}}{\operatorname{arg\min}}f(x^k-t_k\nabla f(x^k))

1.4 最速下降法的缺点

在一般情况下，当用最速下降法寻找极小点时，其搜索路径呈直角锯齿状，在开始几步，目标函数下降较快；但在接近极小点时，收敛速度就不太理想。特别是当目标函数的等值线为比较扁平的椭圆时，收敛就更慢了。

2 高斯牛顿法

2.1 牛顿法

2.1.1 用牛顿法求解零点

用牛顿法求解方程 f(x)=0f(x)=0 时，将 f(x)f(x) 在 xkx_k 做一阶展开

f(xk)+f′(xk)Δxk=0f(x_k)+f'(x_k)\Delta x_k=0

第 k→k+1k\rightarrow k+1 次迭代中，变化量为

Δxk=−H(xk)−1f(xk)\Delta x_k=-H(x_k)^{-1}f(x_k)

其中 H(xk)≈f′(xk)H(x_k)\approx f'(x_k)。

2.1.2 用牛顿法求解无约束优化

2.1.2.1 理解一

对于无约束优化问题 min⁡x∈RnF(x)\underset{x\in \mathbb{R}^n}{\min}F(x)，转为求解器一阶最优化的必要条件

F′(x∗)=0F'(x^*)=0

便将无约束优化求解“转化”为了方程零点的求解。

由上面部分做的一阶展开，可以得到

H(xk)≈F′′(xk)H(x_k)\approx F''(x_k)

牛顿法求解优化问题时，实际求解的是

F′(xk)+F′′(xk)Δxk=0F'(x_k)+F''(x_k)\Delta x_k=0

当上式和 H(xk)≻0H(x_k)\succ 0 同时满足时，原目标函数在 xkx_k 的二阶展开的局部最小条件也被满足，即

Δxk=argmin⁡⁡ΔxF(xk)+F′(xk)Δx+12ΔxTF′′(xk)Δx\Delta x_k=\underset{\Delta x}{\operatorname{arg\min}}F(x_k)+F'(x_k)\Delta x+\frac{1}{2}\Delta x^TF''(x_k)\Delta x

也就是说，牛顿法单步迭代求解的过程可以理解为对原函数在 xkx_k 的二阶展开求最小值。

每一次迭代过程取增量为

Δxk=−H(xk)F′(xk)=−[∇2F(xk)]−1∇F(xk)\Delta x_k=-H(x_k)F'(x_k)=-[\nabla^2 F(x_k)]^{-1}\nabla F(x_k)

2.1.2.2 理解二

我们对目标函数 F(x)F(x) 在 xkx_k 进行二阶泰勒展开

F(xk+Δxk)≈F(xk)+J(xk)TΔxk+12ΔxkTH(xk)ΔxkF(x_k+\Delta x_k)\approx F(x_k)+J(x_k)^T\Delta x_k+\frac{1}{2}\Delta x_k^TH(x_k)\Delta x_k

式中仅 Δxk\Delta x_k 为变量，上式为一个以 Δxk\Delta x_k 为自变量的二次函数，且二次项系数为正（Hesse 矩阵正定 H(xk)≻0H(x_k)\succ 0）的情况下可以取得最小值，在最小值点满足导数为 0，即

dF(xk+Δxk)dΔxk=0\frac{\mathrm{d}F(x_k+\Delta x_k)}{\mathrm{d}\Delta x_k}=0

解得增量方程

H(xk)Δxk=−J(Δxk)H(x_k)\Delta x_k=-J(\Delta x_k)

第 k+1k+1 轮迭代取

xk+1=xk+Δxkx_{k+1}=x^k+\Delta x_k

2.1.2.3 算法流程

  1. 给定初始值 x0x_0 和精度阈值 ε\varepsilon，并令 k=0k=0；
  2. 计算 H(xk)H(x_k) 和 J(xk)J(x_k)；
  3. 若 ∣∣J(xk)∣∣<ε||J(x_k)||<\varepsilon，终止迭代输出 x∗≈xkx^*\approx x_k，否则确定抖索方向为 Δxk=−H(xk)−1J(xk)\Delta x_k=-H(x_k)^{-1}J(x_k)；
  4. 计算新的迭代点 xk+1=xk+Δxkx_{k+1}=x_k+\Delta x_k，并令 k=k+1k=k+1，转 2。

2.1.3 阻尼牛顿法

2.1.3.1 与牛顿法的区别

相比于牛顿法，阻尼牛顿法会对每一轮迭代的进行最优步长搜索，步长 λk\lambda_k 为

λk=argmin⁡⁡λ∈RF(xk+λΔxk)\lambda_k=\underset{\lambda\in\mathbb{R}}{\operatorname{arg\min}}F(x^k+\lambda \Delta x_k)

此时，第 k+1k+1 轮迭代取

xk+1=xk+λkΔxkx_{k+1}=x^k+\lambda_k\Delta x_k

2.1.3.2 算法流程

  1. 给定初始值 x0x_0 和精度阈值 ε\varepsilon，并令 k=0k=0；
  2. 计算 H(xk)H(x_k) 和 J(xk)J(x_k)；
  3. 若 ∣∣J(xk)∣∣<ε||J(x_k)||<\varepsilon，终止迭代输出 x∗≈xkx^*\approx x_k，否则确定抖索方向为 Δxk=−H(xk)−1J(xk)\Delta x_k=-H(x_k)^{-1}J(x_k)；
  4. 利用 λk=argmin⁡⁡λ∈RF(xk+λΔxk)\lambda_k=\underset{\lambda\in\mathbb{R}}{\operatorname{arg\min}}F(x^k+\lambda \Delta x_k) 确定搜索步长；
  5. 计算新的迭代点 xk+1=xk+λkΔxkx_{k+1}=x_k+\lambda _k\Delta x_k，并令 k=k+1k=k+1，转 2。

2.2 高斯牛顿法

高斯牛顿法针对最小二乘问题，采用一定的方法对牛顿法中的 Hesse 矩阵 H(xk)H(x_k) 进行近似，从而简化了计算量。

考虑最小二乘法问题

min⁡x∈RnF(x)=12∣∣f(x)∣∣2\underset{x\in\mathbb{R}^n}{\min}F(x)=\frac{1}{2}||f(x)||^2

当 xx 在 xkx_k 处，具有增量 Δxk\Delta x_k，优化目标函数取值为：

min⁡x∈RnF(xk+Δxk)=12∣∣f(xk+Δxk)∣∣2\underset{x\in\mathbb{R}^n}{\min}F(x_k+\Delta x_k)=\frac{1}{2}||f(x_k+\Delta x_k)||^2

不同于牛顿法，我们不对优化目标函数进行泰勒展开，我们对优化目标函数中的一部分，即 f(x)f(x) 进行一阶泰勒展开。展开后优化目标函数为：

F(xk+Δxk)=12∣∣f(xk)+J(xk)TΔxk∣∣2F(x_k+\Delta x_k)=\frac{1}{2}||f(x_k)+J(x_k)^T\Delta x_k||^2

注意，不同于牛顿法，高斯牛顿法的雅可比矩阵 J(xk)J(x_k) （列向量）是函数 f(x)f(x) 的，而不是目标函数。

式中仅 Δxk\Delta x_k 为变量，上式为一个以 Δxk\Delta x_k 为自变量的二次函数，且二次项系数为正（Hesse 矩阵正定 H(xk)≻0H(x_k)\succ 0）的情况下可以取得最小值，在最小值点满足导数为 0，即

d12∣∣f(xk)+J(xk)TΔxk∣∣2dΔxk=0\frac{\mathrm{d}\frac{1}{2}||f(x_k)+J(x_k)^T\Delta x_k||^2}{\mathrm{d}\Delta x_k}=0

解得

J(xk)J(xk)TΔxk=−J(xk)f(xk)J(x_k)J(x_k)^T\Delta x_k=-J(x_k)f(x_k)

记

H(xk)=J(xk)J(xk)T, g(xk)=−J(xk)f(xk)H(x_k)=J(x_k)J(x_k)^T,\ g(x_k)=-J(x_k)f(x_k)

因此得到高斯牛顿法的增量方程（搜索方向）

H(xk)Δxk=g(xk)H(x_k)\Delta x_k=g(x_k)

这里 H(xk)H(x_k) 是对 Hesse 矩阵的一种近似。

2.3 高斯牛顿法计算流程

  1. 给定初始值 x0x_0，令 k=0k=0；
  2. 对于第 kk 次迭代，求出当前雅可比矩阵 J(xk)J(x_k) 和 f(xk)f(x_k) ；
  3. 求解增量方程 H(xk)Δxk=g(xk)H(x_k)\Delta x_k=g(x_k)；
  4. 若 Δxk\Delta x_k 足够小，则停止，输出 x∗≈xkx^*\approx x^k。否则，令 xk+1=xk+Δxkx_{k+1}=x_k+\Delta x_k，返回 2。

2.4 高斯牛顿法的缺点

  • 求解 Δxk\Delta x_k 时需要求解 (J(xk)J(xk)T)−1(J(x_k)J(x_k)^T)^{-1}，但是 J(xk)J(xk)TJ(x_k)J(x_k)^T 只有半正定性，极有可能是奇异矩阵或病态矩阵，导致增量不稳定，算法无法收敛。
    • 奇异矩阵：不满秩矩阵。
    • 病态矩阵：求解方程组时，对数据较小的扰动就会导致结果较大的波动。
  • 若求得的 Δxk\Delta x_k 过大就会使得泰勒展开的局部近似不准确，导致不收敛。

3 L-M（列文伯格-马夸尔特）法

3.1 L-M 法思想

针对高斯牛顿法的不足，L-M 法做了两点改进：

  • 当求解增量 Δxk\Delta x_k 时，对其设置置信区间；
  • 在求解增量 Δxk\Delta x_k 时，对其近似效果进行了量化，并根据量化结果对信赖区域进行调整，再重新计算增量Δxk\Delta x_k，直到近似效果量化结果达到阈值。

3.2 置信区间 μ\mu 的设置及对增量方程的影响

3.2.1 置信区间

增加置信区间后，问题从无约束的最小二乘问题变为有约束的最小二乘问题：

min⁡x∈RnF(xk+Δxk)=12∣∣f(xk+Δxk∣∣2s.t.∣∣DΔxk∣∣2≤μ\underset{x\in\mathbb{R}^n}{\min}F(x_k+\Delta x_k)=\frac{1}{2}||f(x_k+\Delta x_k||^2\\ s.t.||D\Delta x_k||^2\le \mu

式中 ∣∣DΔxk∣∣2≤μ||D\Delta x_k||^2\le \mu 为置信区间的条件，DD 是系数矩阵，为 nn 维行向量，μ\mu 为置信半径。

3.2.2 增量方程

L-M 法的增量方程为

(H(xk)+λDTD)Δxk=g(xk)(H(x_k)+\lambda D^TD)\Delta x_k=g(x_k)

其中，H(xk),g(xk)H(x_k), g(x_k) 和高斯牛顿法中相同。

3.2.3 近似程度的量化

定义指标 ρ\rho 衡量近似程度：

ρ=f(xk+Δxk)−f(xk)J(xk)TΔxk\rho=\frac{f(x_k+\Delta x_k)-f(x_k)}{J(x_k)^T\Delta x_k}

  • 当 ρ\rho 接近 1 时，近似程度最好；
  • 当 ρ\rho 太小时，实际减小的值远小于近似函数减小的值，近似效果差，需要缩小近似范围 μ\mu；
  • 当 ρ\rho 较大时，实际减小的值大于近似函数减小的值，近似效果差，需要增大近似范围

3.3 L-M 法

3.3.1 列文伯格法

设置系数矩阵 D=InD=I_n。

3.3.2 马夸尔特法

系数矩阵 DD 为 J(xk)J(xk)TJ(x_k)J(x_k)^T 对角元素的平方根。

3.3.3 算法流程

3.3.3.1 算法

  1. 给定初始值 x0x_0 以及初始置信半径 μ\mu，令 k=0k=0；
  2. 对于第 kk 次迭代，在高斯牛顿法的基础上，添加置信区间约束 ∣∣DΔxk∣∣2≤μ||D\Delta x_k||^2\le \mu；
  3. 计算 ρ\rho；
  4. 若 ρ>34\rho>\frac{3}{4}，设置 μ=13μ\mu=\frac{1}{3}\mu，若 ρ<14\rho<\frac{1}{4}，设置 μ=2μ\mu=2\mu；
  5. 令 xk+1=xk+Δxkx_{k+1}=x_k+\Delta x_k；
  6. 判断算法是否收敛，若不收敛，则回到 2，否则结束输出 x∗≈xkx^*\approx x_k。

3.3.3.2 置信域调整方案一

  • 如果 ρ>34\rho>\frac{3}{4}，设置 μ=13μ\mu=\frac{1}{3}\mu；
  • 如果 ρ<14\rho<\frac{1}{4}，设置 μ=2μ\mu=2\mu。

3.3.3.3 置信域调整方案二

  • 如果 ρ>0\rho>0，则 μ=μ×max⁡{13,1−(2ρ−1)3}, ν=2\mu=\mu\times\max\{\frac{1}{3},1-(2\rho-1)^3\},\ \nu=2；
  • 否则，μ=μ×ν, ν=2ν\mu=\mu\times\nu,\ \nu=2\nu。

对于错误的 Δx\Delta x，也就是 ρ<0\rho<0 的情况，将会迅速增大阻尼系数，使下一步的 Δx\Delta x 趋向于更小。比方案一多了一个参数 ν\nu。

3.4 L-M 法与高斯牛顿法对比

列文伯格-马夸尔特在一定程度上解决了线性方程组系数矩阵的非奇异和病态问题，但算法收敛速度相较于高斯牛顿更慢。因此，当问题性质较好时，选择高斯牛顿方法；问题接近病态时，选择列文伯格-马夸尔特方法。

4 共轭梯度法

共轭梯度法用于解决无约束凸二次规划

min⁡f(x)=12xTQx+qTx+c\min f(x)=\frac{1}{2}x^TQx+q^Tx+c\\

式中 QQ 为对称正定矩阵，q∈Rnq\in\mathbb{R}^n，c∈Rc\in\mathbb{R}。

目标函数在最优点 x∗x^* 有一阶导数为 0，即

Qx∗+q=0Qx^*+q=0

这也意味着求解上述规划等价于求解方程

Qx=−qQx=-q

4.1 确定步长

根据迭代公式，有

f(xk+1)=f(xk+λkΔxk)=12(xk+λkΔxk)TQ(xk+λkΔxk)+qT(xk+λkΔxk)+cf(x_{k+1})=f(x_k+\lambda_k\Delta x_k)=\frac{1}{2}(x_k+\lambda_k\Delta x_k)^TQ(x_k+\lambda_k\Delta x_k)+q^T(x_k+\lambda_k\Delta x_k)+c

求导

df(xk+1)dλk=ΔxkTQ(xk+λkΔxk)+ΔxkTq=ΔxkTQxk+ΔxkTQλkΔxk+ΔxkTq=0ΔxkT(Qxk+q)=−ΔxkTQλkΔxkλk=−ΔxkT(Qxk+q)ΔxkTQΔxk=−Δxk∇f(xk)ΔxkTQΔxk\frac{\mathrm{d}f(x_{k+1})}{\mathrm{d}\lambda_k}=\Delta x_k^TQ(x_k+\lambda_k\Delta x_k)+\Delta x_k^Tq=\Delta x_k^TQx_k+\Delta x_k^TQ\lambda_k\Delta x_k+\Delta x_k^Tq=0\\ \Delta x_k^T(Qx_k+q)=-\Delta x_k^TQ\lambda_k\Delta x_k\\ \lambda_k=-\frac{\Delta x_k^T(Qx_k+q)}{\Delta x_k^TQ\Delta x_k}=-\frac{\Delta x_k\nabla f(x_k)}{\Delta x_k^TQ\Delta x_k}

4.2 寻找最优方向

4.2.1 共轭方向组

若 Q∈Rn×nQ\in \mathbb{R}^{n\times n} 为正定对称矩阵，d1,d2,...,dm∈Rnd_1,d_2,...,d_m\in \mathbb{R}^n，当 i≠ji\ne j 时有 diTQdj=0d_i^TQd_j=0，则 d1,d2,...,dmd_1,d_2,...,d_m 为关于矩阵 QQ 共轭，称其为共轭方向组。

  • 性质：共轭方向组线性无关。

4.2.2 最优方向

只需在共轭方向上寻找步长即可找到全局最优解（证明略）。将方向写为

Δxk=−∇f(xk)+γkΔxk−1\Delta x_k=-\nabla f(x_k)+\gamma_k\Delta x_{k-1}

将共轭方向组性质 Δxk−1TQΔxk=0\Delta x_{k-1}^TQ\Delta x_k=0 带入得

Δxk−1TQ(−∇f(xk)+γkΔxk−1)=0γk=Δxk−1TQ∇f(xk)Δxk−1TQΔxk−1\Delta x_{k-1}^TQ(-\nabla f(x_k)+\gamma_k\Delta x_{k-1})=0\\ \gamma_k=\frac{\Delta x_{k-1}^TQ\nabla f(x_k)}{\Delta x_{k-1}^TQ\Delta x_{k-1}}

4.2.3 证明 Δxi\Delta x_i 和 Δxj(j=0,1,...,i−1)\Delta x_j(j=0,1,...,i-1) 满足共轭关系

这里，为方便，改为证明 Δxi+1\Delta x_{i+1} 和 Δxj(j=0,1,..,i)\Delta x_j(j=0,1,..,i) 满足共轭关系：

ΔxjTQΔxi+1=ΔxjTQ(−∇f(xi+1)+γiΔxi)=ΔxjTQ(−∇f(xi+1))+γiΔxjTQΔxi\Delta x_j^TQ\Delta x_{i+1}=\Delta x_j^TQ(-\nabla f(x_{i+1})+\gamma_i\Delta x_i)=\Delta x_j^TQ(-\nabla f(x_{i+1}))+\gamma_i\Delta x_j^TQ\Delta x_i

令 ΔxjTQΔxi=0\Delta x_j^TQ\Delta x_{i}=0 可以得到 ΔxjTQΔxi+1=ΔxjTQ(−∇f(xi+1))\Delta x_j^TQ\Delta x_{i+1}=\Delta x_j^TQ(-\nabla f(x_{i+1}))，同乘最优步长 λi\lambda_i 得到

λkΔxjTQΔxi+1=λiΔxjTQ(−∇f(xi+1))=(−∇f(xi+1))TQλiΔxj\lambda_k\Delta x_j^TQ\Delta x_{i+1}=\lambda_i\Delta x_j^TQ(-\nabla f(x_{i+1}))=(-\nabla f(x_{i+1}))^TQ\lambda_i\Delta x_j

根据 xj+1=xj+γiΔxjx_{j+1}=x_j+\gamma_i\Delta x_j 得到：

(−∇f(xi+1))TQλiΔxj=(−∇f(xi+1))TQ(xj+1−xj)=(−∇f(xi+1))T(∇f(xj+1)−∇f(xj))(-\nabla f(x_{i+1}))^TQ\lambda_i\Delta x_j=(-\nabla f(x_{i+1}))^TQ(x_{j+1}-x_j)=(-\nabla f(x_{i+1}))^T(\nabla f(x_{j+1})-\nabla f(x_j))

根据 Δxj+1=−∇f(xj+1)+γjΔxj\Delta x_{j+1}=-\nabla f(x_{j+1})+\gamma_j\Delta x_j 得到

∇f(xj+1)−∇f(xj)=(γjΔxj−Δxj+1)−(γj−1Δxj−1−Δxj)(−∇f(xi+1))TQλiΔxj=(−∇f(xi+1))T((γjΔxj−Δxj+1)−(γj−1Δxj−1−Δxj))\nabla f(x_{j+1})-\nabla f(x_j)=(\gamma_j\Delta x_j-\Delta x_{j+1})-(\gamma_{j-1}\Delta x_{j-1}-\Delta x_j)\\ (-\nabla f(x_{i+1}))^TQ\lambda_i\Delta x_j=(-\nabla f(x_{i+1}))^T((\gamma_j\Delta x_j-\Delta x_{j+1})-(\gamma_{j-1}\Delta x_{j-1}-\Delta x_j))

由于 Δxi(i=0,1,...,n−1)\Delta x_i(i=0,1,...,n-1) 与 ∇f(xn)\nabla f(x_n) 正交（证明略），Δxi(i=0,1,...,n−1)\Delta x_i(i=0,1,...,n-1) 构成共轭方向组，对于共轭方向组 {Δxj,Δxj+1}\{\Delta x_j,\Delta x_{j+1}\}，有

(−∇f(xi+1))T(γjΔxj−Δxj+1)=0(−∇f(xi+1))T(γj−1Δxj−1−Δxj)=0(-\nabla f(x_{i+1}))^T(\gamma_j\Delta x_j-\Delta x_{j+1})=0\\ (-\nabla f(x_{i+1}))^T(\gamma_{j-1}\Delta x_{j-1}-\Delta x_j)=0

因此得到 ΔxjTQΔxi=0\Delta x_j^TQ\Delta x_{i}=0，从而证明共轭。

4.3 算法流程

  1. 初始化 x0, Δx0=∇f(x0), k=0x_0,\ \Delta x_0=\nabla f(x_0),\ k=0；
  2. 若 ∇f(xk)=0\nabla f(x_k)=0，返回 x∗≈xkx^*\approx x_k，结束算法，否则进入 3；
  3. 若 ∇f(xk)≠0\nabla f(x_k)\ne0，则令 xk+1=xk+λkΔxk, k=k+1x_{k+1}=x_k+\lambda_k\Delta x_k,\ k=k+1，返回 2。

5 拟牛顿法

为了解决牛顿法中海森矩阵 HH 计算复杂的问题，拟牛顿法提供了另外一种解决思路：通过使用不含有二阶导数的矩阵 UU 代替牛顿法中的 HH，根据矩阵 UU 构造的不同，具有不同的拟牛顿法。

5.1 拟牛顿法的基本原理

牛顿法搜索方向为

Δxk=−H(xk)−1g(xk)\Delta x_k=-H(x_k)^{-1}g(x_k)

现在构造一个矩阵 UU 来近似矩阵 HH，当 f(x)f(x) 为二次函数的时候，HH 为常数矩阵，且任意两点 xkx_k 和 xk+1x_{k+1} 之间的梯度差为：

∇f(xk+1)−∇f(xk)=H⋅(xk+1−xk)xk+1−xk=H−1⋅(∇f(xk+1)−∇f(xk))\nabla f(x_{k+1})-\nabla f(x_k)=H\cdot(x_{k+1}-x_k)\\ x_{k+1}-x_k=H^{-1}\cdot(\nabla f(x_{k+1})-\nabla f(x_k))

因此，当 f(x)f(x) 不是二次函数的时候，仿照上式，构建矩阵 UU 且要求满足下面的关系：

xk+1−xk=Uk+1⋅[∇f(xk+1)−∇f(xk)]x_{k+1}-x_k=U_{k+1}\cdot[\nabla f(x_{k+1})-\nabla f(x_k)]

此时可以按照下式来计算增量（也称为拟牛顿条件）：

Δxk=Uk+1−1⋅Δg(xk)\Delta x_k=U_{k+1}^{-1}\cdot\Delta g(x_k)

5.2 DFP

为方便区分，下面把 UU 称为 DD。

5.2.1 DFP 结果

在已知 DkD_k 的情况下，使用增量形式求取 Dk+1=ΔDk+DkD_{k+1}=\Delta D_k+D_k，构造 ΔDk\Delta D_k 如下：

ΔDk=ΔxkΔxkTΔg(xk)TΔxk−DkΔg(xk)Δg(xk)TDkΔg(xk)TDkΔg(xk)\Delta D_k=\frac{\Delta x_k\Delta x_k^T}{\Delta g(x_k)^T\Delta x_k}-\frac{D_k\Delta g(x_k)\Delta g(x_k)^TD_k}{\Delta g(x_k)^TD_k\Delta g(x_k)}

5.2.2 DFP 算法流程

  1. 给定初始点 x0x_0，允许误差 ε\varepsilon，令 D0=In, k=0D_0=I_n,\ k=0；
  2. 计算搜索方向 Δxk=−Dk−1⋅g(xk)\Delta x_k=-D_k^{-1}\cdot g(x_k)；
  3. 从 xkx_k 出发沿 gkg_k 做一维搜索，获得最优步长 λk=argmin⁡⁡λ∈Rf(xk+λΔxk)\lambda_k=\underset{\lambda\in\mathbb{R}}{\operatorname{arg\min}}f(x_k+\lambda\Delta x_k)，更新 xk+1=xk+λkΔxkx_{k+1}=x_k+\lambda_k\Delta x_k；
  4. 若 ∣∣g(xk+1)∣∣<ε||g(x_{k+1})||<\varepsilon 则停止迭代，输出 x∗≈xk+1x^*\approx x_{k+1}，否则转 5；
  5. 计算 Δg(xk), Δxk\Delta g(x_k),\ \Delta x_k，用上面的式子计算 Dk+1D_{k+1}；
  6. 令 k=k+1k=k+1，转 2。

5.3 BFGS

为方便区分，下面把 UU 称为 BB。

5.3.1 BFGS 结果

ΔBk=Δg(xk)Δg(xk)TΔxkTΔg(xk)−BkΔxkΔxkTBkΔxkTBkΔxk\Delta B_k=\frac{\Delta g(x_k)\Delta g(x_k)^T}{\Delta x_k^T\Delta g(x_k)}-\frac{B_k\Delta x_k\Delta x_k^TB_k}{\Delta x_k^TB_k\Delta x_k}

即将 DFP 中 Δx\Delta x 和 Δg\Delta g 调换位置。

5.3.2 谢尔曼 · 莫里森公式

对于任意非奇异方阵 AA，nn 维向量 u, vu,\ v，若 1+vTA−1u≠01+v^TA^{-1}u\ne0，则：

(A+uvT)−1=A−1−(A−1u)(vTA−1)1+vTA−1u(A+uv^T)^{-1}=A^{-1}-\frac{(A^{-1}u)(v^TA^{-1})}{1+v^TA^{-1}u}

通过使用两次谢尔曼 · 莫里森公式，可以迭代求取 BB 的逆：

Bk+1−1=(In−ΔxkΔg(xk)TΔxkTΔg(xk))Bk−1(In−Δg(xk)ΔxkTΔxkTΔg(xk))+ΔxkΔxkTΔxkTΔg(xk), B0=InB_{k+1}^{-1}=\Bigg(I_n-\frac{\Delta x_k\Delta g(x_k)^T}{\Delta x_k^T\Delta g_(x_k)}\Bigg)B_k^{-1}\Bigg(I_n-\frac{\Delta g(x_k)\Delta x_k^T}{\Delta x_k^T\Delta g(x_k)}\Bigg)+\frac{\Delta x_k\Delta x_k^T}{\Delta x_k^T\Delta g(x_k)},\ B_0=I_n

5.3.3 算法流程

  1. 给定初始点 x0x_0，允许误差 ε\varepsilon，令 B0=In, k=0B_0=I_n,\ k=0；
  2. 计算搜索方向 Δxk=−Bk−1⋅g(xk)\Delta x_k=-B_k^{-1}\cdot g(x_k)；
  3. 从 xkx_k 出发沿 gkg_k 做一维搜索，获得最优步长 λk=argmin⁡⁡λ∈Rf(xk+λΔxk)\lambda_k=\underset{\lambda\in\mathbb{R}}{\operatorname{arg\min}}f(x_k+\lambda\Delta x_k)，更新 xk+1=xk+λkΔxkx_{k+1}=x_k+\lambda_k\Delta x_k；
  4. 若 ∣∣g(xk+1)∣∣<ε||g(x_{k+1})||<\varepsilon 则停止迭代，输出 x∗≈xk+1x^*\approx x_{k+1}，否则转 5；
  5. 计算 Δg(xk), Δxk\Delta g(x_k),\ \Delta x_k，用上面的式子计算 Bk+1−1B_{k+1}^{-1}；
  6. 令 k=k+1k=k+1，转 2。

需要注意的是，BFGS 从头至尾都不用计算矩阵的逆。

5.4 L-BFGS

5.4.1 L-BFGS 原理

在BFGS 中需要保存 B−1B^{-1}，而实际上，仅保存每轮的 Δx\Delta x 和 Δg\Delta g 即可递归的计算当前的 B−1B^{-1}，从而减小内存占用

Bk+1−1=(∏i=k0ViT)B0−1(∏i=0kVi)+∑j=0k(∏i=kj+1ViT)(ρjΔxjΔxjT)(∏i=j+1kVi)ρk=1ΔxkTΔg(xk)Vk=In−ρkΔg(xk)ΔxkTB_{k+1}^{-1}=\Bigg(\prod_{i=k}^{0}V_i^T\Bigg)B_0^{-1}\Bigg(\prod_{i=0}^{k}V_i\Bigg)+\sum_{j=0}^k\Bigg(\prod_{i=k}^{j+1}V_i^T\Bigg)(\rho_j\Delta x_j\Delta x_j^T)\Bigg(\prod_{i=j+1}^{k}V_i\Bigg)\\ \rho_k=\frac{1}{\Delta x_k^T\Delta g(x_k)}\\ V_k=I_n-\rho_k\Delta g(x_k)\Delta x_k^T

实际应用中，可以仅使用最近 mm 个 Δx\Delta x 和 Δg\Delta g 来计算

Bk−1=(∏i=k−1k−mViT)B0−1(∏i=k−mk−1Vi)+∑j=k−1k−m(∏i=kj+1ViT)(ρjΔxjΔxjT)(∏i=j+1kVi)B_{k}^{-1}=\Bigg(\prod_{i=k-1}^{k-m}V_i^T\Bigg)B_0^{-1}\Bigg(\prod_{i=k-m}^{k-1}V_i\Bigg)+\sum_{j=k-1}^{k-m}\Bigg(\prod_{i=k}^{j+1}V_i^T\Bigg)(\rho_j\Delta x_j\Delta x_j^T)\Bigg(\prod_{i=j+1}^{k}V_i\Bigg)

5.4.2 算法流程

  1. 给定初始点 x0x_0，允许误差 ε\varepsilon，预定保留最近 mm 个向量，设置 B0−1, t=0B_0^{-1},\ t=0；
  2. 计算搜索方向 Δxk=−Bk−1⋅g(xk)\Delta x_k=-B_k^{-1}\cdot g(x_k)；
  3. 从 xkx_k 出发沿 gkg_k 做一维搜索，获得最优步长 λk=argmin⁡⁡λ∈Rf(xk+λΔxk)\lambda_k=\underset{\lambda\in\mathbb{R}}{\operatorname{arg\min}}f(x_k+\lambda\Delta x_k)，更新 xk+1=xk+λkΔxkx_{k+1}=x_k+\lambda_k\Delta x_k；
  4. 若 ∣∣g(xk+1)∣∣<ε||g(x_{k+1})||<\varepsilon 则停止迭代，输出 x∗≈xk+1x^*\approx x_{k+1}，否则转 5；
  5. 判断 k>mk>m，删掉存储的 Δxk−m\Delta x_{k-m} 和 Δg(xk−m)\Delta g(x_{k-m})；
  6. 计算 Δg(xk), Δxk\Delta g(x_k),\ \Delta x_k，令 k=k+1k=k+1，转 2。

6 直线拟合——最小二乘

给定一些列 xi,yi(i=1,2,...,N)x_i,y_i(i=1,2,...,N)，假定它们具有线性关系，既可以用 y=kx+by=kx+b 的方式进行拟合，任务是找到最优的 kk 和 bb 使得 xi,yix_i,y_i 和直线y=kx+by=kx+b 最接近（一元回归）。此时问题就变为了求解目标规划

min⁡f=∑i=1N(yi−kxi−b)2\min f=\sum_{i=1}^N(y_i-kx_i-b)^2

使用线性代数的表达形式：

{y1=kx1+by2=kx2+b...yN=kxN+b\left\{ \begin{matrix} y_1=kx_1+b\\ y_2=kx_2+b\\ ...\\ y_N=kx_N+b \end{matrix} \right.

转为矩阵形式

[y1y2...yN]=[x11x21......xN1][kb]\left[ \begin{matrix} y_1\\ y_2\\ ...\\ y_N \end{matrix} \right] = \left[ \begin{matrix} x_1&1\\ x_2&1\\ ...&...\\ x_N&1 \end{matrix} \right] \left[ \begin{matrix} k\\ b \end{matrix} \right]

令

Y=[y1y2...yN], K=[kb], X=[x11x21......xN1]Y= \left[ \begin{matrix} y_1\\ y_2\\ ...\\ y_N \end{matrix} \right],\ K= \left[ \begin{matrix} k\\ b \end{matrix} \right],\ X= \left[ \begin{matrix} x_1&1\\ x_2&1\\ ...&...\\ x_N&1 \end{matrix} \right]

则可以写成矩阵表达式

Y=XKY=XK

此时，为求解 KK，若 XX 可逆，那么可以直接求得，然而若不可逆，则可以通过如下的方式求解：

\Rightarrow& Y=XK\\ \Rightarrow& X^TY=X^TXK\\ \Rightarrow& (X^TX)^{-1}X^TY=(X^TX)^{-1}X^TXK