數值分析系列四：超越多項式的插值

這是數值分析系列第四篇，主題是超越全局多項式的插值，重點放在徑向基底函數（Radial Basis Function，RBF）。

想像氣象學家由散佈歐洲的測站重建溫度場。資料不在規則網格，而且有兩個以上空間維度；一維高次多項式的工具不能直接搬過來。RBF 的吸引力，在於它只需要點與中心之間的距離，能自然處理 scattered data。

傳統插值作為起點

線性插值

對 $(x_0,y_0)$ 與 $(x_1,y_1)$ ，

p_1(x) =y_0+\frac{y_1-y_0}{x_1-x_0}(x-x_0).

若 $f\in C^2$ ，區間長度 $h=|x_1-x_0|$ ，則

\|f-p_1\|_\infty \leq\frac{h^2}{8}\|f''\|_\infty.

多維 multilinear interpolation 通常需要結構網格；單一 cell 有 $2^d$ 個角點，但這不代表所有線性或 simplex 方法都固定需要 $2^d$ 點。

多項式插值

Lagrange 形式為

p_n(x)=\sum_{i=0}^{n}y_i \prod_{j\ne i}\frac{x-x_j}{x_i-x_j},

Newton 形式則以 divided differences 建立：

p_n(x)=a_0+a_1(x-x_0) +a_2(x-x_0)(x-x_1)+\cdots.

高次等距節點可能出現 Runge 現象，monomial Vandermonde 基底亦可高度病態。不過誤差並非對所有函數都「隨 $n$ 指數增長」，Chebyshev 節點亦不是唯一穩定方法；節點、基底、函數正則性與算法共同決定表現。

RBF 插值

給定中心 $X=\{x_1,\ldots,x_n\}\subset\mathbb R^d$ ，RBF 插值寫成

s(x)=\sum_{j=1}^{n}c_j\phi(\|x-x_j\|),

並要求 $s(x_i)=y_i$ 。常見核包括：

Gaussian： $\phi(r)=e^{-(\varepsilon r)^2}$ ；
inverse multiquadric： $\phi(r)=(1+(\varepsilon r)^2)^{-1/2}$ ；
Matérn kernels；
polyharmonic splines，如 $r^k$ 或 $r^k\log r$ ，視維度與階數而定。

它不需要 mesh，適合天氣站、污染感測器、地形、點雲及不規則量度。

RBF 並沒有消除維度詛咒

RBF 免除了 tensor grid，但高維空間要保持相同 fill distance，樣本數仍會快速增加。建立 dense kernel matrix 需 $O(n^2)$ 記憶體，直接求解通常需 $O(n^3)$ 時間，而 $m$ 個新點的 dense 評估需 $O(mn)$ 。距離計算本身亦隨維度增加。

因此，RBF 改變了資料幾何與離散方式，不代表複雜度「不受維度影響」。

由 variational problem 到 native space

在所有穿過資料的函數中，可選擇 native-space norm 最小者：

\min_f\|f\|_{\mathcal N_\phi} \quad\text{subject to}\quad f(x_i)=y_i.

對適當正定平移不變核，其 native norm 可透過 Fourier transform 表示：

\|f\|_{\mathcal N_\phi}^2 =\int_{\mathbb R^d} \frac{|\widehat f(\omega)|^2} {\widehat\phi(\omega)}\,d\omega,

忽略依 Fourier convention 而異的常數。Representer property 令最小 norm 解落在核截面的 span：

s(x)=\sum_{j=1}^nc_j\phi(\|x-x_j\|).

插值矩陣為

A_{ij}=\phi(\|x_i-x_j\|), \qquad Ac=y.

若 $\phi$ 對互異點 strictly positive definite， $A$ 正定，解唯一。某些 conditionally positive definite RBF 則需要額外多項式項與 moment constraints，不能直接套用同一結論。

「 $H^1$ norm 產生 splines、 $L^2$ norm 產生 polynomials」是過度簡化；具體 minimizer 取決於導數 semi-norm、邊界條件及 null space。

Power function 與誤差

令

k_X(x)= \begin{bmatrix} \phi(\|x-x_1\|)&\cdots& \phi(\|x-x_n\|) \end{bmatrix}^{T}.

