0%

非线性最小二乘(Nonlinear Least Squares)

本文对非线性最小二乘的原理、思路和应用做了简单介绍。

非线性最小二乘

最小二乘是一种数学最优化建模方法,它的一种常见解释是对残差做了满足正态分布的极大似然估计。最小二乘分两类:线性最小二乘和非线性最小二乘,区别是残差函数是否为线性函数,线性最小二乘比较简单,可以直接求解析解,非线性最小二乘比较复杂,需要用迭代式的最优化方法求解。

最小二乘回顾

假设有观测数据集{xi,yi}i=1n,某个参数为w的函数f(x;w)描述了观测数据的因变量和自变量的映射关系,即: y^i=f(xi;w) 但数据拟合时一定不是严丝合缝的,会有误差即残差,所以有: yi=y^i+ϵ 假设这个残差服从偏差为0的正态分布,则有: ϵN(0,σ2) 于是有: p(yi|xi,w)=12πσe(yif(xi;w))22σ2

对于这n个独立的观测样本就会产生n个方程,采用极大似然估计有: p(y|x,w)=i=1n12πσe(yif(xi;w))22σ2=1(2πσ)nei=1n(yif(xi;w))22σ2

由贝叶斯公式可知: p(w|x,y)=p(y|x,w)p(w|x)p(y|x)=p(w|x)p(y|x)p(y|x,w)p(y|x,w) 所以,利用极大似然求参有: argmaxwp(w|x,y)=argminwln(p(y|x,w))=nlnσ+i=1n(yif(xi;w))22σ212σ2i=1n(yif(xi;w))2 最后有: w=argminwi=1n(yif(xi;w))2

非线性最小二乘

当残差函数非线性时,方便起见,我们定义残差函数为:r(w)=yf(x;w),对n个观测样本有n个方程: F(w)=[r1(w),r2(w),......,rn(w)]T 非线性最小二乘问题变成: min||r(w)||2=r1(w)2+r2(w)2+......rn(w)2 如果此时还有其他约束条件,例如,等式约束q(w)=0的话,以上问题变成了约束非线性最小二乘问题: min||r(w)||2=r1(w)2+r2(w)2+......rn(w)2s.t.q1(w)=0......qm(w)=0

拉格朗日乘数法

考虑求解通用的约束最优化问题,以最小化问题为: minf(x)s.t.q1(x)=0......qm(x)=0 几个定义:

1、xn维向量,当x满足以上所有约束条件时,被叫做可行解;

2、所有可行解中,当满足f(x)f(x),则x被称作全局最优解;

3、当可行解只有在一个邻域R中时才满足f(x)f(x),则x被称作局部最优解,其中||xx||R

为求解上述约束问题,需要将其转变为无约束问题,在机器学习与人工智能技术分享-第四章 最优化原理中有过介绍,定义拉格朗日函数: L(x,λ)=f(x)+λ1q1(x)+λ2q2(x)+...+λmqm(x)=f(x)+λTq(x) 其中:λ=(λ1,λ2,...,λm)m维的拉格朗日乘子向量,q(x)=(q1(x),q2(x),...,qm(x))m维的约束函数向量。

拉格朗日函数的梯度为: L(x,λ)=[xL(x,λ)λL(x,λ)] 其中: xL(x,λ)=f(x)+λ1q1(x)+...+λmqm(x)=f(x)+Dq(x)TλλL(x,λ)=q(x)

如果可行解x为局部最优解,则必须满足拉格朗日条件(first-order necessary conditions),即此条件是必要条件: f(x)+Dq(x)Tλ=0q(x)=0

其中: 𝐷q(𝑥)=[q1(x)x1q1(x)xnqm(x)x1qm(x)xn]=[q1(x)Tqm(x)T]

需要注意,如果在可行解x下,𝐷q(𝑥)={q1(x),...,gm(x)}是线性相关的,则拉格朗日乘子不存在,因此无法通过拉格朗日乘数法求解,但此时局部最优解可能是存在的; 反之,如果𝐷q(𝑥)线性无关,则x被称为正则点(regular point),举个例子: minf(x1,x2)=x1+x2s.t.q1(x1,x2)=(x11)2+x221=0q2(x1,x2)=(x12)2+x224=0

