Syslab使用手册
⋮
2026a

# tpaps

薄板平滑样条

函数库: TyCurveFitting

# 语法

st, P = tpaps(x,y)
st, P = tpaps(x,y,p)

# 说明

st, P = tpaps(x, y) 返回对给定数据位点 x[:,j] 以及给定数据值 y[:,j] 的 st 型的薄板平滑样条 f 以及最后的样条结果中使用的平滑参数的值，无论有没有指定 p。P 值在从 st,P = tpaps(x, y)开始得到对 p 的好的初始估计的实验中有用。x[:,j] 必须为平面上的不同点，其值可以为标量、向量、矩阵或多维数组，数据值的个数必须与数据位点的个数相同。

薄板平滑样条 f 为加权和

的唯一最小化子，其中 E(f) 为误差测度

R(f) 为粗糙测度

这里积分在整个 上取， 表示 z 的所有元素的平方和， 表示 f 对其第 i 个变元的偏导，因此该积分包含了 f 的二阶偏导。该函数选择平滑参数 p 使得 (1-p)/p 与矩阵 A 的对角线元素的均值相同，其中 A+(1-p)/p*eye(n) 为要构造的光滑样条的 n 个系数的线性系的系数方程，这保证了结果位于插值的两个极端之间（即 p 接近 1，系数矩阵基本上是 A)以及完全平滑（即 p 接近 0，稀疏矩阵基本上是单位阵的积）。这为 p 提供了好的初始假设。示例

st, P = tpaps(x, y, p) 也指定了平滑参数 p，该数介于 0 到 1 之间。当平滑参数从 0 变到 1 时，平滑样条就从对数据的最小二乘估计的线性多项式（p 为0）变到对数据的薄板样条插值（p 为1）。示例

# 示例

从噪声数据复原准确的光滑值

接下来的代码从一个光滑函数在31个任意选取的数值位点得到其取值，将这些值加上任意噪声，然后使用tpaps复原准确的光滑值。为了展现tpaps在这种情况下表现多好，代码除了平滑样条外还绘制了精确值（用黑球表示）以及从光滑值指向对应噪声值的箭头。

using TyPlot
using TyCurveFitting
using TyMath
rng = MersenneTwister(1234)
nxy = 31
xy = 2 .*(rand(rng,2,nxy).-.5)
vals = sum(xy.^2, dims=1)
noisyvals = vals + (rand(rng,size(vals,1), size(vals,2)).-.5)./5
st = tpaps(xy, noisyvals)[1]
fnplt(st)
hold("on")
avals = fnval(st, xy)
plot3(xy[1,:], xy[2,:], vals[:], "wo", markerfacecolor="k")
quiver3(xy[1,:], xy[2,:],avals[:],zeros(1,nxy)[:],zeros(1,nxy)[:],(noisyvals-avals)[:],color="r", length=5)
hold("off")
使用薄板样条插值构造地图

接下来的代码对向量值数据值使用薄板样条插值来构造了地图，从平面到平面，将单位球 近似地映到单位圆盘 。

using TyPlot
using TyCurveFitting
n = 64
t = LinRange(0, 2 * pi, n + 1)
t = t[1:end-1]'
values = [cos.(t); sin.(t)]
centers = values ./ repeat(maximum(abs.(values), dims=1), 2, 1)
st = tpaps(centers, values, 1)[1]
fnplt(st)
axis("equal")

注意此处选择平滑参数为 1 以得到插值。

# 输入参数

x — 数据位点
一维数组 | 元组

拟合数据值 y 的数据位点，单变元时指定为一维数组，多变元时指定为元组。样条 f 使用每个数据位点 x 作为节点，且使得对于所有的 j，满足 f(x[j])=y[:,j]。

在多变元情形应为网格数据，可以指定x为对每个变元维度指定数据位点的元组：f(x1[i],x2[j],...,xn[k])=y[:,i,j,...,k]。

数据类型： Int | Float

y - 拟合数据
向量 | 矩阵 | 多维数组

