B样条插值奇异矩阵问题：为何需调整Cox-de-Boor条件？

问题背景 #

根据Wolfram Mathworld等B样条相关文献，Cox-de-Boor递归函数的0阶条件定义为：当阶数d_=0时，若knots_[k_] ≤ t_ < knots_[k_+1]则返回1.0，否则返回0.0。但在构建B样条插值线性系统时，该实现会产生奇异矩阵（例如4个点的案例中，矩阵右下角元素为0而非预期的1.0），触发LinAlgError。将条件改为knots_[k_] ≤ t_ ≤ knots_[k_+1]后，矩阵恢复正常，样条插值结果正确。需要解释为何要调整这个条件才能正确遍历到最后一个元素并得到正确结果。

完整示例代码 #

import numpy as np
import math
from geomdl import knotvector

def cox_de_boor( d_, t_, k_, knots_):
    if (d_ == 0):
        if ( knots_[k_] <= t_ <= knots_[k_+1]):
            return 1.0
        return 0.0
    
    denom_l = (knots_[k_+d_] - knots_[k_])
    left = 0.0    
    if (denom_l != 0.0):
        left = ((t_ - knots_[k_]) / denom_l) * cox_de_boor(d_-1, t_, k_, knots_)
    
    denom_r = (knots_[k_+d_+1] - knots_[k_+1])
    right = 0.0
    if (denom_r != 0.0):
        right = ((knots_[k_+d_+1] - t_) / denom_r) * cox_de_boor(d_-1, t_, k_+1, knots_)

    return left + right

def interpolate( d_, P_, n_, ts_, knots_ ):
    A = np.zeros((n_, n_))
    
    for i in range(n_):
        for j in range(n_):
            A[i, j] = cox_de_boor(d_, ts_[i], j, knots_)
    
    control_points = np.linalg.solve(A, P_)
    return control_points

def create_B_spline( d_, P_, t_, knots_):
    sum = MVector()
    for i in range( len(P_) ):
        sum += P_[i] * cox_de_boor(d_, t_, i, knots_)
    return sum
    
def B_spline( points_ ):
    d     = 3
    P     = np.array( points_ )
    n     = len( P )
    ts    = np.linspace( 0.0, 1.0, n )
    knots = knotvector.generate( d, n ) # len = n + d + 1
    
    control_points = interpolate( d, P, n, ts, knots)
    
    crv_pnts = []
    for i in range(10):
        t = float(i) / 9
        crv_pnts.append( create_B_spline(d, control_points, t, knots) )
    return crv_pnts


control_points = [ [float(i), math.sin(i), 0.0] for i in range(8) ]
cps = B_spline( control_points )

原因解析 #

1. 边界节点的特殊性：代码中使用的是夹紧型节点向量（Clamped Knot Vector）（由knotvector.generate(d, n)生成），这类向量的首尾节点会重复d+1次（比如n=4、d=3时，节点向量为[0,0,0,0,1,1,1,1]）。这种重复设计是为了让B样条曲线在端点处与第一个/最后一个控制点重合，满足插值的端点约束。
2. 原定义的边界遗漏：数学定义中的左闭右开区间knots_[k_] ≤ t_ < knots_[k_+1]是针对内部非重复节点的通用情况，目的是避免相邻0阶基函数在节点处同时取1。但在重复的边界节点处，当t_等于最后一个节点值（比如1.0）时，对于最后一个基函数j = n-1，knots_[k_+1]就是最后一个重复节点，此时t_ < knots_[k_+1]不成立，导致0阶基函数返回0，递归计算出的高阶基函数值也为0。
3. 矩阵奇异的直接原因：插值系统矩阵A的元素A[i,j]是第j个基函数在插值参数ts_[i]处的值。当ts_包含最后一个节点值时，原条件会让A的最后一行最后一列元素为0，破坏了矩阵的满秩性，导致矩阵奇异无法求解。
4. 修正条件的合理性：将条件改为包含等号后，当t_等于重复的边界节点时，最后一个基函数的0阶条件成立（knots_[k_] ≤ t_ ≤ knots_[k_+1]），递归计算出的高阶基函数在端点处会正确取1，确保矩阵A的对角元素非零，矩阵满秩可解。这种调整并非偏离数学定义，而是针对夹紧型节点向量的边界情况做的必要适配——数学定义的左闭右开是通用规则，而重复节点的边界需要特殊处理才能满足插值的端点需求。

内容的提问来源于stack exchange，提问作者Constantinos Glynos