基于广义逆矩阵理论的分段线性拟合方法2

1年前浏览1606

简介

在《基于广义逆矩阵理论的分段线性拟合方法》一文中，我们基于广义逆矩阵理论，提出了一种分段线性拟合方法。在已知待拟合数据点的情况下，仅需指定分段断点数量以及位置，即可自动计算分段线性拟合直线的系数，并保证相邻拟合直线在分段断点处的连续性。

然而，通过数值试验我们发现，虽然分段线性拟合比单一直线拟合在某些情况下能够取得更好的逼近效果，但分段断点的布置情况同样对拟合效果有很大的影响。鉴于在目前的方法中，对于分段断点的布置仍然具有一定随意性，因此我们希望在指定分段断点数量的情况下，由程序自动给出最优的分段断点位置，以取得最佳的拟合效果。在前文中我们已经提到，这是一个以断点位置作为优化变量，分段拟合效果作为优化目标的典型优化问题。

分段断点位置对拟合效果的影响

在讨论如何求解该优化问题之前，我们首先通过一个简单的示例来研究分段断点位置对拟合效果的影响。这里，我们选取均方根误差（Root Mean Square Error, RMSE）作为评价拟合效果的指标：

式中：f_i为拟合值，y_i为观测值（待拟合数据点），n为观测值的数量。

采用均方根误差作为评价指标的原因为：

（1）均方根误差能够衡量观测值与真值之间的偏差；

（2）它的意义在于开根号后，误差的结果就与数据是同一个级别，可以更好地描述数据；

（3）均方根误差对于一组测量测量中的特大值或特小值反映非常敏感（易受异常值的影响），因此能够更好地反映拟合精度。

此时，该优化问题变为寻找一组特定范围内的分段断点坐标值，使得分段线性拟合后的均方根误差取得最小值，即有：

式中：s₁,…,s_M为分段断点的坐标，M为分段断点数量，x₁,…,x_N为待拟合数据点的横坐标，N为待拟合数据点数量。

为了构造一组确定的待拟合数据点，本文通过如下所示的表达式生成待拟合数据点：

利用该表达式生成的待拟合数据点如图1所示。

图1 待拟合数据点示例

首先考虑最简单的一种情况，即仅有一个分段断点，此时分段断点的取值范围为(0, 2π)，通过计算可以得到拟合均方根误差随分段断点位置的变化规律如下图2所示。

图2 均方根误差随断点位置变化规律（单个断点）

从图中可以看到，当仅有一个断点时，均方根误差随着断点位置的改变呈现反复波动的变化规律，这可能与待拟合数据点本身的特性有关，因为本文采用了正弦函数来构造待拟合数据点。当断点取值约为5时，可以看到拟合均方根误差取得最小值，约为4.0，表明此时可以取得最优的拟合效果。如果继续增加断点数量，同样可以得到当布置两个分段断点时，均方根误差随分段断点的变化规律如下图所示。

图3 均方根误差随断点位置变化规律（两个断点）

继续增加分段断点不便于给出均方根误差的图形化展示，因此本文不再给出其余情况的结果。总结图2和图3中的规律可以发现，均方根误差随断点位置的变化应该是连续的，但在某些特殊位置（使得分段断点数量退化的位置），均方根误差可能会出现突变，因此我们可以猜测均方根误差是一个连续但局部不光滑的函数，并且在函数图像中存在较多的波谷，随着分段断点数量的增加，这些波谷的数量也会随之增加。这会增加优化过程的难度，因为对于采用梯度下降方法之类的优化算法，这会使得优化结果对于断点初值的选取非常敏感，不恰当的断点初值极易使得均方根误差最终落入局部最优解。

基于粒子群算法的分段断点位置优化

为了优化分段断点位置，使得分段线性拟合的均方根误差取得最小值。本文采用了目前应用较为广泛的粒子群优化（Particle Swarm Optimization，PSO）算法。该算法与遗传算法类似，均属于基于生物学原理的算法。该算法最初是受到飞鸟集群活动的规律性启发，进而利用群体智能建立的一个简化模型。粒子群算法在对动物集群活动行为观察基础上，利用群体中个体对信息的共享使整个群体的运动在问题求解空间中产生无序到有序的演化过程，从而获得最优解。与众多优化算法相比，粒子群算法需要调整的参数较少，算法结构简单，并且具有更高效的并行计算能力，可以在较短时间搜寻到全局最优解。并且，粒子群算法对种群数量的变化不是十分敏感，寻优性能稳定，因此算法稳健性较强。当然，目前所有的优化算法均不能保证一定能够找到全局最优解。因此在优化完成后仍然需要对优化结果进行判断。

