基于改进移动最小二乘法的数据拟合*

王红平，史 明

(长春理工大学机电工程学院机电一体化实验室，长春 130022)

0 引言

最小二乘(LS)方法在数据拟合中的应用最为广泛，其常用的基函数是多项式[1]、有理函数[2]、曲线拟合中的高斯样条、指数样条、平滑样条、B样条、非均匀有理B样条、Bezier曲面和曲面构造中的径向基函数。同时，LS的变形为递归最小二乘(RLS)、总最小二乘(TLS)、偏最小二乘(PLS)、加权最小二乘 (WLS)、广义最小二乘(GLS)[3]、等式约束最小二乘(LSE)。然而，所有基于LS的方法全部是全局逼近格式，不适用于数据量大、分布不规则和分布分散的情况。因此，在测量数据处理中，提出了一种局部逼近的移动最小二乘(MLS)方法。

Shepard在最低阶情形下引入移动最小二乘(MLS)作为逼近方法[4]，Lancaster和Salkauskas将其推广到更高的程度[5]。目前，移动最小二乘(MLS)逼近应用越来越广泛,已用于几何模型的生成[6]、粒子电场的重构[7]和正则化回归中的学习性能[8]。MLS近似在文献中已有大量记载，并被许多学者用于工程实际问题的优化，如:Fu Fangyan等采用最小二乘粒子流体力学(MLSPH)方法求解Burgers方程[9]，Hosseininia M等提出一种基于移动最小二乘(MLS)形状函数无网格方法提高了求解模型的精度[10]。MLS逼近不只应用于工程问题，在图像处理中也得到了大量应用。Hwang Y等提出一种具有空间约束的概率移动最小二乘法，用于图像之间的非线性颜色转移[11]。Yu Chong等基于运动积分最小二乘法进行轮廓图像的变形[12]，在图形处理器上利用移动最小二乘法加速多维插值等[13]。

虽然MLS得到了大量的应用但仍有很大的提升空间，Bayona Victor进行了移动最小二乘法(MLS)与RBF+poly法两种无网格法对插值与导数逼近的比较，证明MLS在拟合高阶多项式时的局限[14]。Zhang L等针对移动最小二乘未考虑全局变量的缺陷，将基于奇异值分解的参数λ引入局部逼近，提出了一种MLS的优化方法用于生成测量数据的曲线和曲面[15]。但上述方法均未对加权函数进行研究，也未对其中影响参数进行分析。

本文基于MLS通过对加权函数的分析与对比,提出了改进的移动最小二乘法(IMLS)。对移动最小二乘的缺陷进行了改进，提升了移动最小二乘的拟合性能，增强了其实用性。

1 移动最小二乘法

1.1 传统移动最小二乘

在MLS近似下，试验函数可以表示为[16]：

(1)

其中，pi(x)，i=1,2,…,m,是单项基函数，m是基函数中的项数，而ai(x)是基函数的系数。一般而言基函数pT(x)=[p1(x),p2(x),…pm(x)]可以是多项式，切比雪夫多项式，勒让德多项式，三角函数，径向基函数等。

通常，基函数有以下形式：

线性基:

一维:pT=(1,x1)m=2

(2)

二维:pT=(1,x1,x2)m=3

(3)

二次基:

(4)

(5)

由兰卡斯特和萨尔考斯卡斯定义的局部近似为：

(6)

对于函数u(x)的精确局部逼近，必须通过加权最小二乘法使局部逼近uh(x)和函数u(x)之间的差值达到最小化。

定义函数:

(7)

其中，u(xI),I=1,2,…,N,是给定的均匀分布的函数结点。w(x-xI)是一个具有紧凑支持的加权函数，并且xI,I=1,2,…,n,是点x影响域中的离散点，下标“I”意味着加权函数的中心位于xI处。图1显示了MLS中加权方案的示意图[17],权值施加在拟合值与给定值之间的平方误差上,通过调节节点的权重值达到拟合的最优化。

图1 加权方案示意图

式(7)的矩阵形式可以表示为：

(8)

(9)

(10)

(11)

a=(a(x),a(x),...,a(x))T

(12)

对于式(8)关于系数a(x)求偏导数，可以得到以下表达式：

A(x)=pTwp,B(x)=pTw

因此系数a(x)为：

a(x)=A-1(x)B(x)

(13)

uh(x)=PT(x)a(x)=PT(x)A-1(x)B(x)u=Φ(x)u

(14)

这里的Φ(x)被称为MLS中的形状函数。

1.2 移动最小二乘的优化

1.2.1 总体最小二问题

求解线性系统AX=B中设矩阵A的误差矩阵为E，向量b的误差向量为e，则表达式等价于

(A+E)·X=(B+e)

