SciPy 平滑样条：GCV、自动节点与加权误差

作者：SciPy 文档贡献者；来源：Smoothing splines。本文按2026年10月5日读取时标示为 SciPy 1.18.0 的在线教程翻译整理，覆盖一维平滑、参数曲线、旧 FITPACK 接口与二维曲面。代码只是静态审查后的阅读示例，本次没有执行拟合。

插值曲线要求穿过每一个给定数据点。如果数据带有噪声，这样做会把随机波动也当成真实结构。平滑样条改为寻找一条接近观测值、又不过度弯曲的曲线。SciPy 提供的两条主要路线在惩罚项、节点选择和参数含义上不同，不能只看到“数值越大越平滑”就把参数互换。

经典平滑样条：用 λ 惩罚曲率

make_smoothing_spline 是对 H. Woltring 的经典 gcvspl 算法的净室重实现。对观测 (x_j, y_j)，它寻找一个自然三次样条 g，最小化：

Σ w_j |y_j − g(x_j)|² + λ ∫[x_1,x_n] (g″(u))² du。

第一项惩罚偏离数据，第二项惩罚二阶导数带来的弯曲。节点位于数据的横坐标上；给定非负的 λ 后，算法优化样条系数。λ=0 得到自然插值样条；当惩罚相对数据项足够大时，曲线趋近于直线。这个“大”取决于数据和坐标尺度，不能把某个固定数值当作普适阈值。

接口返回可调用的 BSpline 对象。lam=None 或省略 lam 时，使用广义交叉验证 GCV 选择惩罚参数。其思想与留一交叉验证有关，但采用更高效的计算方式；这不是对所有数据分布或区间外预测准确性的保证。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import make_smoothing_spline

x = np.arange(0, 2*np.pi + np.pi/4, 2*np.pi/16)
rng = np.random.default_rng(20261005)  # 编辑改动：固定随机种子
y = np.sin(x) + 0.4 * rng.standard_normal(size=len(x))
xnew = np.arange(0, 9/4, 1/50) * np.pi

for lam in [0, 0.02, 10, None]:
    spl = make_smoothing_spline(x, y, lam=lam)
    plt.plot(xnew, spl(xnew), "-.", label=f"lambda={lam}")
plt.plot(x, y, "o")
plt.legend()
plt.show()

原文展示了插值跟随噪声、较大惩罚压平曲线、GCV 接近其示例底层正弦曲线的图。这里不复述为本次实测结果：固定种子后样本已不同，曲线也不会与原图逐点一致。make_smoothing_spline 要求至少5个数据点；整理真实数据时还应核对横坐标顺序、有限值、重复点处理以及权重合法性。

该接口允许多维 y，用 axis 指定数据轴，要求 y.shape[axis] == len(x)，其他轴作为批次轴。这与下面自动选择节点的方法不同。

自动选择节点：用 s 限制拟合误差

make_splrep 为一维函数构造平滑样条，make_splprep 则构造多维空间中的参数曲线。这一系列重实现 P. Dierckx 的 FITPACK 方法。以三次样条为例，在内部节点 t_j 处计算三阶导数的跳跃：

D_j = g‴(t_j + 0) − g‴(t_j − 0)。

算法最小化 Σ |D_j|²，同时满足 Σ [w_j (g(x_j) − y_j)]² ≤ s。对 k 次样条，使用的是第 k 阶导数的跳跃。如果所有跳跃为零，整段曲线就成为一个多项式。因而非常大的 s 允许结果退化为单个最佳拟合多项式；它不像经典曲率惩罚在大 λ 下趋于直线。

s=0 对应插值；增大 s 通常允许更平滑、节点更少的结果。默认边界条件为 not-a-knot，而经典接口是自然三次样条。λ 是目标函数中的惩罚系数，s 是残差预算；二者不是同一个量，不能照抄数值。