基本理论

下面，本文对粒子群算法的基本原理做简要介绍。假设在一个D维的目标搜索空间（D个待优化变量）中，有N个粒子组成一个群落，其中第i个粒子表示为一个D维的向量：

第i个粒子的“飞行”速度也是一个D维的向量，记为：

在第t代的第i个粒子向第t+1代进化时，根据如下式子更新：

式中：w为惯性权重，表示群体引起的惯性效应；c₁表示个体学习系数，p_ij(t)代表个体的极值；c₂表示全局学习系数，p_gj(t)代表全体的极值。

本质上说，目前的优化算法都是在满足约束条件的情况下，随机生成大量的优化变量来搜索全局最优解，但为了保证每次迭代时优化变量的取值均能够使得目标函数值尽可能地更趋近于全局最优解，则需要合理地制定搜索策略，而不是盲目地进行搜索。而以粒子群为代表的仿生优化算法则主要参考了自然界中群体动物趋利避害的特点来制定搜索策略。对于上述公式，我们可以以鸟群寻找食物的过程来进行解释。例如，假定鸟群中某一只鸟的当前位置为C，这只鸟认为可能存在食物的位置（个体极值位置）为B，而整个鸟群认为最有可能存在食物的位置（可能的全局最优位置）为A，则这只鸟下一秒的飞行方向如下图所示。

图4 粒子群算法示意图

图中灰色区域为飞鸟的活动范围，代表了优化变量的约束范围，B对应了位于C点的粒子i（飞鸟i）发现的个体极值p_i，A为全局极值p_g。黄色向量为飞鸟当前的飞行方向，绿色向量为飞向个体极值的飞行方向。受到惯性的作用，飞鸟i有保持原有飞行方向的趋势，也有飞向个体极值的自我意识，因此最终的飞行速度方向为黑色向量指向的方向。而蓝色向量为飞向当前全局最优的飞向方向，受到群体意识的影响，飞鸟i也有趋向群体规律的趋势，因此最终的运动方向应该是这三个速度向量的合成，即图中的红色向量。可以看到，如果飞鸟朝红色向量指示的方向飞行，则将更接近于全局最优解。

上面各式中的参数w、c₁和c₂即为控制每个速度向量大小的参数，其取值越大，则飞鸟向该方向运动的趋势也越大。这种搜索策略既考虑了飞鸟的个体意识，也同时考虑了飞鸟群的群体意识，使得在搜索全局最优解的过程中，既能够使得所有粒子逐步趋于全局最优解，也能够各自探索更大的范围，增加寻找到全局最优解的概率。

基于粒子群算法进行寻优的过程可以概括如下：

（1）初始化粒子群，设置粒子群规模、惯性权重、学习系数以及最大迭代次数，并随机生成优化变量初始值、粒子群更新速度；

（2）根据目标函数，计算每个粒子的适度值，寻找个体极值；

（3）粒子间共享个体极值，找到当前粒子群的全局最优解；

（4）每个粒子根据当前速度、位置和全局最优值调整下一步的速度和位置；

（5）反复迭代寻找全局最优值，直至达到最大迭代次数或两次迭代间全局最优值的相对变化小于设定的容差，寻优过程结束。

MATLAB全局优化工具箱

由于MATLAB的全局优化工具箱（Global Optimization Toolbox）已经内置有粒子群优化算法，因此本文通过调用MATLAB的全局优化工具箱来进行分段断点的位置优化。

在MATLAB的全局优化工具箱中，粒子群优化算法的使用语法如下所示。

[x,fval,exitflag,output]=particleswarm(fun,nvars)
[x,fval,exitflag,output]=particleswarm(fun,nvars,lb,ub)
[x,fval,exitflag,output]=particleswarm(fun,nvars,lb,ub,options)

其中，输入参数包括：

fun

目标函数（objective function），以函数句柄（function handle）或函数名的形式来指定。注意，目标函数只能接受一个行矢量并返回一个标量，该行矢量通常由若干个设计变量构成，因此其维度为设计变量的数量（nvars）。

nvars

变量的数量，以正整数的形式来指定，求解器将会把维度为nvars的一个矢量传递给目标函数fun。

