2026a

# 使用自然三次样条进行最小二乘逼近

最小二乘逼近的构造通常需要待逼近数据空间的基。而正如“自然”三次样条空间的示例所展示的，这个基的显示构造并不总是那么直接。

本文表明一个显示的基实际上是不需要的，我们其实只需要在逼近空间的某种方式的插值就可以。而这就必需要曲线拟合工具箱样条函数对于向量值样条的支持了。

本文讨论了以下使用“自然”三次样条进行最小二乘插值的相关问题。

问题
通常解法
对基映射的需求
“自然”三次样条的基映射
单行解法
对正确外插的需求
正确的单行解法
使用三次样条的最小二乘插值

# 问题

你想要使用给定断点 b[1]<...<b[l+1] 的“自然” 三次样条组成的空间 S 构造给定数据 (x,y) 的最小二乘逼近。

# 通常解法

如果给定断点序列 b 的所有“自然”三次样条组成的线性空间 S 的一个基 (f1,f2,...,fm)，那么就应该寻求 c[1]f1+c[2]f2+...+c[m]fm 形式的最小二乘逼近，其中 c 是线性系统 A*c=y 的最小二乘解，稀疏矩阵 A 由下式给出：

换句话说，c=A\y。

# 对基映射的需要

通常解法似乎要求你知道一个基。然而，为了构造系数序列 c，你只需要知道矩阵 A。对此，知道一个基映射就足够了。基映射是指函数 F 使得 F(c) 返回由特定加权和 c[1]f1+c[2]f2+...+c[m]fm 指定的样条。因为有了这样的函数，实际上对于 j=1:m，A 的第 j 列就是 fnval(F(ej),x)，其中 ej 是 eye(m) 的第 j 列，而 eye(m) 是阶数为 m 的单位矩阵。

好在曲线拟合工具箱样条函数支持向量值函数，因此你也应该能构造出来向量系数 c[i] 的基映射 F。然而，本工具箱认为的向量值系数指的是列向量，因此序列 c 应该是列向量的行向量，也就是一个矩阵。这样，F(eye(m)) 就表示基 fi 的第 i 个元素，i=1:m。因此，假设数据点向量 x 为行向量，fnval(F(eye(m)),x) 就是一个第 (i,j) 个元素为 fi 在 x[j] 的值的矩阵，也就是你想要的矩阵 A 的转置。另一方面，正如刚刚所指出的那样，你的基映射 F 希望系数序列 c 是个行向量，即向量 A\y 的转置。因此，假设对应的数据值向量 y 是个行向量，那么从 S 得到的数据 (x,y) 的最小二乘逼近就是

F(y/fnval(F(eye(m)),x))

如果你觉得 x,y 可以是任意形式的相同长度的向量，也可以用

F(vec(y)'/fnval(F(eye(m),vec(x)')))

# “自然”三次样条的基映射

那给定断点序列 b[1]<...<b[l+1] 的“自然”三次样条组成的线性空间 S 的基映射 F 是什么呢？假设这个线性空间的维数是 m，则映射 F 应该在 m 长度的向量和 S 之间建立一个线性的一对一的映射关系，而这正是 csape(b,.,"var") 所做的。

考虑以下函数 F:

F(c) = csape(b,c,"var")

对于给定的向量 c（与 b 的长度相同），它提供了*唯一一个*具有断点序列 b，在 b[i] 处值为 c[i]，i=1:l+1 的“自然”三次样条。这里的唯一性是关键，这保证了向量 c 和样条 F(c) 的对应关系是一一的。特别地，m 等于 length(b)。而且，因为函数 f 在一个点 t 的值 f(t) 线性地依赖于 f，这种唯一性还确保了 F(c) 是线性地依赖于 c 的（因为 c 等于 fnval(F(c),b) 而可逆线性映射的逆映射也是线性映射）。

# 单行解法

综上，你就可以得到以下代码

csape(b,vec(y)'/fnval(csape(b,eye(length(b)),"var"),vec(x)',"var")