算法从没有内部节点的样条开始。若误差约束不满足，就在偏差最大的区间增加节点、重新拟合，再检查约束，直到满足要求或达到节点数上限。默认最大节点向量长度是 len(x)+k+1。s=0 时直接构造相应插值节点向量，而非执行完整迭代。下面保留原文的概念性伪代码；它省略了最终平滑优化的细节，不能替代完整库算法：

for t in generate_knots(x, y, s=s):
    g = make_lsq_spline(x, y, t=t)
    if ((y - g(x))**2).sum() < s:
        break

节点生成器可以单独查看：

from scipy.interpolate import generate_knots
xk = np.arange(7)
yk = xk**4
knots = list(generate_knots(xk, yk, s=1))
interpolation_knots = list(generate_knots(xk, yk, s=0))
quadratic_knots = list(generate_knots(xk, yk, s=1, k=2))

原文三次样条示例从端点重复的 [0,0,0,0,6,6,6,6] 开始，随后加入内部节点；s=0 直接给出包含 2,3,4 的插值节点。通常节点来自数据横坐标，但偶数次数、s=0 或达到插值节点上限时，节点可位于数据点之间，例如二次示例中的 1.5、2.5、3.5、4.5。浮点舍入可能造成细微差异，小幅调整 s 也可能令节点数明显变化。

延续前面的同一组带噪数据，可比较两种误差预算：

from scipy.interpolate import make_splrep

plt.plot(xnew, np.sin(xnew), "-.", label="sin(x)")
plt.plot(xnew, make_splrep(x, y, s=0)(xnew), label="s=0")
plt.plot(xnew, make_splrep(x, y, s=len(x))(xnew), label=f"s={len(x)}")
plt.plot(x, y, "o")
plt.legend()
plt.show()

原文图中 s=0 追随随机起伏，较大的 s 更接近底层正弦，但区间外结果随 s 剧烈变化。因此，曲线“更顺”并不意味着可以可靠外推。选择 s 通常需要试验：当权重为观测标准差的倒数时，可从 m−√(2m) 到 m+√(2m) 附近寻找；当权重全为1时，s≈mσ² 是一种起点。它们是带假设的经验范围，不是自动质量认证。

核验勘误：两类接口的权重定义相反

原教程的文字注释把“哪个接口对权重平方”写反了。本文依据同页数学公式以及 make_smoothing_spline API 和 make_splrep API 修正如下：

不能把同一个权重数组原样交给这两类接口后声称保持了同一统计意义。加权残差定义、s 的尺度与数据单位必须一起核对。

二维轨迹不是二维曲面

对于一串平面上的点 p_j=(x_j,y_j)，目标可以是参数曲线 g(u)=(g_x(u),g_y(u))。这里自变量是沿曲线走过的参数 u，而不是必须把 y 写成 x 的单值函数。方法可推广到任意允许维度；惩罚与误差约束都再对坐标分量求和。

make_splprep 接收形状为 (维度数, 点数) 的数据，通常写成多个一维坐标数组的列表，返回 (spl, u)。默认按相邻点之间距离的累计值归一化，得到弦长参数化；也可以显式传入 u。库能提供默认值，却不能替使用者决定哪种参数化符合问题物理意义，不同参数化可能得到明显不同的曲线。

from scipy.interpolate import make_splprep

th = np.linspace(-0.2, np.pi/2 + 0.2, 21)
r = 3*np.sin(th)*np.cos(th) / (np.sin(th)**3 + np.cos(th)**3)
x_curve, y_curve = r*np.cos(th), r*np.sin(th)
rng_curve = np.random.default_rng(20261005)
xn = x_curve + 0.1*rng_curve.uniform(-1, 1, size=len(th))
yn = y_curve + 0.1*rng_curve.uniform(-1, 1, size=len(th))

spl, u = make_splprep([xn, yn], s=0)
smooth, u_smooth = make_splprep([xn, yn], s=0.1)
plt.plot(xn, yn, "o")
plt.plot(*spl(u), "--")
plt.plot(*smooth(u_smooth), "-")
plt.show()