x=(0,0)为唯一的可行解,因此它也是最优解,其拉格朗日函数为: L(x,λ)=x1+x2+λ1((x11)2+x221)+λ2((x12)2+x224)x=(0,0)处的拉格朗日条件为: xL(x,λ)=[11]+2λ1[x11x2]+2λ2[x12x2]=0λL(x,λ)=[(x11)2+x221=0(x12)2+x224=0] 显然,λ1λ2都不存在,x=(0,0)不是正则点(q1(x)q2(x)是线性相关的),无法使用拉格朗日乘数法求解该问题。

约束非线性最小二乘

min||r(w)||2s.t.q(w)=0

其中w为待求参数,r(w)=(r1(w),...,rn(w))n个样本的残差,q(w)=(q1(w),...,qm(w))m个约束条件。

转化为拉格朗日函数:L(w,λ)=||r(w)||2+λTq(w),由拉格朗日条件知,如果w为最优解,则存在λ使得: 2Dr(w)Tr(w)+Dq(w)Tλ=0,q(w)=0

ps:Dq(w)的行向量为线性无关组。

r(w)q(w)都是线性函数时,该问题可以直接求出解析解,否则需要采用迭代法,在机器学习与人工智能技术分享-第四章 最优化原理中我们介绍了由泰勒展开式衍生出的各种一阶和二阶优化算法,例如各种梯度下降法和各种拟牛顿法, 简单地说,优化算法就是在研究怎么找方向和找步长,先找方向再找步长是Line Search方法,先找步长再找方向(也可以看做同时找步长和方向)是Trust Region方法。

Newton Step-based迭代法

前面也说到泰勒展开式是大部分最优化方法的“根”,如果展开到二阶,则可以变成各种牛顿法,简单回顾: 求解无约束最小化问题: minf(x) 在当前迭代xk点下,做泰勒展开到二阶: f(xk+Δxk)=f(xk)+f(xk)TΔxk+12ΔxkT2f(xk+tΔxk)Δxk 其中:t=(0,1)是个缩放因子,迭代式为:xk+1=xk+Δxk

简单起间用fk代替f(xk),用gk代替f(xk),用pk代替Δxk有: f(xk+pk)=fk+gkTpk+12pkT2f(xk+tpk)pk 由于Hessian矩阵Hk=2f(xk)计算成本很高,所以用一个矩阵Bk去估计Hk,于是上面问题变成: mk(pk)=fk+gkTpk+12pkTBkpk 其中:mk(pk)f(xk+pk)之间的误差是||pk||2的高阶无穷小,即:O(||pk||2),当pk很小时这个误差就很小。 由最优值必要条件可知,在第k轮迭代有: mk(pk)=gk+Bkpk=0Bkpk=gk

1、如果Bk=HkHk可逆、正定,上面问题就变成了牛顿法;

2、如果BkHkBk始终保持正定,就是拟牛顿法,它和牛顿法一样,都具有二次收敛速度,即: 当xk靠近最优值x时,有:||xk+1x||<ϵ||xkx||2ϵ<1

ps:只有Bk始终保持正定才能保证收敛,因为,对于最小化问题,始终希望搜索方向是一个下降方向,即满足pkTgk<0,所以: pkTgk<0Bkpk=gkpkTBkpk>0

3、如果Bk=I,显然单位矩阵不是Hessian矩阵的一个好的估计,但它够简单,于是有:pk=gk,这就是妥妥的梯度下降法了,这时的收敛速度只能到线性速度了,即:||xk+1x||<ϵ||xkx||ϵ<1

4、如果Hk非正定,需要重新定义Bk=Hk+λI,总能找到λ使得Bk正定,即Bkpk=gk(Hk+λI)pk=gk,从而总能找到最优pk,换个角度看这个形式:

  • λ很小时,如果Hk正定,上式蜕变为牛顿法;
  • λ很大时,出现两个现象: 1、上式蜕变为梯度下降法 2、||pk||会被减小

