考虑一个近似问题：假设有一个函数f(x)，x ∈ [-1, 1]^p，它的计算代价很高；因此，你希望构建一个新的函数f’，这个函数计算代价更低，但能够在[-1, 1]^p上准确地近似f，并且尽可能少地计算f。
这样的近似在统计学和其他领域有很多应用。例如，f可能表示一个概率密度函数，对其计算需要使用复杂的积分，而有效的近似可以为你提供计算可信区间的机制。
下面，我将介绍两种构建函数近似的强大技术：切比雪夫插值和自适应稀疏网格。
切比雪夫插值是近似的基本组成部分。如果f有一点平滑性（例如，利普希茨连续），那么切比雪夫插值将会收敛；如果f是高阶可微的或解析的，切比雪夫插值将会极快地收敛（比其他技术如三次样条要好得多）。
利普希茨连续（Lipschitz continuity）是描述函数光滑程度的一种数学概念。具体来说，一个函数如果在整个定义域内满足某个固定的“最大变化速率”，那么这个函数就是利普希茨连续的。
将切比雪夫插值与自适应稀疏网格结合起来，可以提供一种高效的机制来近似高阶函数。
让我们先从单维情况开始。
单变量函数的近似
有一种误解常常被提起，那就是高阶多项式插值存在问题。实际上，高阶多项式插值往往是完成任务的最佳工具，我们只需正确选择插值点。
虽然在等间距点上进行插值可能会带来灾难性的结果，但在切比雪夫节点上进行插值却具有优良的收敛特性。
- 在区间[-1, 1]上的切比雪夫节点示例。这些节点对应于半圆上的等间距点，因此在区间的端点附近密度更大。
示例1：龙格函数（Runge function）
考虑龙格函数：
等间距插值对近似这个函数的效果非常差（被称为龙格现象）。但如我们所见，切比雪夫插值能很好地处理这个函数。
龙格现象是指在用高次多项式进行插值时，特别是在等间距节点上进行插值时，插值多项式在区间边界处产生较大误差和震荡的现象。这个现象得名于德国数学家卡尔·龙格（Carl Runge），他在研究正弦函数的插值时首次发现了这一问题。
- 龙格函数的切比雪夫插值，使用9、17和33个点。
切比雪夫插值非常稳健。如果f在[-1, 1]上是利普希茨连续的，其切比雪夫插值将会收敛。
对于平滑函数，切比雪夫插值收敛得非常快。如果f有v个导数，并且第v个导数具有有界变差V，那么
其中 p_n 表示f的n个点的切比雪夫插值，并且
当f是解析函数时，收敛性更好
对于某个 ρ>1。
切比雪夫插值的收敛速度远远优于三次样条等方法，三次样条的收敛速度通常是O(n^-4)。正如Lloyd Trefethen所写的那样：
引入样条函数是一个误导：样条函数的真正优势不是它们在多项式失败时仍能收敛，而是它们更容易适应不规则点分布并且更具局部性。
多变量函数的处理
切比雪夫插值有一个自然的扩展，可以应用于多维情况。设X^n表示[-1, 1]上的n个切比雪夫节点的集合。那么它们的笛卡尔积
确定了一个唯一的多元切比雪夫插值p_{n^p}。与单维情况一样，当n足够大时，多元切比雪夫插值将很好地近似f；然而，就所需的计算次数而言，这将变得非常“昂贵”。将近似扩展到多维的简单方法并不是很高效。幸运的是，我们可以做得更好。
设
那么X^i定义了称为切比雪夫-高斯-洛巴托节点的嵌套点序列，并作为更高效插值的基本构建块。
- 从X^1到X^5的切比雪夫-高斯-洛巴托节点。
设ψ_j^i表示唯一的(m_i-1)阶多项式，其中
设V^i表示由基函数{ψ_j^i}（j=1,…,m_i）生成的向量空间。
- 基函数示例。
设ΔX?={}，ΔX^i=X^i \ X^{i-1}，并且
我们将使用向量空间的元素来近似函数
其中I = {i : i_1 + ··· +i_p ≤r}，r是控制近似精度水平的整数。在点上的f值
将唯一确定W^I中的插值。
- 从r=3到r=8的稀疏网格。
稀疏网格最初由俄罗斯数学家Sergey Smolyak发现，并用于开发更高效的数值积分规则。假设维度为p的函数f可微v次。使用密集正交规则，积分误差可以表示为
其中N_l是评估点的数量。当使用稀疏网格时，误差改进为
类似的性能改进也适用于内插。
自适应近似
虽然稀疏网格比密集网格更高效地近似多变量函数，但它们仍然远未达到最优。在大多数实际问题中，函数在某些区域需要更多的评估点以达到目标精度。对Smolyak稀疏网格进行一些修改，使其更具自适应性，可以使算法更加高效。
首先，考虑指标集I。我们不再要求I = {i : i_1 + ··· +i_p ≤r}，而是可以施加较弱的要求，即指标集只需要是可接受的：如果i∈I且i_k > 0，则i-e_k∈I，其中e_k表示单位向量，其(e_k)_k = 1。放宽到可接受性的要求允许网格的某些维度比其他维度具有更大的细化。
接下来，通过跟踪连续细化之间的差异，我们可以近似估计近似中的局部误差，并且可以调整稀疏网格算法，仅在尚未达到所需目标精度的区域添加细化点。
我们可以将这两个修改结合到一个贪心算法（greedy algorithm）中，不断细化在最不佳近似函数的区域，直到所有区域都达到所需的精度。
让我们看几个示例。我们将使用Python包bbai来拟合自适应稀疏网格插值。
示例2：局部自适应性
考虑这个函数：
该函数沿着半圆x2+y2=0.3绘制出一个锐利的脊线，最大值为10。让我们拟合一个稀疏网格来近似。
import numpy as np
def f(x, y):
t = np.abs(0.3 - x*x - y*y)
return 1/(t + 0.1)
# Fit a sparse grid to approximate f
from bbai.numeric import SparseGridInterpolator
ranges = [(0, 1), (0, 1)]
interp = SparseGridInterpolator(tolerance=1.0e-2, ranges=ranges)
interp.fit(f)
# Test the accuracy at some random points
N = 100
np.random.seed(0)
sample_pts = np.random.uniform(0, 1, size=(N, 2))
errs = [np.abs(interp.evaluate(x, y) - f(x, y)) for x, y in sample_pts]
print('mean_error: ', np.mean(errs))
# prints:
# mean_error: 0.0019457748570626515
print('max_error: ', np.max(errs))
# prints:
# max_error: 0.029750550242468213我们可以通过在函数的等高线图上叠加插值点来可视化构建的网格。正如我们所见，该算法自适应地在脊线周围的区域添加了更多的插值点。
示例3：按维度自适应
接下来让我们看一个函数，其中一个维度比另一个维度影响更大：
我们将再次拟合该函数，并在一些随机点上进行测试。
import numpy as np
def f(x, y):
dx = x - 0.5
dy = y - 0.5
mult = 0.01
return np.exp(-2 * dx*dx - mult * dy * dy)
# Fit a sparse grid to approximate f
from bbai.numeric import SparseGridInterpolator
ranges = [(0, 1), (0, 1)]
interp = SparseGridInterpolator(tolerance=1.0e-5, ranges=ranges)
interp.fit(f)
# Test the accuracy at some random points
N = 100
np.random.seed(0)
sample_pts = np.random.uniform(0, 1, size=(N, 2))
errs = [(np.abs(interp.evaluate(x, y) - f(x, y))) / f(x, y) for x, y in sample_pts]
print('mean_error: ', np.mean(errs))
# prints:
# mean_error: 3.265838801494769e-10
print('max_error: ', np.max(errs))
# prints:
# max_error: 1.7325188938498375e-09绘制网格点后，可以看到x维度的细化程度比y维度更高。
示例4：拟合10维高斯函数
让我们看看当我们尝试拟合一个高维函数时会发生什么。
定义一个10维高斯函数
# Define a Gaussian of 10 dimensions
import numpy as np
p = 10
mult = np.array([0.2699884 , 0.35183688, 0.296529 , 0.26805488, 0.20841667,
0.31774713, 0.21527071, 0.43870707, 0.47407319, 0.18863377])
offsets = np.array([0.79172504, 0.52889492, 0.56804456, 0.92559664, 0.07103606,
0.0871293 , 0.0202184 , 0.83261985, 0.77815675, 0.87001215])
def f(*args):
n = len(args[0])
t = np.zeros(n)
for j, arg in enumerate(args):
delta = arg - offsets[j]
t -= mult[j] * delta * delta
return np.exp(t)让我们测量拟合这个函数需要多长时间。
# Fit function
from bbai.numeric import SparseGridInterpolator
domain = [(0, 1)]*p
interp = SparseGridInterpolator(tolerance=1.0e-4, ranges=domain)
import time
t1 = time.time()
interp.fit(f)
t2 = time.time()
print('took {} seconds to fit sparse grid'.format(t2 - t1))
# Prints
# took 1.0177533626556396 seconds to fit sparse grid
print('num_points =', interp.points.shape[1])
# Prints
# num_points = 13256在我的电脑上拟合这个函数只需要一秒钟，只需要大约13k个点。让我们测试插值的准确性。
# Test accuracy at some random points
N = 100
X = np.random.uniform(size=(N, p))
args = [X[:, j] for j in range(p)]
f_true = f(*args)
f_approx = interp.evaluate(*args)
errs = np.abs((f_true - f_approx) / f_true)
print('mean_error:', np.mean(errs))
# Prints
# mean_error: 1.2278939934746143e-05
print('max_error:', np.max(errs))
# Prints
# max_error: 7.69422569392403e-05可以看到拟合非常准确。当然，并非每个高维函数都像高斯函数那样容易处理，但这表明自适应稀疏网格可以扩展到维数较高的情况，而密集网格在这种情况下会变得非常困难。
在统计学中的应用
考虑函数
其中X_1和X_2是均值为μ_ML，标准差为σ_ML的独立正态分布随机变量。
对于上下文，这个函数出现在计算默认的贝叶斯因子，用于对方差未知的正态分布数据的平均值进行假设检验时。
显然，通过采样，该函数易于近似；但通过切比雪夫插值，我们可以设计出更准确和高效的方法。
一些代数运算可以将f简化为
其中
设
那么
根据Cochran定理，\sqrt{2} \bar{z}服从标准正态分布，2 s_z^2服从自由度为1的卡方分布，并且这些分布是独立的。因此，
服从一个非中心t分布，其中心性参数为\sqrt{2} (μ_ML — t) / σ_ML。
设
其中Z表示具有中心性参数t的非中心t分布的随机可达性。给定G，计算f很简单。下面是近似G(因此是f)的步骤：
- 我们可以使用两个双指数exp-sinh正交规则[10]在t的左右两侧对任意给定点的G进行数值计算，以达到很高的精度。
- 使用切比雪夫插值来近似G。因为G是一个奇函数，我们只需要在[0，∞]上近似G。通过变换φ(s) = s/(1-s)和G(φ(1)) = π/2，我们可以很容易地将其转化为有限范围。
- 将把插值转换为等效的切比雪夫级数，然后截断小于最小阈值的系数。这给了我们一个更有效的近似值。
让我们尝试一下，用更简单的模拟方法测试它的准确性。
import numpy as np
# These are the coefficient of the approximation when represented as
# a truncated series of Chebyshev polynomials of the second kind
coefs = [
0.77143906574069909, 0.85314572098378538, 0.024348685879360281,
-0.080216391111436719, -0.016243633646524293, 0.014244249784927322,
0.0083546074842089004, -0.0013585115546325592, -0.0032124111873301194,
-0.00091825774110923682, 0.0007309343106888075, 0.00071403856022216007,
9.4913853419609061e-05, -0.00023489116699729724, -0.00017729416753392774,
-8.6319144348730995e-06, 7.1368665041116644e-05, 5.0436256633845485e-05,
1.8715564507905244e-06, -2.1699237167998914e-05, -1.6200449174481386e-05,
]
def g(s):
# Invert the mapping to go from the domain [0, ∞] -> [-1, 1]
t = s / (1 + s)
t = 2 * t - 1
# evaluate the approximation at t
return np.polynomial.chebyshev.chebval(t, coefs)
# Approximate f using Chebyshev polynomials
def f(mu, sigma, t):
t = np.sqrt(2) * (mu - t) / sigma
mult = 1.0
if t < 0:
mult = -1.0
t = -t
return 0.5 - mult * g(t) / np.pi
# Approximate f using a brute force simulation approach
N = 1000000
def f_sim(mu, sigma, t):
vals = np.random.normal(mu, sigma, size=(N, 2))
X1 = vals[:, 0]
X2 = vals[:, 1]
xbar = (X1 + X2) / 2.0
sx = np.sqrt(((X1 - xbar)**2 + (X2 - xbar)**2) / 2.0)
vals = 0.5 - np.arctan((xbar - t) / sx) / np.pi
return np.mean(vals)
# Compare the approximations against each other
mu_ml = 0.123
std_ml = 1.5
np.random.seed(0)
for t in [0.0, 0.01, 0.1, 0.2, 0.3, 5, 10, 100]:
ft = f(mu_ml, std_ml, t)
fpt = f_sim(mu_ml, std_ml, t)
rerr = np.abs((ft - fpt) / fpt)
print(t, f(mu_ml, std_ml, t), f_sim(mu_ml, std_ml,t), rerr)输出：
0.0 0.47059206926744995 0.4705617410005234 0.0006829181799177236
0.01 0.47297735696592846 0.47284440543778317 0.00044800674755964237
0.1 0.49449429024483676 0.4947273411923894 0.0001393914288991927
0.2 0.518426675093887 0.5184967680183266 2.2982920448911075e-06
0.3 0.5422530672885648 0.5421840161788148 0.0008338542554252333
5 0.943798414491211 0.9438468766277546 7.809732402546117e-05
10 0.9726189450498345 0.9726397857047364 4.339034701734454e-05
100 0.9973033311969053 0.9973052128447453 1.2957899605273032e-06可以看到两种近似方法的差异仅为百分之几的小数。
讨论
问题1：在切比雪夫节点上插值是近似光滑函数最准确的方法吗？
不是。但在实际应用中，它们可能是最佳解决方案，更准确的算法往往不会带来太大的改进。
其他多项式系统，如勒让德多项式，生成的插值精度与切比雪夫多项式相似。切比雪夫多项式的主要优势在于，可以使用快速傅里叶变换在基展开中的插值值和系数之间有效地转换。
Remes 算法可以用于生成更准确的近似多项式，但正如 Trefethen 所指出的，对于解析和光滑函数，最佳近似不会比切比雪夫插值带来显著的精度提升。
问题2：自适应稀疏网格是否解决了“维数灾难”？ 我不认为“维数灾难”是思考许多问题时最有用的框架。
自适应稀疏网格在需要实现相同精度的评估点数量方面会比指数增长要好得多；但它们需要比线性增长更多的评估点数量，因此足够高维度的近似问题最终仍会变得过于“昂贵”。
重要的是，自适应稀疏网格是解决许多实际问题的绝佳工具，它们可以处理密集网格无法实现的维数。
结论
作为一种粗略的思维模型，我喜欢将数值问题的解决方案分为三个级别。
第一级：可以表示为整数、有理数值或以已知常数或特殊函数表示的符号解的精确答案。
第二级：如此接近实际值的答案，以至于在实际应用中可以将其视为精确答案。这包括调用标准数值函数或应用BLAS-LAPACK操作。
第三级：答案足够接近有用，但与实际解决方案有显著差异。这可能包括通过正态分布或许多不确定抽样算法近似具有有限方差的独立同分布随机变量的总和。
特别是在统计学中，我觉得对第一级和第三级方法的强调过多。一个问题通常不会被认为是“解决了”，除非它提供了第一级解决方案；并且经常做出不合理的假设，例如共轭先验，以便答案可以整齐地以第一级形式表示。然而，对于最终用户来说，提供第二级解决方案的软件与获得第一级答案一样好；正如我们在假设检验示例中看到的，使用现代数值技术，如高精度求积规则、稀疏网格和切比雪夫插值，我们可以为许多不适合精确解决方案的常见统计问题提供确定性、高效的第二级解决方案。
参考：
[1]: Trefethen, Lloyd. (2011). Six myths of polynomial interpolation and quadrature. Mathematics Today. 47.
[2]: Trefethen Lloyd N. Approximation Theory and Approximation Practice, Extended Edition. USA: Society for Industrial and Applied Mathematics, 2019.
[3]: Kershaw, D. (1971). A Note on the Convergence of Interpolatory Cubic Splines. SIAM Journal on Numerical Analysis, 8(1), 67–74. http://www.jstor.org/stable/2949522
[4]: S. A. Smolyak. Quadrature and interpolation formulas for tensor products of certain classes of functions. Dokl. Akad. Nauk SSSR, 148:1042–1043, 1963. Russian, Engl. Transl.: Soviet Math. Dokl. 4:240–243, 1963.
[5]: Garcke, Jochen (2012). Sparse Grids in a Nutshell. In Garcke, Jochen; Griebel, Michael (eds.). Sparse Grids and Applications. Springer. pp. 57–80. ISBN 978–3–642–31702–6.
[6]: Gerstner Thomas, Griebel Michael. Dimension–Adaptive Tensor–Product Quadrature. Computing. 09 2003. 71. 65–87
[7]: Jakeman John D., Roberts Stephen G. Local and Dimension Adaptive Sparse Grid Interpolation and Quadrature. 2011.
[8]: Klimke Andreas. Uncertainty Modeling using Fuzzy Arithmetic and Sparse Grids. 01 2006. 40–41.
[9]: Berger, J. and J. Mortera (1999). Default bayes factors for nonnested hypothesis testing. Journal of the American Statistical Association 94 (446), 542–554. postscript
[10]: Takahasi, Hidetosi; Mori, Masatake (1974), Double Exponential Formulas for Numerical Integration, Publications of the Research Institute for Mathematical Sciences, 9 (3): 721–741
[11]: Boyd, John & Marilyn, To & Eliot, Paraphrasing. (2000). Chebyshev and Fourier Spectral Methods.
[12]: Aurentz and Trefethen, Chopping a Chebyshev series, Trans. Math. Softw., 2017
