移动最小二乘法MLS（Moving Lest Squares）简要介绍

1. Overview

移动最小二乘法（MLS, Moving Least Squares）是建立大量离散数据拟合曲线的理想方法。当大量离散数据的分布较为杂乱时， 使用传统的最小二乘法，往往需要对数据进行分段拟合，此外还要避免相邻分段上的拟合曲线不连续不平滑的问题。而MLS法在处理相同问题时则不需要上述这些繁琐的步骤，简单易于实现。

在MLS法中，需要在一组不同位置的节点（node）附近建立拟合曲线，每个节点都有自己的一组系数 ($a_j(x_{node})$) 用于定义该位置附近拟合曲线的形态。因此，在计算某个节点附近的拟合曲线时，只需要计算该点的该组系数值($a_j$)即可。

此外，每个节点的系数（$a_j$）取值只考虑其临近采样点，且距离节点越近的采样点贡献越大，对于未置较远的点则不予考虑。

2. 拟合函数

相比于传统的最下二乘法，MLS法中的拟合函数不是一个多项式，而是一组系数向量函数$a_j(x)$和基函数$p_j(x)$，其中$x$为空间坐标。

某个节点node附近的拟合函数即为$u_{node}(x)$，具体定义为如下公式。
$u_{node}(x)={\sum }_{j=1}^m{a}_j(x_{node})\ast p_j(x)$

其中，$x_{node}$为节点的空间坐标，$x$为节点node附近的某一个位置坐标，$a_j$即为用于定义节点node附近拟合曲线的一组系数，$p_j(x)$则是一组基函数，其常用形式如下：

空间维度	一阶形式	二阶形式
1维： $\mathbf{x}=(x)$	$p_1=1,p_2=x$	$p_1=1,p_2=x,p_3=x^2$
2维：$\mathbf{x}=(x,y)$	$p_1=1,p_2=x,p_3=y$	$$\begin{gathered} p_1=1,p_2=x,p_3=y, \\ {p}_4=xy,p_5=x^2,p_6=y^2 \end{gathered}$$
3维：$\mathbf{x}=(x,y,z)$	$p_1=1,p_2=x,p_3=y,p_4=z$	$$\begin{gathered} p_1=1,p_2=x,p_3=y,p_4=z, \\ {p}_5=xy,p_6=xz,p_7=yz, \\ {p}_8=x^2,p_9=y^2,p_{10}=y^2 \end{gathered}$$

MLS法中，需要通过调整系数$a$，使得节点node附近采样点取值与拟合函数在采样点取值之间差的加权平方和最小，由此可建立最优化模型J。
$$\begin{gathered} J=\sum _{p=1}^{n}w(x_{node}-x_p)[u(x_{node})-u_p{]}^2 \\ =\sum _{p=1}^{n}w(x_{node}-x_p)[\sum _{j=1}^m{a}_j(x_{node})\ast p_j(x)-u_p{]}^2 \\ =\sum _{p=1}^{n}w(x_{node}-x_p)[\sum _{j=1}^m{a}_j^{node}\ast p_j(x)-u_p{]}^2 \end{gathered}$$
其中，$x_{node}$为节点node的空间位置，$p$为采样点，$x_p$为采样点p的空间位置，$u_p$是采样点p的取值，此处$a_j(x_{node})=a_j^{node}$。

上式中的$w$为权函数，该函数保证距离node越近的采样点，对于函数$J$的贡献越大，而对于距离node较远的位置的采样点，则不会对J产生影响。函数$w$的常用形式如下（三次样条曲线）。

其中，当采样点p与节点node的距离大于一定值时（此处为2），权重函数$w$的取值为零。(上图中$|x|$为矢量$x$的模)。当然，该函数的邻域阈值可根据需要设定为合适的值。

最终求解出在问题域上全部节点的系数，即可得到全域上的拟合曲线。

3. 系数计算

如上文所述，每个节点都有自己的一组系数 ($a_j(x_{node})$) 用于定义该位置附近拟合曲线的形态，本节将叙述该系数的具体求解方法。