这是原文笛卡尔叶形线例子的整理版：保留函数和噪声模型，新增固定种子、区分变量名。spl(u) 的两个坐标数组可用星号展开给 plot。评估更密的 u 可以显示曲线，但不能凭更密的绘图点制造更多观测信息。

批次处理：每条曲线可能需要不同节点

make_splrep 不像经典接口那样接受任意批次轴；参数曲线的二维坐标数组也不等于“一批独立曲线”。自动节点数取决于各组观测，批次中的节点向量可能长短不同。原文给出一个 BatchSpline 包装类，内部仍然逐条调用拟合。更直白的写法如下：

# 编辑改写：用显式循环表示独立批次，避免掩盖输出形状差异。
# y_batch.shape == (batch_size, len(x))
splines = [make_splrep(x, row, s=1.0) for row in y_batch]
values = np.stack([spline(xnew) for spline in splines], axis=0)

这里的 y_batch 是读者传入的业务数组，并非前一节已经生成的变量；实际数据需先准备。原包装类在调用时直接访问 x.shape，不能把普通 Python 标量或列表不经转换就当成已验证接口。本文改为显式循环，且没有把此处的任意 s=1.0 当成默认最佳参数。

旧 FITPACK 接口：维护旧代码时怎么读

教程仍介绍 splrep、splprep 和一元样条类，但将它们归为 legacy。官方当前没有计划移除这些接口，新的代码宜使用现代构造器。旧接口默认值不完全相同，尤其 s 会随是否传入权重以及具体函数变化；应明确写出 s。

splrep(x,y,s=...) 返回三元组 (t,c,k)：节点、系数和次数。默认三次；系数至少有 len(t)−k−1 个，一些旧接口补零至节点长度，多余系数不会参与求值。该表示可包装成 BSpline(*tck)。

from scipy.interpolate import splrep, BSpline
tck = splrep(x, y, s=0)
tck_s = splrep(x, y, s=len(x))
interpolated = BSpline(*tck)(xnew)
smoothed = BSpline(*tck_s)(xnew)

splev 可求值或用 der 求导，spalde 返回各阶导数，splint 计算定积分，sproot 可求满足条件的三次样条的根。下面将原文多段重复绘图代码合并，保留其运算；积分输出改成浮点数组，防止整数输入引起截断：

from scipy import interpolate

xs = np.arange(0, 2*np.pi + np.pi/4, 2*np.pi/8)
ys = np.sin(xs)
tck = interpolate.splrep(xs, ys, s=0)
xeval = np.arange(0, 2*np.pi, np.pi/50)

values = interpolate.splev(xeval, tck, der=0)
first_derivative = interpolate.splev(xeval, tck, der=1)
all_derivatives = interpolate.spalde(xeval, tck)
integrals = np.array(
    [interpolate.splint(0, float(v), tck) for v in xeval], dtype=float
) - 1.0
roots = interpolate.sproot(tck)

原文指出，sproot 可能漏掉近似区间边界上显而易见的零点，例如0。把拟合区间扩到 [−π/4, π+π/4] 后，原示例能找回0与π附近的根；这是其数值例子，不是本文的执行输出。应把边界情况纳入实际验证。

splprep([x,y],s=...) 返回 (tck,u)，其中每个坐标方向各有一组系数。转换成 BSpline 时要转置系数：

t_circle = np.arange(0, 1.1, 0.1)
xc = np.sin(2*np.pi*t_circle)
yc = np.cos(2*np.pi*t_circle)
tck_curve, u_curve = interpolate.splprep([xc, yc], s=0)
tt, cc, degree = tck_curve
curve = BSpline(tt, np.asarray(cc).T, degree)
xy = curve(u_curve)
x_back, y_back = xy.T