这就可以用来得到给定断点序列 b 的“自然”三次样条的最小二乘逼近。

# 对合理外插的需要

让我们用一些数据试一下吧，比如人口统计数据，这可以使用以下 Syslab 指令得到

using TyCurveFitting
using TyBase
using TyPlot
include(pkgdir(TyCurveFitting) * "/examples/docs/census.jl")

这提供了 cdate=1790:10:1990 ，表示年份，以及 pop，表示对应的人口数据。使用 1810:40:1970 作为断点序列。

b = 1810:40:1970
s = csape(b,vec(pop)'/fnval(csape(b,eye(length(b)),"var"),vec(cdate)'),"var")
fnplt(s,[1750,2050],2.2)
hold("on")
plot(cdate,pop,"o")
hold("off")

这样就得到了图使用有三个内部断点的“自然”三次样条的最小二乘逼近。这个图使用深蓝色表示得到的逼近，也展示出了给定的数据。

这看上去像是个好的逼近，除了这看上去不像是个“自然”三次样条。一个“自然”三次样条必须在它的第一个断点左侧以及最后一个断点的右侧保持线性，这个插值这两个条件都不满足，这是因为以下的原因。

我们使用 csape 得到的“自然”三次样条插值是 pp 型，该样条的区间来自于其基本区间与给定数据位点区间的拓展。另一方面，在基本区间以外的 pp 型的取值在 Syslab 中是通过数学库的 ppval 或者曲线拟合库的 fnval 完成的，这些是使用 pp 型的相关的端点多项式进行延拓完成的，也即是保阶外插。而在“自然”三次样条的情形下，你实际上想要的是二阶外插。这意味着你想要在第一个端点的左侧看到一个与三次样条在第一个断点处相交并相切的直线。这样的外插是通过 fnxtr 完成的。因为“自然”三次样条在第一个端点处二阶导为零，这样的外插实际上可以说它阶数为三，即它满足三个连续条件。同样地，在三次样条最后一个断点以外，你想要看到一个与样条在最后一个断点处相交并相切的直线，这同样也是 fnxtr 提供的。

使用有三个内部断点的“自然”三次样条的最小二乘逼近

# 正确的一行解法

以下的一行代码提供了使用给定断点序列 b 的“自然”三次样条对数据 (x,y) 的正确最小二乘逼近：

fnxtr(csape(b,vec(y)'/fnval(fnxtr(csape(b,eye(length(b)),"var")),vec(x)'),"var"))

然而实际上这一行有点长了。

接下来的代码使用了上面的正确公式并在之前的图之上使用更细的红线来表示了得到的逼近，这也就是图使用有三个内部断点的“自然”三次样条的最小二乘逼近。

ss = fnxtr(csape(b,vec(pop)'/fnval(fnxtr(csape(b,eye(length(b)),"var")),vec(cdate)'),"var"))
hold("on")
fnplt(ss,[1750,2050],1.2,"r")
grid("on")
hold("off")
legend(["incorrect approximation","population","correct approximation"])

# 三次样条的最小二乘逼近

如果您想要的是给定断点序列 b 的所有三次样条组成的空间 S 的逼近，那当前解法就已经足够了。甚至都可以不用使用曲线拟合工具箱的样条函数，因为可以用数学工具箱的 spline。你知道，如果 c 是一个比 b 多两个元素的序列，spline(b,c) 提供了给定断点序列 b ，对于所有 i 在 b[i] 取值 c[i+1]，在 b[1] 处具有斜率 c[1]，在 b[end] 处具有斜率 c[end] 的唯一三次样条。因此，spline(b,) 是 S 的基映射。

而且，spline(b,c,xi) 提供了插值样条在 xi 处的值。最后，spline 支持向量值数据。因此，接下来的一行代码就构造了具有断点序列 b 的三次样条对数据 (x,y) 的最小二乘逼近：

spline(b,vec(y)'/spline(b,eye(length(b)),vec(x)'))