非线性最优化python 非线性最优化计算方法

转载

代码匠人之心 2023-12-26 17:47:46

文章标签 非线性最优化python 牛顿法迭代最小二乘 文章分类 Python 后端开发

非线性最小二乘

定义：简单的非线性最小二乘问题可以定义为

minx12||f(x)||22minx12||f(x)||22

其中自变量x∈Rnx∈Rn,f(x)f(x)是任意的非线性函数，并设它的维度为mm,即f(x)∈Rmf(x)∈Rm.

对于一些最小二乘问题，我们可以利用目标函数对xx求导并令导数等于0来求解。但是导数

d(12||f(x)||22)dx=0d(12||f(x)||22)dx=0

不一定可以直接求解xx,这个导函数可能是一个复杂的非线性方程。这种情况下一般采用迭代来求解，具体步骤可以表示为 :

(1) 给定一个初试值x0x0
(2) 对于第kk次迭代，寻找一个增量ΔxkΔxk,使得||f(xk+Δxk)||||f(xk+Δxk)||达到极小值
(3) 若ΔxkΔxk足够小，则停止迭代
(4) 否则，令xk+1=xk+Δxkxk+1=xk+Δxk,返回第(2)步骤

这个其实是通过迭代让目标函数一步步下降，直到最终达到收敛为止。一般而言，增量ΔxΔx可通过一阶梯度或二阶梯度来确定。

一阶梯度和二阶梯度法

首先，我们将目标函数在xx附近进行泰勒展开

||f(x+Δx)||22≈||f(x)||22+J(x)Δx+12ΔxTH(x)Δx||f(x+Δx)||22≈||f(x)||22+J(x)Δx+12ΔxTH(x)Δx

这里的J(x)J(x)是f(x)f(x)关于xx的导数（雅可比矩阵），H(x)H(x)是二阶导数(海森矩阵)。我们可以选择保留泰勒公式的一阶导数和二阶导数，如果保留一阶导数，则增量的解就是

Δx∗=−JT(x)Δx∗=−JT(x)

它理解起来就是，沿着沿着梯度相反的方向前进，目标函数下降得最快。通常会在这个方向上计算一个步长λλ,迭代公式表示为

xk+1=xk−λJT(xk)xk+1=xk−λJT(xk)

该方法称为最速梯度下降法。如果保留二阶梯度信息，增量可以表示为

Δx∗=argminΔx||f(x)||22+J(x)Δx+12ΔxTH(x)ΔxΔx∗=arg⁡minΔx||f(x)||22+J(x)Δx+12ΔxTH(x)Δx

对ΔxΔx求导数并令它等于0，则

JT+HΔx=0JT+HΔx=0

于是增量的解为

HΔx=−JTHΔx=−JT

这种方法称为牛顿法，它的迭代公式可以表示为

xk+1=xk−H−1Jxk+1=xk−H−1J

牛顿法和最速梯度下降法思想比较简单，只需要将函数在迭代点附近展开，然后对增量求最小化即可，然后可以通过线性方程直接求的增量的解。这两种方法的主要缺点为：

最速梯度下降过于贪心，容易走出锯齿状，反而增加了迭代次数
牛顿法则需要计算函数的HH矩阵，这在问题规模较大时候非常困难，通常做法是避免去计算HH

下面介绍了高斯牛顿法和列文伯格-马夸尔特法（LM）法更为实用。

高斯牛顿法

高斯牛顿法的思想是对f(x)f(x)进行一阶泰勒展开，注意不是目标函数||f(x)||22||f(x)||22

f(x+Δx)≈f(x)+J(x)Δxf(x+Δx)≈f(x)+J(x)Δx

其中J(x)J(x)是f(x)关于xx的导数，称为雅可比矩阵，维度为m×nm×n.于是我们需要

Δx∗=argminΔx12||f(x)+J(x)Δx||22Δx∗=arg⁡minΔx12||f(x)+J(x)Δx||22

我们继续对ΔxΔx求导并令导函数等于0

∇Δx12||f(x)+J(x)Δx||22=∇Δx12[fT(x)f(x)+2fT(x)J(x)Δx+ΔxTJT(x)J(x)Δx]=JTf(x)+JTJΔx=0∇Δx12||f(x)+J(x)Δx||22=∇Δx12[fT(x)f(x)+2fT(x)J(x)Δx+ΔxTJT(x)J(x)Δx]=JTf(x)+JTJΔx=0

