數值分析系列（五）：高維插值

封存文章

本文由舊站數值分析系列重寫。原稿把若干直觀敘述誤寫成定理，亦列出缺乏一般性的距離與條件數公式；本版本改以清楚的函數類別與假設說明結果。

高維插值出現在電子結構計算、多資產衍生工具、參數化偏微分方程、不確定性量化與 surrogate modeling。若在每個坐標方向放 $m$ 個節點，完整 tensor grid 便有

N=m^d

個節點。當維度 $d$ 增加，這個指數增長很快超出儲存與計算能力。但「高維必然不可解」也過於悲觀：許多實際函數只依賴少數變數、低階交互作用、低秩結構或較平滑的混合導數。現代方法的關鍵不是消除維度，而是找出並利用這些結構。

一、維度災難的準確含義

1. Tensor-product 網格

假設一維插值要以網格寬度 $h$ 達到某個誤差，而每軸約需 $h^{-1}$ 個節點，則 $d$ 維完整網格需要

N\asymp h^{-d}.

若某個具有固定 Lipschitz 常數的函數類別要求 $h=O(\varepsilon)$ 才令 $L^\infty$ 誤差不超過 $\varepsilon$ ，便得到

N=O(\varepsilon^{-d}).

這個結論需要函數類別與 regularity 假設。對「所有連續函數」，不能單憑 $\varepsilon$ 指定統一節點數，因為不同函數有不同 modulus of continuity。

2. 高維幾何

$d$ 維單位球的體積為

V_d=\frac{\pi^{d/2}}{\Gamma(d/2+1)},

在 $d$ 很大時趨近零。更直接的 shell concentration 是

\frac{\operatorname{Vol}(B_d(1-\delta))} {\operatorname{Vol}(B_d(1))} =(1-\delta)^d\longrightarrow0

對任何固定 $0<\delta<1$ 成立；因此均勻分布於球內的點，大部分落在靠近邊界的薄殼。

原文曾寫

\frac{\max_{x\in B_d}\|x\|} {\min_{x\in B_d}\|x\|}\to1,

但因為球包含原點，分母是零，該式沒有意義。常見的「距離集中」結果是針對指定概率分布下隨機點之間的距離，且需要說明正規化與概率收斂模式；它不是所有高維資料的普遍定律。

3. 資料稀疏

當樣本數固定而維度上升，fill distance

h_{X,\Omega} =\sup_{x\in\Omega}\min_{x_i\in X}\|x-x_i\|

通常很難保持小。局部插值器找不到真正鄰近的點，kernel length scale、噪聲和邊界效應都變得更敏感。這是幾何覆蓋問題，而不只是「電腦不夠快」。

二、Sparse grid

Smolyak sparse grid 不保留所有 tensor-product level 組合，而只保留總層級較低的組合。對具有足夠 bounded mixed derivatives 的函數，它可把節點數由完整網格的指數量級降低為近似

N_L=O(2^L L^{d-1}),

並取得帶 logarithmic factor 的誤差界。一個常見形式是

\|f-A_Lf\| \le C\,N^{-k}(\log N)^{(d-1)(k+1)},

但範數、smoothness space、邊界處理與 $k$ 的定義會隨定理版本而異。常數亦可能隨 $d$ 顯著增加，所以 sparse grid 是減輕、不是神奇消除維度災難。

它最適合：

各坐標方向有可分層的一維 rule；
mixed regularity 足夠；
維度中等而函數評估昂貴；
可用 hierarchical surplus 做 adaptivity。

若函數在斜向低維 manifold 上變化，固定坐標 sparse grid 未必有效，此時需 anisotropic weights、旋轉坐標或 dimension reduction。

三、Kernel 與 radial basis interpolation

給定節點 $x_1,\ldots,x_N$ ，kernel interpolant 可寫為

s(x)=\sum_{i=1}^{N}\alpha_iK(x,x_i),

係數由

K_X\alpha=y,\qquad (K_X)_{ij}=K(x_i,x_j)

決定。Gaussian、Matérn 和 compactly supported RBF 都可處理 scattered data，而且公式本身不需要 tensor grid。

但計算成本仍與樣本數有關：

dense factorization 通常需 $O(N^3)$ ；
儲存 kernel matrix 需 $O(N^2)$ ；
length scale 過大或節點過近會令矩陣病態；
在高維中，若資料沒有低維結構，取得足夠覆蓋仍需大量樣本。

因此「公式不顯式含 $d$ 」不代表 sample complexity 與維度無關。大型問題會使用 local kernels、partition of unity、迭代線性求解、 Nyström／random-feature 近似或 inducing points。

若加入正則化，kernel ridge regression 解

\min_{f\in\mathcal H_K} \sum_{i=1}^{N}(f(x_i)-y_i)^2 +\lambda\|f\|_{\mathcal H_K}^2,

它與精確插值不同： $\lambda>0$ 允許對噪聲資料平滑，亦改善數值條件。

四、有效維度

1. Functional ANOVA

若輸入分布是乘積測度，平方可積函數可分解為

f(x) =f_\varnothing +\sum_i f_i(x_i) +\sum_{i<j}f_{ij}(x_i,x_j) +\cdots+f_{1\cdots d}(x).