另一处核验修正：教程旧接口一节把 splprep 的默认 u 说成等距序列。实际 API 文档 给出的默认算法是累计点间距离后归一化，即弦长参数化。本文采用 API 定义；需要等距参数时应自己传入。

面向对象接口中，UnivariateSpline 可用非零 s 平滑；InterpolatedUnivariateSpline 总是插值；LSQUnivariateSpline 让调用方指定内部节点。对象本身可调用，也提供积分、导数和求根等方法。原文的非均匀节点示例如下：

exact = interpolate.InterpolatedUnivariateSpline(xs, ys)
internal_knots = [np.pi/2 - 0.1, np.pi/2 + 0.1,
                  3*np.pi/2 - 0.1, 3*np.pi/2 + 0.1]
custom = interpolate.LSQUnivariateSpline(xs, ys, internal_knots, k=2)
exact_values = exact(xeval)
custom_values = custom(xeval)

二维平滑曲面：区分散点与矩形网格

曲面是 z=g(x,y)，可看作一维样条的张量积。散点数据中，x[i]、y[i]、z[i] 成对表示一个观测；约束为 Σ [w_i(g(x_i,y_i)−z_i)]² ≤ s。矩形网格则由一维 x 和 y 的所有组合定义，二维 z[i,j] 对应 (x[i],y[j])，残差同时对两个网格索引求和。

FITPACK 自动在边界盒内布置节点，节点不一定落在原始数据点上。虽然形式上 s=0 对应插值，二维散点 FITPACK 过程并非通用插值器，可能产生意外结果；原文建议一般散点插值优先考虑 griddata，规则网格插值考虑 RegularGridInterpolator。坐标单位相差许多数量级时还应先考虑缩放，避免数值伪影。

散点接口 bisplrep 与 bisplev

bisplrep 接收一维 x,y,z，默认 kx=ky=3，返回 [tx,ty,c,kx,ky]。bisplev 在两个一维求值坐标的笛卡尔积上计算曲面，也支持指定偏导阶数。bisplrep 的默认 s=m−√(2m)；若权重为标准差倒数，可在 m±√(2m) 附近试选。

原文先在20×20网格取函数样本，再在70×70网格求值。以下保留步骤，并显式展平散点接口的输入：

x_edges, y_edges = np.mgrid[-1:1:21j, -1:1:21j]
xp = x_edges[:-1, :-1] + np.diff(x_edges[:2, 0])[0] / 2
yp = y_edges[:-1, :-1] + np.diff(y_edges[0, :2])[0] / 2
zp = (xp + yp) * np.exp(-6.0 * (xp*xp + yp*yp))

new_x_edges, new_y_edges = np.mgrid[-1:1:71j, -1:1:71j]
new_x = new_x_edges[:-1, :-1] + np.diff(new_x_edges[:2, 0])[0] / 2
new_y = new_y_edges[:-1, :-1] + np.diff(new_y_edges[0, :2])[0] / 2
surface_tck = interpolate.bisplrep(xp.ravel(), yp.ravel(), zp.ravel(), s=0)
new_z = interpolate.bisplev(new_x[:, 0], new_y[0, :], surface_tck)

plt.pcolormesh(new_x_edges, new_y_edges, new_z, shading="flat",
               cmap="RdBu_r", vmin=-0.25, vmax=0.25)
plt.colorbar()
plt.show()

这仍是对原文“强制插值”的教学演示，应读取拟合警告，不应据此推定所有散点都可顺利插值。此类二维方法也不适合直接用于大量图像像素的样条表示；图像任务应考虑 scipy.signal 或 scipy.ndimage 的相应算法。

SmoothBivariateSpline：不要隐藏警告

SmoothBivariateSpline 是散点曲面的对象接口，功能是过程式接口的一个子集，默认单位权重、s=m 和 kx=ky=3。原文对已知函数采样，比较默认平滑和 s=0。在原例中默认值过度平滑，s=0 更接近采样函数，却出现 FITPACK 内部困难的警告；原文通过全局忽略警告使绘图继续。本文明确删除了这一做法：