Power function 定義為

P_X(x)^2 =\phi(0)-k_X(x)^TA^{-1}k_X(x).

對 $f\in\mathcal N_\phi$ ，

|f(x)-s_f(x)| \leq P_X(x)\|f\|_{\mathcal N_\phi}.

$P_X$ 由點集與核決定，不需要知道資料值，適合協助 adaptive sampling；但它是 native-space worst-case bound，不等於觀測雜訊下的機率置信區間。

基本實作：

import numpy as np
from scipy.spatial.distance import cdist
from scipy.linalg import solve

def rbf_matrix(X, Y, phi):
    return phi(cdist(X, Y))

def rbf_interpolation(X, y, X_eval, phi):
    A = rbf_matrix(X, X, phi)
    coeffs = solve(
        A, y,
        assume_a='pos',
        check_finite=True
    )
    K = rbf_matrix(X_eval, X, phi)
    return K @ coeffs

有雜訊時，精確穿過所有點可能過度擬合，通常加入 regularization：

(A+\lambda I)c=y.

點集幾何與收斂

Fill distance 與 separation distance 為

h_X=\sup_{y\in\Omega}\min_{x\in X}\|x-y\|, \qquad q_X=\frac12\min_{i\ne j}\|x_i-x_j\|.

$h_X$ 衡量域內最差「空洞」， $q_X$ 衡量點是否過度聚集。Mesh ratio $h_X/q_X$ 反映點集均勻程度。

對有限平滑核，常見誤差形式是

\|f-s\|_{L^\infty(\Omega)} \leq C h_X^\beta\|f\|_{\mathcal N_\phi},

其中 $\beta$ 與核、維度、域及函數空間有關。對 Gaussian 及足夠規則的函數與點集，可取得非常快甚至指數型界線，例如示意性的

\|f-s\|_\infty\leq Ce^{-c/h_X},

但不是所有平滑函數、域、 $\varepsilon$ 與點集都共享同一公式。

存在、唯一與最優恢復

若 $\phi$ strictly positive definite，對互異點及任意資料值，插值矩陣正定：

v^TAv =\sum_{i,j}v_iv_j\phi(\|x_i-x_j\|)>0 \quad(v\ne0).

因此係數唯一。插值函數亦是所有符合資料約束的 native-space 函數中 norm 最小者：

\|s_f\|_{\mathcal N_\phi} =\min\{ \|g\|_{\mathcal N_\phi}:g(x_i)=f(x_i) \}.

對 Matérn kernel，其 native space 與某個 Sobolev／Bessel potential space 等價；常用參數化下大致為 $H^{\nu+d/2}(\mathbb R^d)$ 。原文把特定 $(1+\varepsilon r)e^{-\varepsilon r}$ 同時當成任意 $\nu$ 的公式並不準確；該閉式對應特定半整數 smoothness。

形狀參數與穩定—準確取捨

Gaussian 的 $\varepsilon$ 控制基底寬度。 $\varepsilon$ 很小時核很「平」，近似可能非常準確，標準基底的 $A$ 卻極度病態； $\varepsilon$ 太大時基底太窄，亦可能失去整體近似能力。

Condition number 的實際增長依點集與尺度，不應把

\kappa(A)=O( \varepsilon^{-2d}e^{C/\varepsilon^2} )

當作普遍精確定律。重要訊息是 flat limit 需要穩定基底及高品質線性代數。

普通 pivoted QR 只是在既有病態矩陣上作較穩定分解，不等於 RBF-QR。真正 RBF-QR／RBF-GA 方法透過解析基底轉換處理 flat limit；不能用以下三行普通 QR 取代：

Q, R, piv = scipy.linalg.qr(
    A, pivoting=True
)

大型問題可考慮 compactly supported RBF、partition of unity、domain decomposition、iterative solvers、fast summation 或 multilevel 方法。

PDE collocation

對

Lu=f\quad\text{in }\Omega, \qquad Bu=g\quad\text{on }\partial\Omega,

可令

u_h(x)=\sum_jc_j\phi(\|x-x_j\|)