在標準 side conditions 下，各分量互相正交並且分解唯一。若大部分變異只來自主效應與少數二階交互作用，就可建立低階 surrogate，而不用解析所有 $2^d$ 個分量。

這個想法亦連接 Sobol sensitivity indices：先辨識重要變數與交互作用，再把採樣資源集中到它們。

2. Active subspace

若 $f$ 可微，定義 gradient second-moment matrix

C =\mathbb E_\rho[ \nabla f(X)\nabla f(X)^\mathsf T].

特徵分解

C=W\Lambda W^\mathsf T

若在第 $r$ 個特徵值後出現明顯 gap，便可嘗試以

f(x)\approx g(W_r^\mathsf Tx)

取代原 $d$ 維函數。這只捕捉重要的線性組合；對彎曲 manifold、不連續函數或 noisy gradient，效果未必理想。

3. 低秩與可分結構

另一條路是尋找

f(x_1,\ldots,x_d) \approx \sum_{\ell=1}^{R} \prod_{j=1}^{d}g_{\ell j}(x_j),

或使用 tensor train、proper orthogonal decomposition 等低秩表示。若所需 rank $R$ 隨精度溫和增長，高維問題便可能可解；若 rank 本身指數增長，低秩方法同樣失效。

五、自適應採樣

固定網格在平滑區浪費節點，在尖峰或邊界層又可能不足。Adaptive interpolation 會以 error indicator 或 posterior uncertainty 選擇下一個點，例如：

由少量 space-filling design 開始；
擬合 surrogate；
用 cross-validation、hierarchical surplus 或 acquisition function 找出不確定區域；
在該處評估昂貴函數；
重複至誤差或預算達標。

可靠的 stopping rule 應用獨立驗證點或 problem-specific quantity of interest，而不是只看 training residual。

六、與機器學習的關係

RBF network 形式上是

\hat f(x) =\sum_{i=1}^{M}w_i\, \phi\!\left( \frac{\|x-c_i\|}{\ell_i} \right),

中心 $c_i$ 、尺度 $\ell_i$ 與權重可由資料學習。它與 RBF interpolation 的差別在於 $M$ 不必等於樣本數，而且通常以經驗風險加正則化訓練。

Gaussian process、kernel regression、neural operator 與深度網絡都可以作高維 surrogate，但能否 generalize 取決於資料分布、目標函數的結構與模型 bias。「深度模型處理高維」不是對維度災難的普遍數學豁免。

七、簡單實驗設計

比較方法時，不應只在一個隨機函數上報 RMSE。可以建立三組測試：

可分函數： $f(x)=\sum_j a_jg_j(x_j)$ ；
低階交互作用： $f(x)=\sum_{j<k}a_{jk}g_{jk}(x_j,x_k)$ ；
真正高階耦合： 例如窄 ridge、旋轉後的尖峰或高 rank tensor。

對每一組，改變維度與樣本預算，記錄：

held-out $L^2$ 與 $L^\infty$ 誤差；
function evaluation 次數；
fitting time 與 peak memory；
matrix condition estimate；
對 noise 和 hyperparameter 的敏感度。

這樣才可看到一種方法究竟利用了甚麼結構，而不是把個別成功例子外推成普遍性能。

八、應用

不確定性量化： 用 sparse polynomial／grid surrogate 傳播多個輸入參數；
參數化 PDE： 近似解場或 quantity of interest 對材料、邊界與幾何參數的映射；
金融計算： 對多資產、隨機波動率或風險因子建立近似，但需保持 no-arbitrage 與 tail accuracy；
量子化學： 高維 potential-energy surface 常依靠 symmetry、locality、 permutational invariance 與低階 many-body expansion，而不是一般黑箱插值；
設計最佳化： 以 surrogate 代替昂貴 simulator，再用 active learning 精化 Pareto front 附近。

結論

高維插值沒有單一萬用方法。完整 tensor grid 的成本是 $m^d$ ，但 sparse grid、kernel、ANOVA、active subspace、低秩表示與 adaptive sampling 都可在特定結構下減少有效複雜度。選方法前應先問：

函數具有哪一類 smoothness？
重要性集中於少數變數、低階交互作用，還是低維 ridge？
樣本是否有噪聲？函數評估成本是多少？
所需誤差是全域、局部，還是某個 quantity of interest？
驗證集是否代表實際輸入分布？

維度災難不是一句「點數太多」便結束的警告，而是一個迫使我們明確描述函數類別、資料幾何與可利用結構的建模問題。

參考資料

Bellman, R. (1957). Dynamic Programming. Princeton University Press.
Bungartz, H.-J., & Griebel, M. (2004). Sparse grids. Acta Numerica, 13, 147–269.
Ledoux, M. (2001). The Concentration of Measure Phenomenon. AMS.
Buhmann, M. D. (2003). Radial Basis Functions: Theory and Implementations. Cambridge University Press.
Fasshauer, G. E. (2007). Meshfree Approximation Methods with MATLAB. World Scientific.
Schaback, R. (1995). Error estimates and condition numbers for radial basis function interpolation. Advances in Computational Mathematics, 3, 251–264.
Fornberg, B., & Flyer, N. (2015). A Primer on Radial Basis Functions with Applications to the Geosciences. SIAM.
Constantine, P. G. (2015). Active Subspaces: Emerging Ideas for Dimension Reduction in Parameter Studies. SIAM.