数值分析 笔记 6

函数逼近与函数插值

这两种方法本质上都是在用简单函数去近似复杂函数的过程，不同的是：

• 逼近要求整体上误差最小
• 插值要求在一些点上简单函数与原函数的值相等

函数逼近

函数逼近是指，对于给定函数 $f(x)$，希望在某个函数空间 $\varPhi$ 中找到一个函数 $p(x)$，使得在某种度量意义下，误差函数 $p(x) - f(x)$ 达到最小

通常 $f(x)$ 是一个有表达式的复杂函数或者表格函数，而 $\varPhi$ 通常有：多项式、有理分式、指数函数、三角函数、分段多项式等简单函数类

函数逼近的例子有很多，例如常见的 Taylor 展开、Fourier 变换、数据点拟合等等

我们通常使用函数空间上的范数来度量误差，$C[a, b]$ 上常见的范数包括：

$$\begin{align*} ||f(x)||_{\infty} &= \max\limits_{x\in[a, b]}|f(x)| \\ ||f(x)||_{1} &= \int_{a}^{b}|f(x)|\mathrm{d}x \\ ||f(x)||_{2} &= \Big[\int_{a}^{b}f^{2}(x)\mathrm{d}x\Big]^{1/2} \\ \end{align*}$$

其中 2-范数又称内积范数，定义见下方

在这些范数中，使用 $\infty$-范数来度量误差是最优的，因为他规定了在整个定义域上两个函数都非常接近，这种逼近方式被称为最佳一致逼近

但是最佳一致逼近非常难以求解，而 1-范数实际上是考察两个函数图像之间的面积，误差较大，因此通常使用 2-范数来度量逼近误差，这被称为最佳平方（最小二乘）逼近

函数内积

函数内积定义为：

$$\langle u, v\rangle = \int_a^b u(x)\overline{v(x)}\mathrm{d}x$$

对于内积有 Cauthy-Schwarz 不等式：

$$|\langle u, v\rangle|^{2} \leq \langle u, u\rangle\langle v, v\rangle$$

对于实内积空间 $S$ 中的 $n$ 个函数 $u_1, \dots, u_n$，其线性无关的充要条件是 Gram 矩阵 $\boldsymbol{G}$ 非奇异，Gram 矩阵定义为：

$$\boldsymbol{G}_{ij} = \langle u_i, u_j\rangle$$

由内积的对称性和正定性可以得到，非奇异的 Gram 矩阵是对称正定的

上述内积定义中，两个函数的地位是等价的，也就是内积的对称性，但是有的时候这并不是想要的，内积本质上是一种广义上的求和操作，因此根据加权求和的思想，引入加权内积

权函数 $\rho(x)$ 定义为满足以下条件的函数：

• $\forall x\in [a, b],\quad\rho(x) \geq 0$
• $\forall k\in \mathbb{N},\quad \int_a^b x^{k}\rho(x)\mathrm{d}x$ 是存在的
• 对于任意非负连续函数 $g(x)$，若 $\int_a^b g(x)\rho(x)\mathrm{d}x = 0$，则 $g(x) \equiv 0$，即$\rho(x)$ 没有局部恒为 $0$ 的情况

则定义加权内积为：

$$\langle u, v\rangle = \int_a^b \rho(x)u(x)\overline{v(x)}\mathrm{d}x$$

最佳平方逼近

考虑函数类 $\varPhi = \text{span}\{\varphi_{1}, \dots, \varphi_{n}\}$，则我们需要找到 $S(t) = \sum\limits_{i=1}^{n}x_i\varphi_{i}(t)$，使得 $||S(t) - f(t)||_{2}^{2}$ 最小

根据内积的双线性性与对称性，有：

$$\begin{align*} F &= ||S(t) - f(t)||_{2}^{2} \\ &= \langle\sum\limits_{i=1}^{n}x_i\varphi_{i} - f, \sum\limits_{i=1}^{n}x_i\varphi_{i} - f\rangle \\ &= \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}x_{i}x_{j}\langle\varphi_{i}, \varphi_{j}\rangle - 2\sum\limits_{i=1}^{n}x_i\langle f, \varphi_{i}\rangle + \langle f, f\rangle \end{align*}$$

因此：

$$\begin{align*} \frac{\partial F}{\partial x_{k}} = 2\sum\limits_{i = 1}^{n}x_{i}\langle \varphi_{i}, \varphi_{k}\rangle - 2\langle f, \varphi_{k}\rangle \end{align*}$$