这个优化框架其实就是原始的Levenberg-Marquardt迭代法,即:当前值在距离最优值比较远时通过梯度下降法优化,虽然速度很慢但保证收敛,接近最优值时通过牛顿法以二次收敛速度逼近。 实际上可以证明,即使Bk不正定,只要对||pk||做合适的限制,就能做到全局收敛,这个思想就是Trust Region(信任域)法。 回看这个问题: minmk(pk)=fk+gkTpk+12pkTBkpkBk=Hk时,相当于泰勒展开式展开到3阶,此时误差为O(||pk||3),当||pk||很小时,这个误差非常非常小,所以确定||pk||这一步可以转化为一个子优化问题: minmk(pk)=fk+gkTpk+12pkTBkpks.t.||pk||Δk 其中Δk>0,是Trust Region的半径,定义||.||为欧氏距离,则最优解p在一个半径为Δk的球体内(把约束变个形式就清楚了:||pk||ΔkpkTpkΔk2),特别是当Hessian矩阵性质不太好(如出现奇异或者病态)时,Trust Region法却能很好的求解。

ps:由于目标函数和约束条件都是二次方形式,所以求解的计算成本比较低。

总结对比下Levenberg-Marquardt与Trust Region的区别

当牛顿法不好使时:

Levenberg-Marquardt:通过对Hk做扰动,求解(Hk+λI)pk=gk中的pk,从而做迭代;

Trust Region:通过求解子优化问题minmk(pk)=fk+gkTpk+12pkTBkpks.t.||pk||Δk得到最优pk,从而做迭代,让当前点xk迅速进入”兴趣域“。

所以Trust Region法可以看做是Levenberg-Marquardt法演化后得到的算法,另外,前者可以很好的处理目标函数是Negative Curvature的情况,后者不行,所以不要将两类算法混为一谈。

到此为止,Trust Region还有个遗留问题是Δk该怎么取? 给定pk,定义: ρk=f(xk)f(xk+pk)mk(0)mk(pk)

观察上述定义:

1、当ρk<0时,由于f(xk+pk)大于当前值,所以目标函数值上升了,此时的pk必须丢弃;

2、如果ρk接近1,说明mk估计的比较好,此时就可以放心的再下次迭代时扩展信任域;

3、如果ρk>0且远小于1,则下次迭代保持信任域不变;

4、如果ρk接近0,则需要在下次迭代时缩减信任域。

Trust Region伪码如下:

Trust Region Algorithm给定上界Δ^>0,Δ0(0,Δ^),η[0,14)for k in range(0,..)求解minmk(pk)=fk+gkTpk+12pkTBkpks.t.||pk||Δk得到最佳pk计算ρkif ρk<14Δk+1=14Δkelse if ρk>34 且 ||pk||=ΔkΔk+1=min(2Δk,Δ^)else Δk+1=Δkif ρk>ηxk+1=xk+pkelse xk+1=xkend 

Trust Region框架应用

1、问题变换

回到之前的问题,为方便后续推理,去掉等式约束(因为通过拉格朗日乘数法可以很容易变成无约束问题),加入迭代轮数标识k,并给目标函数加个系数,即: min12||r(wk)||2 由泰勒展开式可知: r(wk+pk)=12||r(wk+pk)||212||Dr(wk)pk+r(wk)||2

方便起见用符号Jk表示Jacobian矩阵Dr(wk)rk表示r(wk)qk表示q(wk)r(wk+pk)12||Jkpk+rk||2

即待求解的无约束问题从 min12||r(wk)||2 变成: min12||Jkpk+rk||2

利用Trust Region 框架,如果pk为上述问题的最优解,则一定存在λ0,且满足以下条件: (JkTJk+λI)pk=JkTrkλ(Δk||pk||)=0

换个角度看,(JTJ+λI)pk=JTrk其实就是求解下面这个最小二乘问题: minp12||[JkλI]p+[rk0]||2

2、Trust Region框架下描述解法