与最小二乘法仅关注测量向量误差e不同，总体最小二乘法综合考虑了系数矩阵E和观测向量e的误差，具有较高的计算精度和可靠性。与递推最小二乘法相比，该方法无需设定初值和迭代运算，计算简单且便于实现，不需要大量的样本数据。总体最小二乘法的基本思想是使误差E和e小化，即令矩阵[E⋮e]的F范数最小化，可采用奇异值分解法来解算。

构造增广矩阵C，并进行奇异值分解为：

C=[AB]=UΣVH

(15)

其中，Σ=diag(σ1,σ2,…,σn+d),设C的奇异值σ1≥σ2≥…≥σn+d，并进行分块

(16)

当且仅当V22是非奇异矩阵时且σn≠σn+1，TLS问题有解，TLS解为：

(17)

其中，V12为k阶矩阵，k为待求解维数。

1.2.2 移动最小二乘优化

到所有变量的误差，将在A和B之间设置参数λ的TLS方法引入移动最小二乘法的局部逼近中，对移动最小二乘(MLS)法进行优化[16]。

在移动最小二乘法的优化中，为x处的局部近似值定义Cx：

Cx:=Wx[DAxTBx]=UxΣxVxT

(18)

其中，

Wx=diag(w(x-xI),…,w(x-xn))
D=diag(λ,…,λ)n×n、
T=diag(1-λ,…,1-λ)n×n

在x的影响域中，局部系数近似值由以下关系式解出

a(x)=[a1(x),a2(x),…,am(x)]=-Vx12Vx22

(19)

其中，Vx12为P阶矩阵:P=m。

根据文献[16]所述方法确定参数λ：为了节省计算量和防止导致病态方程，经大量数值实验将λ定为定值λ=0.5。

2 改进的移动最小二乘

在MLS方法中，选择拟合性能较好的高斯函数作为加权函数，但对加权函数没有深入研究。本节对权函数进行了详细的分析，通过对正态权函数和两种典型权函数(高斯权函数和四次样条函数)的数值例子比较，证明了正态权函数在加权方面的优越性，并着重研究了正态权函数及其参数，实现了MLS方法的改进。

2.1 加权函数分析

权函数的选取在MLS 近似中具有重要的作用，对计算结果影响很大，它的选取一般应该满足 4 个条件：非负性、紧支性、单调递减性，光滑性。目前常见的权函数有：高斯型、指数型、样条型以及径向基函数等。在第3节中，当选择特定的基函数时，基函数矩阵A(x)和观测值矩阵B(x)就被确定了，所以它们都与未知点无关，它们只是形成形状函数的常数。因此，只有矩阵w定义在x上，并且形状函数主要由加权函数决定。

加权函数的基本要求是紧支撑、非负定的、连续的、具有较高导数的，以保证系数的唯一性。紧凑的支撑特性是MLS的本质，而加权函数的紧凑性是由支撑半径dmi决定的，从图1可以明显看出，相对较大的支持域意味着在计算中涉及更多的节点，各节点的权函数都是紧支函数，因而在某计算点x处的MLS近似只需利用对x有影响的节点xI(I=1,2,…,N)(节点xI的权函数在x点处不等于零或者说该节点的支撑域半径大于它到x点之间的距离)来构造。节点xI的支撑域一般为半径dmi的圆形区域(二维问题)如图2所示或球形区域(三维问题)。如何选择合适的支承半径取决于拟合误差、平滑度和问题本身的特点。

图2 二维平面支撑域示意图

2.2. 加权函数的选择

(1)正态加权函数:

(20)

σ为参数。

其中，r=|x-xI|/dmi是相对距离，dmi是影响域半径，σ是形状参数，根据影响域与支撑域的性质及加权函数仅在影响域中定义，所以加权函数是一个紧凑的支持函数。正态加权函数如图3所示。此外，拟合的平方误差被加权，只是作为一个移动窗口，简单来说就是对每个样本点计算一次加权最小二乘法，然后对该样本的自变量xi求函数值f(xi)，算出来的(xi,f(xi))就是平滑的结果，所以MLS与IMLS可以看做是WLS与SLS的结合。此外|x-xI|为想要求自变量x到附近样本自变量的欧几里德距离:

图3 正态加权函数图

(2)高斯加权函数为:

(21)

β为参数。

在高斯加权函数中，β是形状参数，r=|x-xI|/dmi是相对距离，dmi是影响域半径。

(3)四次样条函数:

(22)

可以看出在正态加权函数与高斯加权函数中有两个调节参数，分别为dmi,σ和dmi，β，而在四次样条加权函数中只有dmi一个调节参数，图4显示了高斯样条函数和四次样条函数，数值试验表明，β=1.9时，高斯函数与四次样条相似甚至近乎一致。相应地，正态加权函数和高斯加权函数对于IMLS是非常灵活的，都可对形状参数β和σ进行调节从而调节局部逼近效果，但哪一个更适用于IMLS，我们将在下面的具体数值算列中进行分析。

图4 四次样条与高斯加权函数