令所有的偏导数等于 $0$，可以得到一个线性方程组 $\boldsymbol{Gx = b}$，其中 $\boldsymbol{G}$ 是 Gram 矩阵，$\boldsymbol{b} = [\langle f, \varphi_{1}\rangle, \dots, \langle f, \varphi_{n}\rangle]^{T}$

如果这些函数线性无关，那么 $\boldsymbol{G}$ 是对称正定的，这个线性方程组有唯一解，并且可以证明这个解可以让误差度量达到最小值，几何意义上来说，这代表：

$$S^{*} - f \perp \Phi$$

多项式最佳平方逼近

可以取：$\Phi = \mathbb{P}_{n - 1}$，即次数不超过 $n - 1$ 的所有多项式的集合

对于 $[0, 1]$ 区间，此时有：

$$\langle\varphi_{i}, \varphi_{j}\rangle = \int_0^1 t^{i + j - 2} \mathrm{d}t = \frac{1}{i + j - 1}$$

即 $\boldsymbol{G} = \boldsymbol{H}_{n}$

根据 Weierstrass 定理，多项式可以对连续函数进行任意好的逼近，但是根据上述结论，其 Gram 矩阵会随着次数增大而极具病态，并且是稠密的，计算复杂度高

如果我们将基函数取成一组正交基，那根据正交的特性很容易得到，Gram 矩阵会变成对角阵！这种情况下解方程是非常容易的

可以通过 Gram-Schmidt 正交化的方式从 $\{1, t, \dots, t^{n - 1}\}$ 得到一组正交多项式，具体来说：

$$\begin{align*} \varphi_{1}(t) &= 1 \\ \varphi_{k}(t) &= t^{k - 1} - \sum\limits_{i = 1}^{k - 1}\frac{\langle t^{k - 1}, \varphi_{i}(t)\rangle}{\langle \varphi_{i}(t), \varphi_{i}(t)\rangle}\varphi_{i}(t) \\ \end{align*}$$

对于正交基，可以很方便地求出逼近结果：

$$S^{*} = \sum\limits_{k = 1}^{n}\frac{\langle f(t), \varphi_{k}(t)\rangle}{\langle \varphi_{k}(t), \varphi_{k}(t)\rangle}\varphi_{k}(t)$$

在选取不同的定义域与不同的权函数时，Gram-Schmidt 正交化的结果可以表示成不同的递推形式，如下表：

表中的递推式可以帮助我们方便的求出 $\varphi_{k}$

如果定义域不满足上表中的条件，可以使用线性映射的方式修改，例如对于 Legendre 多项式，从 $[-1, 1]$ 映射到 $[a, b]$：

$$s = \frac{b - a}{2} t + \frac{b + a}{2}$$

因此映射后的 Legendre 多项式为：

$$\tilde{P_{k}}(s) = P_{k}(\frac{2s - (b + a)}{b - a})$$

注意在求和的时候计算系数与基函数都应该是映射后的

可以证明使用多项式来逼近的时候，逼近误差一致不大于 $\frac{\varepsilon}{\sqrt{n}}$

在 $n$ 充分大的时候度量方式为 2-范数，如果原函数二阶导连续，则可以扩展为 无穷范数

最小二乘法

与最佳平方逼近的思路大致相同，不同的是这种方法通常用于给定数据点对的表格，需要最优化的函数不再是内积范数，而是：

$$\sum\limits_{i=1}^{m}[S(t_{i}) - f(t_{i})]^{2}$$

本质上来说，定义在离散点处的表格函数也构成了一个线性空间，离散情况下，从连续区间上的函数变成了离散的向量，范数也从积分变成了求和，因此最小二乘是一种最佳平方逼近

法方程法

将表格函数定义为向量 $\boldsymbol{f} = [f(t_{1}), \dots, f(t_{m})]^{T}$，则需要找到向量 $\boldsymbol{x}$ 使得 $||\boldsymbol{Ax - f}||_{2}$ 最小，其中 $\boldsymbol{A} \in \mathbb{R}^{m \times n}$ 定义为：

$$\boldsymbol{A}_{ij} = \varphi_{j}(t_{i})$$

$\boldsymbol{A}$ 的每一列代表的是一个基函数，每一行代表的是某一个输入在不同基函数下的函数值，因此可以写成：

$$\boldsymbol{A} = [\boldsymbol{\varphi}_{1}, \dots, \boldsymbol{\varphi}_{n}]$$

考虑到离散函数可以用向量来表示，而向量的范数就是对应下标元素求和的形式，因此：