lb

设计变量的下界（Lower bounds），以列向量的形式指定。

ub

设计变量的上界（Upper bounds），以列向量的形式指定。通过lb和ub可以使得设计变量满足约束条件lb≤x≤ub。

option

用于更改粒子群优化算法以及优化过程中的一些默认参数，由于参数众多，因此这里只列出一些较为实用的参数，有关这些选项更为详细的介绍可以参见MATLAB的帮助文档。

而输出参数包括：

x

输出满足约束条件，并且使得目标函数取得极小值的最优解。

fval

设计变量取值为x时目标函数的取值。

exitflag

当程序满足终止条件后，程序将会终止，并输出标志符exitflag指明程序停止的原因，下表同样仅给出了exitflag的一些常见取值以及对应的程序终止原因。

有关particleswarm函数的更加详细的介绍可以参考MATLAB的帮助文档。

分段断点位置优化程序的实现

为了调用particleswarm函数，我们需要将前文给出的分段线性拟合程序封装为目标函数，使其接受以行向量形式给出的分段断点位置，并以一个标量的形式输出拟合的均方根误差值。然而，需要注意的是，在分段线性拟合时同样还需要将待拟合数据点传递至目标函数中，由于目标函数不接受其他参数的输入，因此需要将这些额外的变量以其他方式传递到目标函数中。这里，本文提供两种传递方法：全局变量法和匿名函数法。

全局变量法即将待拟合数据点以全局变量的形式传递到目标函数中。我们需要创建一个主程序文件和目标函数文件。在主程序文件中，我们将待拟合数据点定义为全局变量，并调用particleswarm函数进行优化计算：

global Data_x Data_y
[Data_s]= particleswarm(@pwlf_objfunc,nvars,ib,ub,options)

因此，我们在目标函数文件中同样需要将待拟合数据点引用为全局变量：

function RMSE=pwlf_objfunc(Data_s)
global Data_x Data_y
%这里输入定义RMSE的代码
end

全局变量法使用起来较为繁琐，并且一旦定义了全局变量，则该变量名就不能在其他函数中被重复使用。下面，本文介绍一种更为简单的参数传递方法，即匿名函数（Anonymous Functions）法。首先，我们定义一个包含有多个参数的目标函数：

function RMSE= pwlf_objfunc(Data_s,Data_x,Data_y)
%这里输入定义RMSE的代码
end

然后，在主程序文件中，我们为待拟合数据点赋值，并以匿名函数的形式定义一个函数句柄func_obj，并调用particleswarm函数进行优化计算：

func_obj=@(x) pwlf_objfunc(x,Data_x,Data_y)
[Data_s]= particleswarm(@func_obj,nvars,ib,ub,options)

需要注意的是，这些附加参数是在匿名函数被创建时被传递进入匿名函数的，因此如果在创建匿名函数后再次更改这些附加函数的取值，则在匿名函数中这些附加参数的取值并不会被更新。

除此之外，还需要注意的一点是，由于粒子群优化算法会生成大量的随机优化变量，即分段断点位置，则某些分段断点位置不可避免地将会落在边界区域或出现重合的情况，落在边界区域将导致无法正常计算广义逆矩阵，而出现重合的情况可以视为分段断点数量退化。因此需要在定义目标函数的函数体中考虑这些特殊情况，对于传递到目标函数中的分段断点位置进行预处理，以剔除落在边界区域的断点。具体的处理方法可以参见本文附录给出的源代码，这里不再赘述。

计算示例

为了验证本文提出的分段线性拟合方法的稳健性，本文仍然采用如下所示的表达式生成待拟合数据点：

这里将x的范围增大为[0,5π]，这将使得待拟合数据点具有更多的波峰和波谷，导致优化过程很容易落入局部最优解，不易取得较好的拟合效果。

本文取分段断点数量为5个，可以得到分段拟合直线和优化迭代过程中全局最优值（RMSE）的变化如图5所示。

图5 分段拟合直线以及全局最优值随迭代次数变化（5个断点）

可以看到，随着迭代次数的增加，拟合的均方根误差迅速下降。虽然粒子群不一定收敛于全局最优解，但此时已经取得了非常好的拟合效果。此外，由于无需确定分段断点的位置，因此只需对分段断点数量做收敛性分析即可判断取得最佳拟合效果的分段断点数量。例如，对于本例，图6给出了上图中最佳均方根误差随断点数量的变化。