我们已经获得了参数为系数$a_j^{node}$的最优化模型$J$，当$J$的取值最小时，则可获得节点node附近的拟合曲线$u_{node}(x)$。

求解最优化模型最直接的方式，即是对模型求导，因此我们可得到如下方程。

$$\frac{\partial J}{\partial a_i^{node}}=2\sum _{p=1}^{n}w(x_{node}-x_p)[\sum _{j=1}^m{a}_j^{node}\ast p_j(x)-u_p]p_i(x)$$
当上式得零时，$J$取得最小值，因此可得：

$$\sum _{p=1}^{n}w(x_{node}-x_p)[\sum _{j=1}^m{a}_j^{node}\ast p_j(x)-u_p]p_i(x)=0$$

$$\sum _{p=1}^{n}w(x_{node}-x_p)\sum _{j=1}^m{a}_j^{node}\ast p_j(x)=\sum _{p=1}^{n}w(x_{node}-x_p)u_p\ast p_i(x)$$

我们引入下述记号以简化表达（内积）：
$$\begin{gathered} (p_i,p_j)=\sum _{p=1}^{n}w(x_{node}-x_p)\ast p_i(x)\ast p_j(x) \\ (u_p,p_j)=\sum _{p=1}^{n}w(x_{node}-x_p)\ast u_p\ast p_j(x) \end{gathered}$$

因此，最终我们可将$\frac{\partial J}{\partial a_i^{node}}=0$整理成如下的线性方程组。该线性方程组也被称为法方程。
$$\left[\begin{array}{ccccc}(p_1,p_1)& ...& (p_1,p_j)& ...& (p_1,p_m)\\ ...& & & & \\ (p_i,p_1)& ...& (p_i,p_j)& ...& (p_i,p_m)\\ ...& & & & \\ (p_m,p_1)& ...& (p_m,p_j)& ...& (p_m,p_m)\end{array}\right]\ast \left[\begin{array}{c}{a}_1^{node}\\ ...\\ a_j^{node}\\ ...\\ a_m^{node}\end{array}\right]=\left[\begin{array}{c}(u_p,p_1)\\ ...\\ (u_p,p_j)\\ ...\\ (u_p,p_m)\end{array}\right]$$
该线性方程组可写成矩阵形式 ：$M\vec{a}=\vec{b}$，其中$M$为由基函数构成的系数矩阵，$\vec{a}$是未知数向量，也就是我们要求节点node上的那一组系数。因此未知数向量$\vec{a}$也可表示为： $\vec{a}=M^{-1}\vec{b}$。

解出上述线性方程组，即可得到在节点node上的改组系数$a_j^{node}$。通过系数$a_j^{node}$与基向量$p_j(x)$就可建立出节点node附近的拟合函数。

进一步的，我们可以将任意位置上的拟合函数写成如下形式（即节点node的位置可为任意点$\mathbf{x}$）：

$$u(\mathbf{x})=\sum _{j=1}^m{a}_j(\mathbf{x})\ast p_j(\mathbf{x})=P\cdot \vec{a}$$
将 $\vec{a}=M^{-1}\cdot \vec{b}$带入后，可得：
$$u(x)=P\cdot \vec{a}=P\cdot M^{-1}\cdot \vec{b}$$
其中：
$$M=\left[\begin{array}{ccccc}(p_1,p_1)& ...& (p_1,p_j)& ...& (p_1,p_m)\\ ...& & & & \\ (p_i,p_1)& ...& (p_i,p_j)& ...& (p_i,p_m)\\ ...& & & & \\ (p_m,p_1)& ...& (p_m,p_j)& ...& (p_m,p_m)\end{array}\right]$$

$$b=\left[\begin{array}{c}(u_p,p_1)\\ ...\\ (u_p,p_j)\\ ...\\ (u_p,p_m)\end{array}\right]$$

因此我们可以在任意点$\mathbf{x}$上计算出拟合函数及其取值，并且空间上相邻的两点的取值是很平滑的。

4. 注意