给定上界Δ^>0,Δ0(0,Δ^),η[0,14)for k in range(0,..)求解minmk(pk)=12||Jkpk+rk||2+gkTpk+12pkTBkpks.t.||pk||Δk1、当pk满足约束条件||pk||<Δk:则采用Gauss–Newton法有:mk(pk)=0JkTJkpk=JkTrk这种方法有几个好处:a.不像Newton法需要计算Hessian矩阵,从而大大提高计算效率b.用JTJ估计Hessian矩阵使得其收敛速度接近Newton法c.只要Jk满秩且目标函数梯度不为0,pk一定是目标函数的下降方向d.可以利用线性最小二乘法寻找搜索方向2、当pk满足约束条件||pk||=Δk:则有:(JkTJk+λI)pk=JkTrkλ(Δk||pk||)=0λ0得到最佳pk计算ρkif ρk<14Δk+1=14Δkelse if ρk>34 且 ||pk||=ΔkΔk+1=min(2Δk,Δ^)else Δk+1=Δkif ρk>ηxk+1=xk+pkelse xk+1=xkend 

罚函数法

对问题: min||r(w)||2s.t.q(w)=0 不再强制要求约束条件严格遵守:q(w)=0,而是增惩罚系数,变成无约束问题: min||r(w)||2+μ||q(w)||2,μ>0. 在每轮迭代k中,通过求解: min||r(w)||2+μk||q(w)||2,μk>0. 利用Levenberg–Marquardt算法更新wk+1,方便起见用符号Jk表示Jacobian矩阵Dr(wk)rk表示r(wk)qk表示q(wk),根据拉格朗日条件有: 2Jkrk+2μk1DqkTqk=0zk=2μk1qk,则有: 2Jkrk+DqkTzk=0||qk||足够小时算法结束迭代(此时μk1足够大)。

算法缺点:

为了让qk尽快逼近到0,需要惩罚系数μ的值快速增大,这个会带来:

1、求解非线性最小二乘子问题变得困难

2、会导致Levenberg–Marquardt算法需要更多的迭代甚至失败

增广拉格朗日法