2.3 数值算例与分析

在这部分给出了两个例子来分析权函数的选择。

示例1：选择函数

y1=sin(x)·cos(3x)

(23)

利用该函数选择均匀分布的节点(xi,yi),i=1,2,…,n。将正态分布的均值为零的随机误差σx和σy分别加到xi和yi中，形成一组测量点,基函数pT(x)=[1,x,x2]。IMLS方法中总体最小二乘参数λ由上述方法确定。采用正态加权函数和高斯加权函数两种加权方法对被测点进行拟合，根据图5所示,显示正态加权函数在拟合性能上优于高斯加权函数，但还需进一步验证。

图5 函数y1=sin(x)·cos(3x)的两种加权的IMLS拟合

示例2:选择函数

(24)

测量点按与示例1相同的方式生成，基函数pT=[1,x]分别用正态加权和高斯加权进行拟合，拟合效果如图6所示。

图6 函数的两种加权的IMLS拟合

从数值实例中可以看出，正态加权拟合性能更好一些。所以在以下的IMLS拟合数值例子中均采用正态加权函数进行加权。

因此，我们可以在IMLS中得出以下结论：

(1)对于基函数的选择，线性基、二次基或高阶多项式都可以作为候选。虽然随着多项式阶数的增加，得到了较好的光滑拟合，但是计算成本会大幅度增加，甚至会导致病态方程。因此，在二维或三维情况下，只选择低阶多项式。

(2)形状函数主要由加权函数决定，对于正态加权函数中的参数设置，影响半径dmi是关键问题。影响半径越大，支撑域越大，拟合光滑性越好，但是计算量越大。减小支撑域，局部性增强，但平滑度下降。从图3可以看出，正态加权函数中的形状参数σ越小，跟踪快速变化能力越强，使局部性拟合效果增强，但平滑度下降，如何选择合适的dmi与σ，取决于拟合误差、平滑度和问题本身特点。

3 IMLS数值拟合实例

通过三个数值算例对IMLS的拟合性能进行了研究，前两个是周期函数，另一个是三维图形函数，测量点按上一节的方式生成，在实例中还应用了MLS方法进行了比较。他们的公式是:

(25)

(26)

z=x*e-x2-y2

(27)

参数设置：

(1)y1,y2中基函数pT(x)=[1,x,x2],m=3。z中基函数pT(x)=[1,x,y,x2,xy,y2],m=6。

(2)函数y1,y2中步长h=0.01,z中步长h=0.08。

(3)λ=0.5。

从数值图像中可以看出IMLS的拟合效果要优于MLS。为了更好的表达拟合优化的效果，列出了误差表1、2和3，误差用均方根误差(RMSE)表示。

(28)

实例1的拟合曲线如图7所示，实例2的拟合曲线如图8所示，实例3的拟合曲面如图9所示。

图7 函数y1=sin(2x)+e-x/3(1+cos(3x))的MLS与IMLS拟合

表1 函数y1=sin(2x)+e-x/3(1+cos(3x))两种拟合方法的RMSE

图8 函数的MLS与IMLS拟合

表2 函数两种拟合方法的RMSE

图9 函数z=x*e-x2-y2的IMLS拟合

表3 函数z=x*e-x2-y2两种拟合方法的RMSE

从上述例子可以看出，IMLS的拟合结果要优于MLS方法。图10显示了上述三个例子IMLS方法的RMSE随影响参数线性比的变化情况。从图中可以看出在一定范围内当dmi/σ=0.9～1.1时，误差值最小。所以在合理的范围内在选取d与σ时，应尽量使dmi/σ接近于0.9～1.1，使拟合效果达到最优化，同时节省试算时间。

图10 均方根误差RMSE随dmi/σ的变化走势图

4 结论

MLS优化考虑了所有变量的误差,将基于奇异值分解的参数λ引入到局部逼近中，但MLS优化没有对加权函数进行深入研究，因此，本文提出了一种MLS近似的改进方法。数值算例和实测数据拟合表明了IMLS优越的拟合性能，可以得出以下结论：

(1)数值算例表明在相同条件下，正态加权的拟合灵活性不输于高斯加权，且正态加权的拟合性能优于高斯加权，所以选择正态加权作为IMLS的加权函数是一种合理的优化方式。

(2)正态加权函数中的影响域半径dmi越大，支撑域越大，拟合光滑性越好，但是计算量越大，减小支撑域，局部性增强，但平滑度下降。形状参数σ越小跟踪快速变化能力越强，使局部性拟合效果增强，但平滑度下降。所以与IMLS方法相比MLS方法具有更大的灵活性。

(3)相同条件下IMLS法与MLS法对离散数据的曲线和曲面拟合结果表明，IMLS方法比MLS方法具有更好的拟合性能，验证了本文提出的IMLS方法的有效性。因此，IMLS是一种有效的数据处理方法。