背景
最近在研究functional 回归,发现有一些smoothing信号处理方法,跟我以前的一些肤浅的想法居然有一些共性,看来不是想不到,而是不敢想,想得不够深入的问题。这种算法提出已经比较久了,其中比较流行的一种平滑处理算法是基于B-spline。样条插值,作为一种插值或者函数逼近,无论是做图形图像还是数值分析,老早就接触过了,所谓了Spline Basis Function,第一反应无非是将函数切割,分段拟合罢了,保证分割处的n次导数相同即可。深入进行发现,这个基函数概念使得问题描述变得更加抽象,和之前的理解不太一样,下面记述一下自己对它的理解。看了一些不错的资料,但感觉描述的过于详细,反而觉得耽误了那些只想了解基本原理的人。所以,本文只想交代最重要的部分,细节可以参考文末给出的连接。
Spline Basis Function
首先,对于目标曲线进行分段,如何确定最优的分割点,据目前的了解也是相当复杂。这里,作为测试,分割点是自己直接设定的,详情可以参照后面的代码。
假设m+1个分割点,
u
0
≤
u
1
≤
u
2
≤
.
.
.
≤
u
m
,
u
i
u_0 \leq u_1 \leq u_2 \leq ... \leq u_m,u_i
u0≤u1≤u2≤...≤um,ui称为节点(knots)。这样,原来的曲线分为m段区间。
下面给出Cox-de Boor递归公式,也就是Spline Basis Function的生成公式
N
i
,
0
(
u
)
=
{
1
,
u
i
≤
u
<
u
i
+
1
0
,
o
t
h
e
r
w
i
s
e
N
i
,
p
(
u
)
=
u
−
u
i
u
i
+
p
−
u
i
N
i
,
p
−
1
(
u
)
+
u
i
+
p
+
1
−
u
u
i
+
p
+
1
−
u
i
+
1
N
i
+
1
,
p
−
1
(
u
)
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)
Ni,0(u)={1,0,ui≤u<ui+1otherwiseNi,p(u)=ui+p−uiu−uiNi,p−1(u)+ui+p+1−ui+1ui+p+1−uNi+1,p−1(u)
我们定义
N
i
,
p
N_{i,p}
Ni,p为第
i
i
i个
p
p
p次
B
−
B-
B−样条基函数,定义有一些拗口,它的递推公式也有点让人一头雾水,没关系,很快就清楚了。
关于它的生成方式,可以引用下图来表示。第一列是分段区间,第二列为0次基函数,依次类推。
我们可以发现两个重要特点
- 区间的数量决定了基函数的最大次数 。区间数量越多,基函数的最大次数越大,不用说,函数的逼近效果也越好,但也容易过拟合。
- 基函数次数越大,它的作用范围也越大。如下图部分,
N
1
,
3
N_{1,3}
N1,3它在蓝色虚线部分的区间
[
u
1
,
u
5
)
[u_1,u_5)
[u1,u5)取值为非0,其他全部为0.
总结一下,我们可以得到,对于 N i , p N_{i,p} Ni,p而言,它的作用范围是 [ u i , u i + p + 1 ) [u_i,u_{i+p+1}) [ui,ui+p+1)。
那么可以知道, N i , p − 1 N_{i,p-1} Ni,p−1的作用范围 [ u i , u i + p ) [u_i,u_{i+p}) [ui,ui+p), N i + 1 , p − 1 N_{i+1,p-1} Ni+1,p−1的作用范围 [ u i + 1 , u i + p + 1 ) [u_{i+1},u_{i+p+1}) [ui+1,ui+p+1),
现 在 可 以 解 释 \color{red}{现在可以解释} 现在可以解释 N i , p \color{red}N_{i,p} Ni,p 的 递 推 公 式 \color{red}{的递推公式} 的递推公式
N i , p ( u ) = u − u i u i + p − u i N i , p − 1 ( u ) + u i + p + 1 − u u i + p + 1 − u i + 1 N i + 1 , p − 1 ( u ) 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) Ni,p(u)=ui+p−uiu−uiNi,p−1(u)+ui+p+1−ui+1ui+p+1−uNi+1,p−1(u)
N i , p N_{i,p} Ni,p实际上就是 N i , p − 1 N_{i,p-1} Ni,p−1和 N i + 1 , p − 1 N_{i+1,p-1} Ni+1,p−1的加权和,加权系数是怎么确定的呢,是根据u在两者对应的作用区间的位置决定的。我们可以看到,当 u 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];