对问题: min||r(w)||2s.t.q(w)=0 其最优值满足拉格朗日条件: 2Dr(w)r(w)+Dq(w)Tλ=0,q(w)=0 定义增广拉格朗日函数: Lμ(w,λ)=L(w,λ)+μ||q(w)||2=||r(w)||2+q(w)Tλ+μ||q(w)||2=||r(w)||2+μ||q(w)+12μλ||214μ||λ||2μ>0 根据拉格朗日条件知最优值w满足: 2Dr(w)r(w)+2μDq(w)T(q(w)+12μλ)=02Dr(w)r(w)+Dq(w)T(2μq(w)+λ)=0λ=λ+2μq(w) 有: 2Dr(w)r(w)+Dq(w)Tλ=0 显然,根据拉格朗日条件可知,如果w是最优解,必须满足q(w)=0,当它的值较大时,可以通过上述λ的迭代式进行迭代。完整表述如下: 1.利用LM算法在wk点开始迭代求解问题:min||r(wk)||2+μk||q(wk)+12μkλk||2,得到wk+12.更新迭代式:λk+1=λk+2μkq(wk+1)3.更新惩罚系数:μk+1={μk=if||q(wk+1)||<14||q(wk)||2μk=otherwise 迭代初始点可以选择:λ1=0,μ1=1,和罚函数法不同,μ只在需要时增加,增长速度缓慢,当q(wk)足够小或者达到最大迭代次数后算法结束。

汽车非线性控制问题

阿克曼转向技术

很长时间人们在车辆(如:马车)转向上采用单铰链转向技术,如下图,这种方式两个前轮的转向角是相同的,内外转向轮指向同一圆心,结构很简单,但转向摆动空间大时会影响稳定性,转弯角度大时也会发生卡死,就像小时候开玩具车时,转弯时一个不小心就人仰马翻。
为了避免上述问题,人们发明了一种转向结构,让内侧转向轮的转角比外侧转向轮大了约2~4度,并让四个轮子路径的圆心大致交会于后轴的延长线上的同一个中心点,从而能够平滑稳定的转向,后来这项技术被英国代理商阿克曼申请了专利,后被称作阿克曼转向技术(Ackermann Steering)如下图:

Simple Car模型

对阿克曼转向做简化,可以构建对车辆控制的简单模型:

车辆在任意一个时刻t的运动模型为: f(p1(t),p2(t),v(t),ϕ(t),θ(t)) 其中pi(t),i{1,2}t时刻位置,v(t)t时刻车辆后轴中点的瞬时速度,ϕ(t)t时刻前轴中点转向角,θ(t)t时刻车身朝向。

内外侧轮转向角的关系为: 内侧转向角:tanϕi=LRW2外侧转向角:tanϕo=LR+W2前轴中点转向角:tanϕ=LR

t时刻后轴中点运动位置的微分方程为:

dp1dt=v(t)cosθ(t)dp2dt=v(t)sinθ(t)

由弧长和角度关系可知,假设st时刻车辆行驶距离,则有微分方程: ds=Rdθ=Ltanϕdθvdt=Ltanϕdθdθdt=vtanϕL

所以t时刻总的的微分方程为: dp1dt=v(t)cosθ(t)dp2dt=v(t)sinθ(t)dθdt=v(t)Ltanϕ(t) 以上方程在本文叫做“车辆运动模型”。

Dubins Vehicle模型

我们思考一个问题:停车场有A、B两个位置,一个人站在A位置,在速度等约束条件同等的情况下,怎么走能最快到达B位置?答案是显而易见的,两点之间线段最短!如果把人换成车答案又是什么?由于车辆有长度和车身朝向,从一个位置到另一个位置需要依据转向角转向后才能向前行进,所以行驶的最短路径由直线段和圆弧段组成,这里最经典的方法是Dubin’s Car,这个方法最大的优点是求解最短路径规划用最简单的几何方法就能搞定,无需复杂的矩阵计算,而且Dubins在1957年完成了数学理论上的完备性证明。Dubin’s Car的控制只有三种操作:L,左转前进;R,右转前进;S,向前直行,则控制组合总共有6种: RSR、LSL、RSL、LSR、RLR、LRL,如果把L和R统一拿C(Curve)来表示的话,组合只有两种:CSC(RSR、LSL、RSL、LSR)和CCC(RLR、LRL)。

不管哪种控制方式,都有一个要遵循的独特切线。这条切线构成了“S”部分的轨迹。直线与圆相切的点是车辆必须经过的点,没有它们车辆无法完成其轨迹。所以解决这些轨迹的问题就可以归结为如何正确计算这些切线。 计算切线的方法有几何法和向量法两种,其中向量法效率更高,这里介绍下:

1、如何求两个圆的切点和切线

已知:两个圆C1C2及其圆心位置分别为p1=(x1,y1)p2=(x2,y2),以及它们的半径r1r2,它们不一定相等

求:两个圆的切点和切线

整个过程分以下几步:

  • 画出两个圆C1C2
  • 连接两个圆的圆心p1p2,得到向量V1=(x2x1,y2y1)||V1||=D
  • 连接两个圆的切点pot1pot2,得到切线V2=pot2pot1
  • 画出向量V2的单位方向向量,用n^表示
  • 以上几个向量关系为:V2n^=0,n^n^=1
  • 求切线的关键点是如何求得n^
  • 修改向量V1为向量V1,使得它平行于向量V2V1=V1(r2r1)n^,于是有关系: n^(V1(r2r1)n^)=0n^V1r2+r1=0n^V1=r2r1n^V1D=r2r1DV1^=V1Dc=r2r1D则有: n^V1^=cn^n^=1V1^V1^=1 由点积的定义: A^B^=||A||×||B||cos(θ) 可知,向量V1^n^的夹角为常量c,即: cos(θ)=csin(θ)=1c2 这里的未知量只有:n^,问题就变成了:把已知向量V1^=(v1x,v1y)旋转角度θ后得到向量n^=(nx,ny)。 旋转向量操作比较简答,即: nx=v1x×cv1y×1c2ny=v1x×1c2+v1y×c
  • 计算得到n^后,可以很方便的求出两个切点:从C1的圆心沿着n^方向,距离为r1的点为C1的切点pot1,然后从pot1出发,沿着V2方向,与C2圆相交,得到切点pot2

2、计算CSC曲线

CSC包括RSR、LSL、RSL、LSR,计算切点就是为了找到CSC曲线的“S”部分,计算CSC路径就变成了:从起点以最小转半径转弯到达第一个切点,直行到达第二个切点后,以最小转弯半径再次转弯,到达终点

以RSR为例,假设起点为:s=(x1,y1,θ1),终点为:g=(x2,y2,θ2)C1C2圆心计算方法为: pc1=(x1+rmincos(θ1π2),y1+rminsin(θ1π2))pc2=(x2+rmincos(θ2π2),y2+rminsin(θ2π2))

得到两个圆心后,就可以计算两个切点,从而得到车辆控制的三段轨迹:起始点到pot1的圆弧,pot1pot2的线段和pot2到终点的圆弧。

假设车辆速度恒定为1,起点为ps=(x,y,θ),终点为peT是从起点到达终点的时间,则根据车辆运动模型有: x^(t)=cosθ(t)y^(t)=sinθ(t)θ^(t)=1R 其中R为转弯半径。

求解起点到终点的最短路径问题就是对下面问题的一个最优化过程: minLoss(p)=0Tx^(t)2+y^(t)2dt 由于x^(t)2+y^(t)2=1, 即最短路径只与T有关。以RSR控制为例,就变成了3组控制:(-最大转向角,T1)、(0,T2)、(-最大转向角,T3),实际实现时采取迭代式: xt+1=xt+δcos(θ)yt+1=yt+δsin(θ)θt+1=θt+δR 其中δ取值一般为:δ=[0.01,0.05],给定弧长AR,有Ti=ARiδ

3、计算CCC

CCC包括RLR和LRL两种,不涉及到切线的计算。

计算方法如下: p1,p2,p3分别为C1,C2,C3的圆心,C3同时与C1C2相切,rmin为最小转向半径, p1p3¯=p2p3¯=2rminp1p2¯=D=(x2x1)2+(y2y1)2V1=p2p1V2=p1p3θ=p2p1p3,有: θ=cos1(D4rmin) 从而有: p3=(x1+2rmincos(θ),y1+2rminsin(θ)) 注意:LRL模式时,θ需要加上atan2(V1),RLR模式时,θ需要减去atan2(V1)。 对V2做缩放,有: V2=V2||V2||rmin 于是有: pt1=p3+V2 pt2的计算完全类似,这样三段弧线就得到了:(x1,y1)pt1pt1pt2pt2(x2,y2)

非线性最小二乘法控制

车辆的连续运动可以通过离散的多次运动做估计,假设有一个很小的时间间隔,车辆运动模型可以变为: p1(t+h)p1(t)+hv(t)cosθ(t)p2(t+h)p2(t)+hv(t)sinθ(t)θ(t+h)θ(t)+hv(t)Ltanϕ(t)

定义输入向量:uk=(v(kh),ϕ(kh)),定义状态向量:xk=(p1(kh),p2(kh),θ(kh)),则运动模型可以表示为: xk+1=f(xk,uk)=[(xk)1+h(uk)1cos((xk)3)(xk)2+h(uk)1sin((xk)3)(xk)3+h(uk)1tan((uk)2)/L]

其中(uk)1代表向量uk的第一个分量,其他的含义类似。于是车辆运动模型就变成了给定初始点、初始车身朝向和停止点,通过一定时间间隔变化产生一个输入序列,从而模拟车辆运行轨迹的过程,由于要求运动的最短路径,所以它是一个非线性约束最优化问题,希望每次迭代以最小的转向角和速度为代价,保持平稳的运动,最终得到代价最小的最优路径,数学描述如下: mink=1N||uk||2+γk=1N1||uk+1uk||2s.t.x2=f(0,u1)xk+1=f(xk,uk),k=2,...,N1xfinal=f(xN,uN) 其中:N为序列长度,u1,...,uN;x1,...,xN为变量。

参考文献

1、《A Comprehensive, Step-by-Step Tutorial on Computing Dubin’s Curves

2、《Constrained nonlinear least squares》

欢迎关注我的其它发布渠道

表情 | 预览
1 评论
Anonymous Chrome 105.0.0.0 Windows 10.0
2022-10-01 回复

老乡,想请问下咱们这个是什么博客系统啊,WordPress吗?

Anonymous Chrome 105.0.0.0 Mac OS 10.15.7
2022-10-16 回复

@Anonymous , 用hexo在腾讯云上搭建,很方便。

Kris Chrome 125.0.0.0 Windows 10.0
2024-05-25 回复

@Anonymous , 妙啊,学到了!😄

Powered By Valine
v1.3.10