並在內點施加 $L\phi_j$ 、邊界點施加 $B\phi_j$ 。矩陣每一列必須由算子作用於以同一組中心建立的基底函數，不能只呼叫 L(construct_rbf_matrix(X_interior))，因為算子通常需要對空間座標求導。

示意接口：

def rbf_pde_system(
    X_interior, X_boundary, centers,
    apply_L, apply_B, rhs, boundary_rhs
):
    A_i = apply_L(X_interior, centers)
    A_b = apply_B(X_boundary, centers)
    A = np.vstack([A_i, A_b])
    b = np.concatenate([
        rhs(X_interior),
        boundary_rhs(X_boundary)
    ])
    return A, b

PDE collocation 還需處理 boundary geometry、operator conditioning、stability、conservation 及 convergence；mesh-free 不等於 analysis-free。

與 machine learning 的連結

RBF 與 kernel ridge regression、Gaussian processes、SVM kernels 及 RBF neural networks 有密切關係，但目標與不確定性解讀不同。修正後的 PyTorch 層為：

class RBFLayer(nn.Module):
    def __init__(self, in_features, out_features):
        super().__init__()
        self.centers = nn.Parameter(
            torch.randn(out_features, in_features)
        )
        self.log_sigmas = nn.Parameter(
            torch.zeros(out_features)
        )

    def forward(self, x):
        diff = (
            x.unsqueeze(1)
            - self.centers.unsqueeze(0)
        )
        dist_sq = torch.sum(diff**2, dim=-1)
        sigma_sq = (
            torch.exp(self.log_sigmas)
            .unsqueeze(0)**2
        )
        return torch.exp(
            -dist_sq / (2 * sigma_sq)
        )

原文使用 init 及 super().init()，並不能正確初始化 PyTorch module。

Shape parameter 選擇

可按 leave-one-out cross-validation、generalized cross-validation、held-out prediction error 或物理尺度選擇 $\varepsilon$ 與 regularization：

def choose_shape_parameter(
    X, y, candidates, score
):
    results = [
        (score(X, y, epsilon), epsilon)
        for epsilon in candidates
    ]
    return min(results)[1]

只最小化 condition number 會偏好數值容易、但可能近似很差的核；準確率與穩定性應共同評估，而且資料標準化要在選參數前固定。

最佳實務

先縮放座標與輸出，說明距離 metric；
檢查重複或近重複點及 $q_X$ ；
在有雜訊資料使用 regularization，不盲目精確插值；
以 cross-validation 或獨立測點選核與形狀參數；
同時報告 residual、held-out error、condition number 與敏感度；
大型問題使用局部／稀疏方法；
不把 power function 當成資料噪音的可信區間；
對 PDE 驗證 boundary conditions、守恆與 mesh refinement。

結語

RBF 插值把 scattered data、kernel Hilbert space 與數值線性代數連在一起。它毋須規則 mesh，能配合複雜幾何，也可延伸至 surrogate modeling 與 PDE。

它並沒有免費消除維度詛咒或計算成本。核、形狀參數、點集幾何、regularization 與穩定基底共同決定結果。真正成熟的 RBF 工作，不只展示一張平滑曲面，還要說明表面在哪裏有資料支持、哪裏對建模選擇敏感。

參考文獻

Kincaid, D. & Cheney, W. (2009). Numerical Analysis: Mathematics of Scientific Computing, 3rd ed. AMS.
Buhmann, M. D. (2003). Radial Basis Functions: Theory and Implementations. Cambridge University Press.
De Marchi, S. & Perracchione, E. (2018). Lectures on Radial Basis Functions. University of Padua.
Duchon, J. (1977). Splines minimizing rotation-invariant semi-norms in Sobolev spaces.
Fasshauer, G. E. (2007). Meshfree Approximation Methods with MATLAB. World Scientific.
Liu, J., Wang, F. & Nadeem, S. (2023). A new type of radial basis functions for problems governed by partial differential equations. PLOS ONE, 18(11).
Barrodale, I. & Zala, C. (1999). Mapping scattered data in three dimensions using radial basis functions.
Meinguet, J. (1979). Multivariate interpolation at arbitrary points made simple.