在样条构造时用于拟合的数据值，为向量、矩阵或数组。数据值 y[:,j]可以为标量、矩阵或 n 维数组。在同一数据位点 x 给定的数据值将会被平均。

数据类型： Int | Float

p - 平滑参数
在 [0,1] 中的标量 | 向量 | 数组元组 | 空数组

光滑参数，单变元情形指定为从 0 到 1 的标量值，多变元情形为值元组。可以通过将 p 以向量形式传入以提供粗糙测度权 的值。为了在多元情形提供粗糙测度权，使用向量元组。如果输入空数组，则函数将会给予数据位点 x 为 p 选择默认值，为粗糙测度权 选取默认值 1。

平滑参数决定了在使得 f 变得光滑以及使 f 与数据接近之间的权重。对于 p=0，f 为数据的最小二乘直线拟合。对于 p=1，f 为变化的，或者说自然的，三次样条插值。当 p 从 0 变到 1 时，平滑样条就从一个极端变到另一个。

p 的范围一般在 附近，其中 h 为数据位点的平均间隔。函数默认会在该范围中选择 p。对于等间距数据，指定 可以得到接近的拟合， 指定 可以得到让人满意的光滑拟合。可以输入 p>1，但这样的选择将会使得平滑样条甚至比变化的三次样条插值还要粗糙。

如果输入 p 为负或者为空，函数将会使用默认的 p。

可以通过以向量形式提供p来指定沿着平滑参数的粗糙测度权 。该向量必须与 x 有同样的大小，其第 i 个元素为 在区间 (x[i-1],x[i]), i=2:length(x) 上的取值。输入向量 p 的第一个元素为希望为平滑参数 p 指定的值。通过提供粗糙测度权，可以使光滑样条在区间的不同部分更加光滑（如果权值较大）或者更加接近数据（如果权值较小）。粗糙测度权必须是非负的。

如果对选择 p 存在困难但是对 y 中的噪声维度有感觉，可以使用 spaps(x, y, tol)。该函数选择 p 以使在误差测度不超过 tol 的基础上让噪声测度越小越好。这种情况下，误差测度常常与 tol 的指定值相同。

数据类型： Float

# 输出参数

st - 样条结构体
样条结构体

样条，以带有以下元素的结构体形式返回。

Form-样条的类型
"st-tp00" | "st-tp10" | "st-tp01" | "st-tp"

样条的类型，返回为"st-tp00"，"st-tp10"，"st-tp01"或"st-tp"。

数据类型： String

Centers-位点序列
矩阵 | 数组

位点序列，单变元时以矩阵，多变元时以数组形式返回。

数据类型： Int | Float

Coefs-多项式系数

每一分段多项式的系数，单变元时为矩阵，多变元时为数组。

数据类型： Int | Float

Ncenters - 中心的数目
标量 | 向量

位点序列的数目。

数据类型： Int

Number - 多项式分段的段数
标量 | 向量

描述样条的多项式分段的数目，单变元时以标量，多变元时以由各变元分段数目组成的向量形式返回。

数据类型： Int

Dim - 维数
标量

目标函数的维数，以标量形式返回。

数据类型： Int

Interv - 基本区间
数组元组

st型的基本区间，其中包含了所有给定的中心，以数组形式返回。

数据类型： Int | Float

P - 平滑参数
标量 | 数组元组

用来计算样条的平滑参数，单变元时以标量，多变元时以标量值的数组元组形式返回。P 在 0 到 1 之间。

数据类型： Float

# 限制

平滑样条的构造包含了求解有与数据点相同个数的位置量的线性系统。因为这个线性系统的矩阵不是稀疏的，这样的求解可能花很长时间，即使现在的算法在超过 728 个数据点时就会使用迭代策略。迭代的收敛速度与 p 的选择强相关，如果 p 越大，收敛越慢。因此，对于大规模问题，推荐仅仅在时间成本可负担时使用插值，即让 p 与 1 相同。

# 另请参阅

csaps | spaps