图6 最佳均方根误差随分段断点数量的变化

从图中可以看到，当分段断点数量超过5个后，最佳均方根误差的变化已经趋于稳定，这说明继续增加分段断点数量对于拟合效果的提升已经非常微弱，反而还会增加优化迭代时间。因此可以认为分段断点数量为5时可以取得最佳的拟合效果。事实上，通过观察待拟合数据点的波峰和波谷数量已经可以直接确定分段断点数量应取为5个。

至此，本文已经改进了前文中提出的分段线性拟合方法，去除了需要预先指定断点位置的前提，仅需指定分段断点数量，即可自动完成最优分段断点位置的计算以及分段线性拟合。

附录

分段线性拟合主程序

Piece_wise_linear_fitting.m

clear,clc
%**************************************************************************
%程序用于分段线性拟合 
%输入参数为待拟合数据以及分段断点数量
%注意 待拟合数据点必须以列向量的形式输入
%**************************************************************************
%***************************待拟合数据及参数输入***************************
Data_x=linspace(0,5*pi,101)';                  %待拟合数据点x
Data_y=sin(Data_x)+0.2*sin(5*pi*Data_x);       %待拟合数据点y     
Split_pt_num=5;                                %分段断点数量
%**************************************************************************
%**************************分段线性拟合计算********************************
%创建目标函数匿名函数句柄用于优化计算
plot_status="false";
func_obj=@(x) pwlf_objfunc(x,Data_x,Data_y,plot_status);
%基于粒子群算法优化断点布置位置
lb=zeros(Split_pt_num,1)+min(Data_x);          %变量的下限
ub=zeros(Split_pt_num,1)+max(Data_x);          %变量的上限
options=optimoptions("particleswarm","PlotFcn","pswplotbestf");
[Data_s,fval,exitflag,output]=...
    particleswarm(func_obj,Split_pt_num,lb,ub,options);
%**************************************************************************
%************************绘制分段线性拟合结果******************************
%创建目标函数匿名函数句柄用于输出分段线性拟合结果
plot_status="true";
[RMSE]=pwlf_objfunc(Data_s,Data_x,Data_y,plot_status);
%**************************************************************************

目标函数程序pwlf_objfunc.m

function[RMSE]=pwlf_objfunc(Data_s,Data_x,Data_y,plot_status)
%**************************************************************************
%Data_s为拟合区段断点 以行向量形式指定
%Data_x和Data_y为待拟合数据点 以列向量形式指定
%**************************************************************************
%***************************待拟合数据及参数输入***************************
Data_num=size(Data_x,1);                        %待拟合数据点总数
%剔除重复断点并按升序排列
Data_s=unique(Data_s)';
%剔除位于拟合点边界的断点(只有可能为断点的最小值和最大值)
if Data_s(1)==Data_x(1)
if Data_s(end)==Data_x(end)
        Interval_num=size(Data_s,1)-1;
        Data_s=Data_s(2:end-1);
else
        Interval_num=size(Data_s,1);
        Data_s=Data_s(2:end);
end
else
if Data_s(end)==Data_x(end)
        Interval_num=size(Data_s,1);
        Data_s=Data_s(1:end-1);
else
        Interval_num=size(Data_s,1)+1;
end
end
%判断剔除操作完成后是否出现无断点的情况
if Interval_num==1
    RMSE=Inf;         %若无断点 假设目标函数取值为无穷大
return;
else
    Data_s=[Data_s;max(Data_x)];
end
%**************************************************************************
%********************确定待拟合数据点所在区段的索引值**********************
Interval_index=zeros(Data_num,1);
Start_index=1;
fori=1:Data_num
forj=Start_index:Interval_num
%区段索引为1的特殊情况
ifj==1
if Data_x(i)<=Data_s(j)
                Interval_index(i)=j;
break;
end
continue;
end
%区段索引不为1的一般情况
if (Data_x(i)>Data_s(j-1)&&Data_x(i)<=Data_s(j))
            Interval_index(i)=j;
            Start_index=j;
break;
end
end
end
%**************************************************************************
%***************************确定等价关系矩阵A******************************
%初始化矩阵A
Matrix_A=zeros(Data_num,2*Interval_num);
%计算矩阵A的分量
fori=1:Data_num                   %行循环
    Matrix_A(i,2*(Interval_index(i)-1)+1)=Data_x(i);
    Matrix_A(i,2*(Interval_index(i)-1)+2)=1;