from scipy.interpolate import SmoothBivariateSpline

train_x, train_y = np.meshgrid(np.arange(-5, 5, 0.5),
                              np.arange(-5, 5, 0.5))
train_x, train_y = train_x.ravel(), train_y.ravel()

def z_func(x, y):
    return np.cos(x) + np.sin(y)**2 + 0.05*x + 0.1*y

train_z = z_func(train_x, train_y)
# 编辑改动：不调用 warnings.simplefilter("ignore")，保留警告。
interp_func = SmoothBivariateSpline(train_x, train_y, train_z, s=0.0)
smth_func = SmoothBivariateSpline(train_x, train_y, train_z, s=len(train_z))
test_x = np.arange(-9, 9, 0.01)
test_y = np.arange(-9, 9, 0.01)
interp_result = interp_func(test_x, test_y).T
smth_result = smth_func(test_x, test_y).T

如果自己的数据触发类似警告，应检查数据尺度、拟合条件与输出状态；原文建议考虑 bisplrep 并调整 nxest、nyest，而不是静默接受结果。示例求值网格大于训练区域，其边界外行为也不能当作学到了真实趋势。原图在该路径表现为边界值常量延拓，这不表示所有样条接口都使用相同外推策略。上例1800×1800求值网格会占用明显内存，可按业务所需分辨率缩小。

规则网格与球面坐标

RectBivariateSpline 针对规则矩形网格：一维 x、y 定义网格，z.shape 必须是 (len(x),len(y))。默认仍为双三次，但默认 s=0。原文建议始终显式指定 s：

from scipy.interpolate import RectBivariateSpline

x_grid = np.arange(-5.01, 5.01, 0.25)
y_grid = np.arange(-5.01, 7.51, 0.25)
xx, yy = np.meshgrid(x_grid, y_grid, indexing="ij")
z_grid = np.sin(xx**2 + 2.0*yy**2)
func = RectBivariateSpline(x_grid, y_grid, z_grid, s=0)
x_dense = np.arange(-5.01, 5.01, 1e-2)
y_dense = np.arange(-5.01, 7.51, 1e-2)
z_dense = func(x_dense, y_dense)
plt.imshow(z_dense)
plt.colorbar()
plt.show()

若数据使用球面坐标 r=r(θ,φ)，对应的 SmoothSphereBivariateSpline 与 RectSphereBivariateSpline 提供球面约束及极点连续性选项。应注意方位角的周期和极点条件，不能把 θ∈[0,π] 与 φ∈[0,2π] 都误解为同样的周期边界；具体约束应查所用类的 API。

实践中的核验顺序

先确认问题是一维函数、参数曲线还是二维曲面；再确认权重的统计意义和数据轴；最后选择边界条件、节点策略和误差预算。观察残差、节点数以及边界附近行为，并把区间外结果单独对待。GCV 是参数选择方法，低残差是拟合性质，都不是系统正确性的完整证据。

本文所列计算没有网络访问、外部命令、文件删除或硬编码秘密；实际风险主要来自错误权重、参数误用、隐藏警告、外推误判以及大网格的资源消耗。代码未执行，原文图和结果没有被冒称为本次实测，静态审查也不构成无漏洞或数值精度保证。

来源与归属：SciPy 文档贡献者；经典方法参考 H. Woltring，FITPACK 方法参考 P. Dierckx。在线文档版权标示 ©2008 The SciPy community。SciPy 项目采用 BSD 3-Clause 许可证；完整许可证声明可下载：SciPy BSD 3-Clause 完整文本。本文译文与配图依据另行授权制作。修改包括修正权重和旧接口参数化说明、固定随机种子、合并重复绘图、显式展平输入、保留警告、使用浮点积分输出；未改变算法为另一套未经说明的方法。