①本文中的$p_i(x)$为基向量，而下标$p$则指的是采样点$p$，如$x_p$；
②本文中，采样点的个数表示为$n$，基向量的个数为$m$;
③网格node上的一组系数表示为$a_i^{node}$， 其中$i=1,...m$;
④每个采样点上携带有两部分信息，一部分是该采样点的空间位置$x_p$，另一部分是采样点取值$u_p$，这两部分信息都是已知的。

5. 扩展

MLS法其实是一种广义的最小二乘法方法。MLS法的核心优化模型J：
$$\begin{gathered} J=\sum _{p=1}^{n}w[u(\mathbf{x})-u_p{]}^2 \\ =\sum _{p=1}^{n}w(\mathbf{x}-{\mathbf{x}}_p)[\sum _{j=1}^m{a}_j(\mathbf{x})\ast p_j(\mathbf{x})-u_p{]}^2 \end{gathered}$$
① 当权函数 $w=1$ 时（如下图曲线），则上述MLS方法退化成标准的标准的最小二乘法。

② 当权函数 $w=w(\mathbf{x})$ 为高斯函数或样条曲线时（如下图曲线），则上述MLS方法退化成固定点最小二乘法(FLS, fixed least-squares)

③ 在MLS法中权函数为：$w=w_p(\mathbf{x})=w(\mathbf{x}-{\mathbf{x}}_p)$，一般也为高斯函数或类似样条函数，且以$\mathbf{x}$点为中心点取最大值，在$\mathbf{x}$其相邻区域上取值逐渐减小。

④ 当权函数为$w=w(\mathbf{x})$的曲线为如下形式时，则上述MLS法退化为多固定点最小二乘法（MFLS，Multiple fixed least-squares）

另附我自己写的代码，若有bug请回复。

import numpy as np  
import random
from matplotlib import pyplot as plt

Consistency = 3  #阶数
d = 0.9	#紧支域半径
dotNumber = 36	#散点数量
dotDinstance = 0.05	#散点间距

test = 0	# 0 为曲线，1 为直线

x = []
v = []
Vp = []
p = np.zeros((dotNumber, Consistency))
M = []


def w(xi,xj,d):
    s = abs(xi-xj)/d
    if(s <= 0.5):
        return (2/3)-4*s+4*s**2+4*s**3
    elif(s<=1):
        return (4/3)-4*s+4*s**2-(4/3)*s**3
    else:
        return 0

def b_compute(i, x, v, d):
    b = []
    for c in range(Consistency):
        t = 0
        for j in range(dotNumber):
            # if abs(x[i] - x[j]) < d:
                t +=  p[j][c] * w(x[i], x[j], d) * v[j]
        b.append(float(t))
    return b
        # for j in range(0, Consistency):
        #     w()


def A_compute(i, x, d):
    A = []
    for ci in range(Consistency):
        vec = []
        for cj in range(Consistency):
            t = 0
            for j in range(dotNumber):
                t +=  w(x[i], x[j], d) * p[j][cj] * p[j][ci]
            vec.append(float(t))
        A.append(vec)
    return A



def polySet(p,x):
    for i in range(dotNumber):
        for j in range(0, Consistency):
            p[i][j] = x[i]**(j)


def Vcontruct(a,x,i):
    value = 0
    for c in range(len(a)):
        value +=  a[c] * x[i]** c
    return value

def main():
    x = np.arange(-0.5*dotNumber*dotDinstance, 0.5*dotNumber*dotDinstance, dotDinstance)

    for i in range(0, dotNumber):
        if test == 0:
            v.append(i + random.randint(-50,50)/10)
            # v.append(random.randint(-50, 50) / 10)
        else:
            v.append(i + 0 / 10)

    polySet(p,x)

    for i in range(dotNumber):
        b = b_compute(i, x, v, d)

        A = A_compute(i, x, d)
        a = np.linalg.solve(A, b)
        Vp.append(Vcontruct(a, x, i))
        # V.append(a[0]+a[1]*x[i])

    plt.title("MLS")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.scatter(x,v,color = "green", s= 10)
    plt.plot(x, Vp)
    plt.show()

if __name__ =='__main__':
    main()