$$\begin{align*} \boldsymbol{G = A}^{T}\boldsymbol{A} \\ \boldsymbol{b = A}^{T}\boldsymbol{f} \end{align*}$$

也即我们将最小二乘法转换为了最佳平方逼近的 $\boldsymbol{Gx = b}$，如果表格基函数线性无关，也即 $\boldsymbol{A}$ 是列满秩的，那 $\boldsymbol{Gx = b}$ 一定存在唯一解

然而，通常最小二乘的使用场景都是拟合大量数据点，即 $m >> n$ 时，但是这种情况下就会导致 $\boldsymbol{Gx = b}$ 并不等价于 $\boldsymbol{Ax = f}$，甚至在列满秩的情况下很容易后者无解，不过可以证明的是 $\boldsymbol{Gx = b}$ 的解满足使得 $||\boldsymbol{Ax - f}||_{2}$ 最小

因此得到

对于连续函数来说，是否线性无关与输入没有关系，而对于表格函数来说，是否线性无关是与观测点的值有关的，并且观测点可能出现相同的（即一对多），有 Haar 条件：

若 $t_{i}$ 中至少有 $n$ 个不同的值，则 $\mathbb{P}_{n - 1}$ 对应的表格基函数是线性无关的

正交变换法

由于法方程法会面临最佳平方逼近的同样问题，因此我们仍然需要选取正交的表格基函数，一种可行的方法是 $\boldsymbol{A}$ 的 QR 变换：

由于正交变换不改变 2-范数，因此：

$$||\boldsymbol{f - Ax}||_{2} = ||\boldsymbol{Q}^{T}\boldsymbol{f - Rx}||_{2}$$

对于 $\boldsymbol{A} \in \mathbb{R}^{m\times n}$，其对应的 QR 分解可以分块为：

$$\boldsymbol{Q} = \begin{bmatrix} \boldsymbol{Q}_{1}, \boldsymbol{Q}_{2} \end{bmatrix}, \quad \boldsymbol{R} = \begin{bmatrix} \boldsymbol{R}_{1} \\ \boldsymbol{O} \end{bmatrix}$$

其中 $\boldsymbol{Q}_{1} \in \mathbb{R}^{m \times n}$，$\boldsymbol{R}_{1} \in \mathbb{R}^{n \times n}$

则有：

$$\begin{align*} ||\boldsymbol{Q}^{T}\boldsymbol{f - Rx}||_{2} &= ||\boldsymbol{Q}_{1}^{T}\boldsymbol{f} - \boldsymbol{R}_{1}\boldsymbol{x}||_{2} + ||\boldsymbol{Q}_{2}^{T}\boldsymbol{f}||_{2} \\ &\geq ||\boldsymbol{Q}_{2}^{T}\boldsymbol{f}||_{2} \end{align*}$$

当且仅当 $\boldsymbol{Q}_{1}^{T}\boldsymbol{f} = \boldsymbol{R}_{1}\boldsymbol{x}$ 时取等

因此得到

由于 $\boldsymbol{R}$ 是三角阵，因此解方程快速且稳定

函数插值

本质上是一种假设数据没有误差的时候做的函数拟合，我们用 $P(x)$ 代表插值结果，则需要满足 $P(x_i) = y_i$

常见的插值函数包括多项式、三角函数、有理分式、分段函数等

多项式插值

对于有 $n + 1$ 个数据点的插值条件，我们可以使用 $n$ 次多项式来进行插值，有代数基本定理可以知道，一定可以找到唯一的一个多项式满足条件

最直接的方式是求解线性方程组：

$$\begin{cases} a_0 + a_1 x_0 &+ \dots + a_n x_0^n = y_0 \\ a_0 + a_1 x_1 &+ \dots + a_n x_1^n = y_1 \\ &\vdots \\ a_0 + a_1 x_n &+ \dots + a_n x_n^n = y_n \\ \end{cases}$$

系数矩阵为 Vandermonde 矩阵：

$$\boldsymbol{A} = \begin{bmatrix} 1 & x_0 & \cdots & x_0^n \\[0.5em] 1 & x_1 & \cdots & x_1^n \\[0.5em] \vdots & \vdots & \ddots & \vdots \\[0.5em] 1 & x_n & \cdots & x_n^n \\[0.5em] \end{bmatrix}$$

但是这个矩阵是稠密的，并且使用这种方法不利于分析，因此引入 Lagrange 插值法

Lagrange 插值

考虑 $n = 1$ 的情况，也就是使用一维方程（直线）去插值两个节点，我们可以使用经典的两点式直线方程得到结果：