end
%**************************************************************************
%***************************生成连续性条件矩阵*****************************
%初始化存放不同区段的连续性条件矩阵的元胞数组
Matrix_C=cell(Interval_num-1,1);
%循环创建不同区段的连续性条件矩阵
forj=1:Interval_num-1
    Matrix_C{j,1}=[-1/Data_s(j) 11/Data_s(j)];
end
%**************************************************************************
%**************************组装转换矩阵分矩阵******************************
%初始化存放不同区段的转换矩阵分矩阵的元胞数组
Matrix_T=cell(Interval_num-1,1);
%考虑i=1的特殊情况[T1]
Matrix_T{1}=zeros(2*Interval_num,2*Interval_num-1);
Matrix_T{1}(1,1:3)=Matrix_C{1};               %填充连续性条件矩阵[C]
Matrix_T{1}(2:end,:)=eye(2*Interval_num-1);   %填充单位矩阵
%考虑i>1的一般情况[Ti]
forj=2:Interval_num-1
    Matrix_T{j}=zeros(2*Interval_num-j+1,2*Interval_num-j);
    Matrix_T{j}(1:j-1,1:j-1)=eye(j-1);          %填充上侧单位矩阵
    Matrix_T{j}(j,j:j+2)=Matrix_C{j};           %填充连续性条件矩阵
    Matrix_T{j}(j+1:end,j:end)=eye(2*Interval_num-2*j+1);   %填充下侧单位矩阵
end
%**************************************************************************
%****************************计算分段拟合直线系数**************************
%计算转换矩阵[Pai_T]
Matrix_Pai_T=Matrix_T{1};
forj=2:Interval_num-1
    Matrix_Pai_T=Matrix_Pai_T*Matrix_T{j};     %通过连乘分矩阵获得转换矩阵
end
%计算独立变量系数计算矩阵[E]
Matrix_E=Matrix_A*Matrix_Pai_T;
%计算独立变量系数计算矩阵[E]的广义逆矩阵
Matrix_E_inv=pinv(Matrix_E);
%计算分段样条曲线的系数矩阵
Matrix_Coeff=Matrix_Pai_T*Matrix_E_inv*Data_y;
%**************************************************************************
%************************计算分段线性拟合的均方根误差**********************
Fitting_y=zeros(Data_num,1);
fori=1:Data_num
    Fitting_y(i)=Matrix_Coeff(2*(Interval_index(i)-1)+1)*Data_x(i)+...
        Matrix_Coeff(2*(Interval_index(i)-1)+2);
end
RMSE=sqrt((sum((Data_y-Fitting_y).^2))/Data_num);
%**************************************************************************
%*************************判断是否需要绘制函数图像*************************
if plot_status=="true"
%初始化用于绘制分段拟合直线的元胞数组
    Fitting_x=cell(Interval_num,1);
    Fitting_y=cell(Interval_num,1);
fori=1:Interval_num
ifi==1
            Fitting_x{i}=linspace(min(Data_x),Data_s(i),51)';
else
            Fitting_x{i}=linspace(Data_s(i-1),Data_s(i),51)';
end
end
%计算分段样条拟合直线
fori=1:Interval_num
%基于分段样条拟合直线的系数进行计算
        Fitting_y{i}=Matrix_Coeff(2*(i-1)+1)*Fitting_x{i}+...
            Matrix_Coeff(2*(i-1)+2);
end
%绘制分段样条直线
figure;
scatter(Data_x,Data_y,70,"Filled","Black");hold on;
fori=1:Interval_num
plot(Fitting_x{i},Fitting_y{i},"LineWidth",2,"Color","red");
end
    xlabel("x");ylabel("y");title("分段线性拟合结果");hold off;grid on;
%输出分段样条直线系数
    fprintf("分段样条直线为:y=ax+b\n");
fori=1:Interval_num
        fprintf("第%d条拟合直线系数为:",i);
        fprintf("a=%10.5f,",Matrix_Coeff(2*(i-1)+1));
        fprintf("b=%10.5f,",Matrix_Coeff(2*(i-1)+2));
ifi==1
            fprintf("取值范围为:[%10.5f,%10.5f]\n",min(Data_x),Data_s(i));
else
            fprintf("取值范围为:[%10.5f,%10.5f]\n",Data_s(i-1),Data_s(i));
end
end
else
return;
end
%**************************************************************************
end