spline function

文章目录

• 背景
• Spline Basis Function
• 代码
• 参考
• 福利

背景

最近在研究functional 回归，发现有一些smoothing信号处理方法，跟我以前的一些肤浅的想法居然有一些共性，看来不是想不到，而是不敢想，想得不够深入的问题。这种算法提出已经比较久了，其中比较流行的一种平滑处理算法是基于B-spline。样条插值，作为一种插值或者函数逼近，无论是做图形图像还是数值分析，老早就接触过了，所谓了Spline Basis Function，第一反应无非是将函数切割，分段拟合罢了，保证分割处的n次导数相同即可。深入进行发现，这个基函数概念使得问题描述变得更加抽象，和之前的理解不太一样，下面记述一下自己对它的理解。看了一些不错的资料，但感觉描述的过于详细，反而觉得耽误了那些只想了解基本原理的人。所以，本文只想交代最重要的部分，细节可以参考文末给出的连接。

Spline Basis Function

首先，对于目标曲线进行分段，如何确定最优的分割点，据目前的了解也是相当复杂。这里，作为测试，分割点是自己直接设定的，详情可以参照后面的代码。
假设m+1个分割点，$u_0 \leq u_1 \leq u_2 \leq ... \leq u_m,u_i$称为节点（knots）。这样，原来的曲线分为m段区间。
下面给出Cox-de Boor递归公式，也就是Spline Basis Function的生成公式
$N_{i,0}(u) = \left\{\begin{matrix} 1,& u_i \leq u < u_{i+1} \\ 0, & otherwise \end{matrix}\right.\\ N_{i,p}(u) = \tfrac{u-u_i}{u_{i+p}-u_i} N_{i,p-1}(u) + \tfrac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}N_{i+1,p-1}(u)$
我们定义$N_{i,p}$为第$i$个$p$次$B-$样条基函数,定义有一些拗口，它的递推公式也有点让人一头雾水，没关系，很快就清楚了。
关于它的生成方式，可以引用下图来表示。第一列是分段区间，第二列为0次基函数，依次类推。

我们可以发现两个重要特点

1. 区间的数量决定了基函数的最大次数 。区间数量越多，基函数的最大次数越大，不用说，函数的逼近效果也越好，但也容易过拟合。
2. 基函数次数越大，它的作用范围也越大。如下图部分，$N_{1,3}$它在蓝色虚线部分的区间$[u_1,u_5)$取值为非0，其他全部为0.
   
   总结一下，我们可以得到，对于$N_{i,p}$而言，它的作用范围是$[u_i,u_{i+p+1})$。
   那么可以知道， $N_{i,p-1}$的作用范围$[u_i,u_{i+p})$，$N_{i+1,p-1}$的作用范围$[u_{i+1},u_{i+p+1})$，
   $\color{red}{现在可以解释}$$\color{red}N_{i,p}$$\color{red}{的递推公式}$
   $N_{i,p}(u) = \tfrac{u-u_i}{u_{i+p}-u_i} N_{i,p-1}(u) + \tfrac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}N_{i+1,p-1}(u)$
   $N_{i,p}$实际上就是$N_{i,p-1}$和$N_{i+1,p-1}$的加权和，加权系数是怎么确定的呢，是根据u在两者对应的作用区间的位置决定的。我们可以看到，当$u$在作用范围的边界时候，基函数的取值都接近于0.

下面看一段代码，来看看具体是怎么一回事情。代码将[0,4]均匀分成4个区间。然后，分别生成0，1，2，3次基函数，刚好4张图。

代码

代码比较简单，不过写得烂也没注释，跑一下就清楚了

clear;
linecolor = ['r','g','b','y','k','c','m','k','r','g','b'];
%knot span
knots = [0,1,2,3,4];
Len = 40;    
x = linspace(0,4,Len);
figure('color','w');
for i = 1:1:4
    bf  = GenerateBasicFunctionCurves(i-1,x,knots);
    subplot(2,2,i);
    drawsubplot(i-1,bf,linecolor,knots);
end

function v  = SplineBasisFunction(i,p,u,knots)
   u_i  = knots(i);    
   u_i1 = knots(i+1);
   if p==0
       if u>=u_i && u<u_i1
           v = 1;
           return;
       end
       v = 0;
       return;
   end
   u_ip1 = knots(i+p+1);
   u_ip = knots(i+p);
   v = (u-u_i)/(u_ip-u_i)*SplineBasisFunction(i,p-1,u,knots)+(u_ip1-u)/(u_ip1-u_i1)*SplineBasisFunction(i+1,p-1,u,knots);
end
function bf = GenerateBasicFunctionCurves(p,x,knots)
    for i = 1:1:length(knots)-1-p
        for j = 1:1:length(x)
            bf(i,j) = SplineBasisFunction(i,p,x(j),knots);
        end
    end
end
function drawsubplot(p,bf,linecolor,knots)
    legendtxt = [];
    for i = 1:1:length(knots)-1-p
        plot(bf(i,:),strcat('-o',linecolor(i)),'LineWidth',1.5);
        hold on;
        txt = 'N';
        txt = strcat(txt,num2str(i-1));
        txt = strcat(txt,',');
        txt = strcat(txt,num2str(p));
        legendtxt  =  [legendtxt;string(txt)] ;
    end
    legend(legendtxt);
    plotsetting
    title(p);
end
function plotsetting()
    a=gca;
    a.FontSize=20; 
    set(gca,'fontweight','bold');
    set(gca,'Linewidth',2);
    xticks([1 ,11, 21,31,40]);
    xticklabels({'0','1','2','3','4' });
end

从上图可以看到，随着基函数的次数p增加，基函数的作用范围增大，基函数的数目减少。基函数的最大值一般靠近中间部位，两头趋近为0。上图有趣的地方在于，0次基函数就是一个门函数，1次基函数像是三角滤波器。

上图展现的是分割点均匀分布的情况，下图来一个不均匀分布的

修改上面代码的结点分割部分即可，其他不用调整

knots = [0,0.3,0.8,3.6,4];

参考

1. CS3621 Introduction to Computing with Geometry Notes
2. B-spline Curves 学习之B样条基函数的定义与性质（2）（中文版）