$$y = \frac{x - x_1}{x_0 - x_1}y_0 + \frac{x - x_0}{x_1 - x_0}y_1$$

这实际上可以看成是一些基函数的线性组合，系数为给定的函数值，即：

$$y = y_0l_0(x) + y_1l_1(x)$$

其中 $l_0(x)$ 与 $l_1(x)$ 都是一次函数，并且满足 $l_0(x_0) = l_1(x_1) = 1, l_0(x_1) = l_1(x_0) = 0$

这就是最基础的 Lagrange 插值

扩展到更高维的形式，则插值函数为：

$$L_n(x) = \sum\limits_{k = 0}^{n}y_kl_k(x)$$

基函数 $l_k(x)$ 的次数也为 $n$，并且满足：

$$l_k(x_i) = \begin{cases} 1 & i = k \\ 0 & i \neq k \end{cases}$$

确定基函数最简单的方式是利用给定的 $n$ 个零点，这可以得到一个成比例的函数簇，之后根据 $x_k$ 处的函数值求出系数即可

我们定义 $R_{n}(x) \equiv f(x) - P_{n}(x)$ 为插值余项，则 Lagrange 插值的余项为：

$$R_{n}(x) = \frac{f^{(n + 1)}(\xi(x))}{(n + 1)!}\omega_{n + 1}(x)$$

其中 $\xi(x) \in (a, b)$

但是 Lagrange 插值有一个巨大的问题，即在插值点有变动的时候，插值函数的变动非常巨大，开销太大，因此引入 Newton 插值

Newton 插值

Newton 插值的主要思想是逐步增加插值点的数目，逐步增加高次项，而不是每一个插值点都使用同样次数的基函数

扩展方式为：

$$P_n(x) = P_{n - 1}(x) + c_{n}\omega_n(x)$$

最终公式为：

$$N_{n}(x) = \sum\limits_{i = 0}^{n}c_{i}\omega_i(x)$$

为了求出系数 $c_i$，我们引入函数差商的概念，函数的 $k$ 阶差商是如下递归定义出来的：

$$\begin{align*} f[x_{0}] &= f(x_{0}) \\ f[x_{0}, x_{1}, \dots, x_{k}] &= \frac{f[x_0, \dots, x_{k - 2}, x_k] - f[x_0, x_1, \dots, x_{k - 1}]}{x_k - x_{k - 1}} \end{align*}$$

可以归纳证明出：