从而

JTJΔx=−JTf(x)JTJΔx=−JTf(x)

记H=JTJ,g=−JTf(x)H=JTJ,g=−JTf(x),则

HΔx=gHΔx=g

上述方程称为高斯牛顿方程或正规方程。对比牛顿法可见，高斯牛顿法采用JTJJTJ牛顿法作为牛顿法HH矩阵的近似，从而避免了复杂的计算。原则上，它要求近似的矩阵HH是可逆的（而且是正定的），而实际计算中得到的JTJJTJ却是半正定的。也就是使用高斯牛顿法会出现JTJJTJ为奇异或者病态情况，此时增量的稳定性较差，导致算法不收敛。即使HH非奇异也非病态，如果求得的ΔxΔx非常大,也会导致我们采用的局部近似不够正确，这样以来可能不能保证收敛，哪怕是还有可能让目标函数更大。

即使高斯牛顿法具有它的缺点，但是很多非线性优化可以看作是高斯牛顿法的一个变种，这些算法结合了高斯牛顿法的优点并修正其缺点。例如LM算法，尽管它的收敛速度可能比高斯牛顿法更慢，但是该方法健壮性更强，也被称为阻尼牛顿法。

LM算法

高斯牛顿法采用二阶泰勒展开来近似，只有在展开点附近才会有比较好的近似效果。LM(Levenberg -Marquard）算法中给变化量ΔxΔx添加一个信赖区域来限制ΔxΔx的大小，并认为在信赖区域里面近似是有效的，否则近似不准确。

确定信赖区域一个好的办法是通过比较近似模型和实际模型的差异来确定，如果差异小，我们就让范围尽可能增大；如果差异太大，就该缩小这个范围。考虑实际模型和近似模型变化量的比值

ρ=f(x+Δx)−f(x)J(x)Δxρ=f(x+Δx)−f(x)J(x)Δx

上面式子可以通过ρρ的值来判断泰勒近似的好坏，其中分子是实际模型的变化量，分母是近似模型的变化量，当

ρρ接近1的时候表明近似模型是非常好的，如果ρρ较小，则实际模型的变化量小于近似模型的变化量，则认为近似模型较差，需要缩小近似范围；反之，当ρρ较大时，说明实际模型变化量更大，我们需要放大近似范围。因此LMLM算法可以表示如下

（1）给定初始迭代值x0x0以及初始优化半径μμ.
（2）对于第kk次迭代，求解优化问题
minΔxk12||f(xk)+J(xk)Δx||2s.t.||DΔxk||2≤μ.minΔxk12||f(xk)+J(xk)Δx||2s.t.||DΔxk||2≤μ.
（3）计算ρρ.
（4）若ρ>3/4ρ>3/4,则μ=2μμ=2μ.
（5）若ρ<1/4ρ<1/4,则μ=0.5μμ=0.5μ.
（6）如果ρρ大于某阈值，则认为近似可行。令xk+1=xk+Δxkxk+1=xk+Δxk.
（7）判断算法是否收敛。如果不收敛，跳回步骤（2），否则结束。

上面算法中，μμ是信赖区域的半径，其实是一个球形的区域，该约束认为只有在球内才是有效的，带上矩阵DD后就是一个椭球。

采用拉格朗日乘子将上述问题转化为一个无约束问题

minΔxk12||f(xk)+J(xk)Δxk||2+12λ||DΔxk||2minΔxk12||f(xk)+J(xk)Δxk||2+12λ||DΔxk||2

对ΔxkΔxk求梯度，得到增量方程

(H+λDTD)Δx=g(H+λDTD)Δx=g

可以看到，增量方程在搞死牛顿法上多了一项λDTDλDTD,如果令D=ID=I，则简化版本表示为

(H+λI)Δx=g(H+λI)Δx=g

当λλ较小时，说明HH占主导地位，说明二次近似在该范围类是比较好的，LM方法更接近于高斯牛顿法；另一方面，当λλ较大时，λλ占主导地位，LM算法更接近一阶梯度下降算法，这说明二次近似不够好。LM算法的求解方式，可以避免线性方程组的矩阵非奇异和病态等问题，提供更稳定、更准确的解法。