$$f[x_0, \dots, x_k] = \sum\limits_{i = 0}^{k}\frac{f(x_i)}{\omega^{'}_{k + 1}(x_i)}$$

直接取 $c_k = f[x_0, \dots, x_k]$ 即可，此时余项为：

$$R_{n}(x) = f[x, x_0, \dots, x_n]\omega_{n + 1}(x)$$

根据该式以及余项的定义，可以证明差商的性质：

在使用 Newton 插值的过程中，可以使用 $N_{k + 1}(x) - N_{k}(x)$ 或者余项来判断是否需要使用新的插值点

分段多项式插值

对于上述的多项式插值有一些很显著的问题：

1. 数值稳定性差，容易被异常数据影响
2. 无法保证曲线的单调性和凸性
3. 容易出现过拟合（Runge 现象：$n$越高对曲线的拟合效果不一定越好）

可以采用分段多项式插值的方式来修改

分段线性插值

定义 $x_{-1} = x_0 < x_{1} < \cdots < x_n = x_{n + 1}$

则基函数定义为：

$$l_i(x) = \begin{cases} \dfrac{x - x_{i - 1}}{x_i - x_{i - 1}} & x_{i - 1} \leq x \leq x_{i} \\[1em] \dfrac{x - x_{i + 1}}{x_i - x_{i + 1}} & x_{i} \leq x \leq x_{i + 1} \\[1em] 0 & \text{others} \end{cases}$$

也即每两个数据点之间是线性的

令

$$I_h(x) = \sum\limits_{i =0}^{n}y_il_i(x)$$

为分段线性插值结果，则每个小区间上的余项都是一个二阶的 Lagrange 余项，则可以证明：

分段 Hermite 插值

Hermite 插值不仅保证每个点处的函数值相同，还保证了每个点处的导数值也相同，即插值目标为：

$$\begin{align*} H(x_i) &= f(x_i) \\ H'(x_i) &= f'(x_i) \end{align*}$$

可以证明 Hermite 多项式一定有解且唯一，结果为：

$$H(x) = \sum\limits_{j=0}^{n}[f_j\alpha_j(x) + f'_j\beta_j(x)]$$

其中 $\alpha_j, \beta_j$ 是基函数，满足：

$$\begin{align*} \alpha_j'(x_i) = \beta_j(x_i) &= 0 \\ \alpha_j(x_i) = \beta_j'(x_i) &= \begin{cases} 1 & i = j \\ 0 & i \neq j \end{cases} \end{align*}$$

他们可以从 Lagrange 插值的基函数推出来：

$$\begin{align*} \alpha_j(x) &= [1 - 2(x - x_j)l'_j(x_j)]l_j^2(x) \\ \beta_j(x) &= (x - x_j)l_j^2(x) \end{align*}$$

如果是特殊的 Hermite 插值，例如只用保证在一个点处的导数值相同，则可以用 Newton 或 Lagrange 插值之后，再在后面添加一项 $C\omega_{n + 1}(x)$ ，求出常数 $C$ 即可

当 $n = 1$ 的时候，Hermite 实际上是在用一条三次曲线去插两个点，这杯称为两点三次插值，利用这种方法可以引出分段 Hermite 插值，即每一段都用两点三次的方法去插，公式很简单，只需要把分段线性插值的基函数代入到上述 $\alpha_j, \beta_j$ 的表达式中即可

保形分段插值

在 Hermite 中我们需要知道插值点处的导数值，但是这在很多情况下是几乎不可能的，例如一组统计数据基本上无法得知其导数是如何的，因此这种情况下我们需要通过某种方法来认为的“设置”导数值，进而得到一个光滑并且保凸的曲线

一种常见的方式是利用数据点与邻近数据点的割线来估计该数据点处的导数，具体来说对于 $x_k$ 处的导数：

首先计算两侧割线斜率：

$$\begin{align*} d_{k - 1} = \frac{f(x_k) - f(x_{k - 1})}{x_k - x_{k - 1}} \\ d_{k} = \frac{f(x_{k + 1}) - f(x_{k})}{x_{k + 1} - x_{k}} \\ \end{align*}$$

如果 $d_{k - 1}d_{k} \leq 0$ 即两侧单调性不同，则认为 $f'_k = 0$，反之对两侧的割线斜率用区间长度做加权平均：

$$\frac{w_{k - 1} + w_{k}}{f'_k} = \frac{w_{k - 1}}{d_{k - 1}} + \frac{w_{k}}{d_{k}}$$

其中：

$$\begin{align*} w_{k - 1} = h_{k - 1} + 2h_k&\quad w_{k} = 2h_{k - 1} +h_{k} \\[1em] h_{k - 1} = x_{k} - x_{k - 1}&\quad h_{k} = x_{k + 1} - x_{k} \end{align*}$$

而端点处的函数值可以做单侧分析得到

样条插值

上述的分段插值有一个显而易见的问题，它只保证了每个插值点处的一阶导数连续，更高阶的导数则无法保证，而我们知道在物理上，如果这条插值曲线的势能达到最小，那曲线一定是二阶导数连续的

满足在整体区间上二阶导函数连续的插值函数称为样条插值函数，如果每个小区间上都是三次函数，则称为三次样条插值函数

每个区间上 $s_j(x)$ 是三次多项式，需要 $4$ 个方程才能唯一确定系数，因此一共需要 $4n$ 个方程，而良定义的方程为：

$$\begin{align*} s_j(x_j) &= f_j \quad j = 0, \dots, n - 1 \\ s_j(x_{j + 1}) &= f_{j + 1} \quad j = 0, \dots, n - 1 \\ s_{j - 1}'(x_j) &= s_j'(x_j) \quad j = 1, \dots, n - 1 \\ s_{j - 1}''(x_j) &= s_j''(x_j) \quad j = 1, \dots, n - 1 \\ \end{align*}$$

一共 $2n + 2(n - 1) = 4n - 2$ 个，因此还需要从端点处设置两个方程才够，有很多不同的设置方法，包括但不限于：

1. 设置端点处的一阶导数值
2. 设置端点处的二阶导数值（如果端点处的二阶导数为 $0$ 则称为自然样条插值）
3. 假设 $x_n - x_0$ 是函数的周期
4. 给定第一个区间和最后一个区间的表达式

在给定足量的方程之后，就需要想办法确定处插值函数了，直接解方程比较复杂，有两种可行的方式：

1. 假设每个插值点的一阶导数已知，于是可以用分段 Hermite 插值法，之后再根据二次导数连续的特性来反解出参数
2. 假设每个插值点的二阶导数已知，先利用插值条件解出表达式，之后反解出参数