计算方法（上）

Posted on 2022-04-07 Edited on 2022-06-12 In notes

Symbols count in article: 7.9k Reading time ≈ 7 mins.

一些杂记

这里的笔记是对于上课内容的一个浓缩

Part 1

不动点迭代

方法

如果$g(r) = r$，那么称$r$为$g$的一个不动点

我们称下面的方法为不动点迭代法：

猜测一个解$x_0$，之后不断的进行迭代：$x_{i+1} = g(x_i)$

性质：对连续函数$g$，如果$\lim x_i$收敛，那么其一定收敛于不动点

注意到$g(r) = g(\lim x_i) = \lim g(x_i) = \lim x_{i+1} = r$

迭代的速率

对于这种迭代方法，记$e_i$为第$i$次迭代的误差，即$e_i = |x_i - r|$

如果$S = \lim e_{i+1} / e_i < 1$，那么说该不动点迭代线性收敛，$S$为其收敛速度

如果$g$连续可导，$r$为不动点，$S = |g'(r)| < 1$，那么对于足够接近$r$的猜测，不动点迭代法可以收敛到$r$，收敛速度为$S$

证明：$(x_{i+1} - r) = g'(\xi_i) (x_i - r)$，当$g'(x) \leq c< 1$时，$x$离$r$的距离一定在减小

误差

如果计算$f$时有误差函数$\epsilon g$，那么实际上我们在解方程$f(r') + \epsilon g(r') = 0$

求导并整理，可以得到$\Delta r = r'-r \approx - \epsilon \frac{g(r)}{f'(r)}$

Part 2

Lagrange Interpolation

方法

给定$n$个点$(x_1,y_1),...,(x_n,y_n)$，求一个过这些点的多项式

根据代数基本定理，恰好过这$n$个点的，次数低于$n$的多项式至多有一个

我们可以给出其存在性：

定义$L_k = A \prod_{i \neq k}(x - x_i)$，其中$A = \prod_{i \neq k} (x_k - x_i)^{-1}$

观察到$L_k(x_k) = 1$，而对于其余的$i \neq k$，有$ L_k(x_i) = 0$

令$L(x) = \sum y_i L_i(x)$就得到了我们希望的多项式

误差

对于插值出来的多项式$P$，及待插值的$n$次可导函数$f$，有

$f(x) - P(x) = \frac{1}{n!} \prod (x- x_i) f^{(n)} (c)$

Bernstein Polynomial

方法

设$B_{n, i} (x) = \binom{n}{i} x^i (1-x)^{n-i}$，并且$B_n(x) = \sum_{i=0}^n f(\frac{i}{n}) B_{n, i}(x), 0 \leq x \leq 1$

那么，随着$n$增大，$B_n(x)$一致地趋近于$f(x)$，或者说$\lim_{n \to \infty} ||B_n(x) - f(x)||_{\infty} = 0$

分析

设$K_n(x) \sim B(n, x)$，那么$Pr(K_n = i) = B_{n,i}(x)$

因此，$E[f(K_n/n)] = \sum_{i} f(\frac{i}{n}) Pr(K_n = i) = \sum_i f(\frac{i}{n}) B_{n,i}(x) = B_n(x)$

我们只需要尝试证明$E[f(K_n / n)]$一致地趋近于$f(x)$，这依赖于下面几个事实

$\lim_{n \to \infty} Pr(|K_n(x) / n - x| > \delta) = 0$ 一致成立（大数定理）
$\lim_{n \to \infty} Pr(|f(K_n(x) / n) - f(x)| > \delta) = 0$ 一致成立（由1及$f$一致连续）
$\lim_{n \to \infty} E[|f(K_n(x) / n) - f(x)|] = 0$一致成立（由2）

更仔细地分析有：

$|f(x) - E[f(K_n(x)/n)]| \leq ||f||_{\infty} / 2n \epsilon^2 + \sup_{|x-y| \leq \epsilon} |f(x) - f(y)|$

逼近定理

$Weierstrass$定理：对于连续函数$f$，多项式可以任意地逼近

使用Bernstein多项式即可证明，但注意Bernstein并不是最佳的一致逼近多项式

Chebyshev Polynomial

引入

现在先不妨假设函数定义在$[-1,1]$上，我们希望选择$x_1, ..., x_n$，使得 \[ \max_{x \in [-1,1]} (x-x_1)(x-x_2)...(x-x_n) \] 最小

这个问题已经被充分的解决了

答案是选择$x_i = \cos (2i - 1) \pi / 2n$时，上面的问题得到最小值$2^{-(n-1)}$

Chebyshev Polynomial

对于上面这个问题，我们选择上述点之后，实际上得到了一个多项式，我们称之为Chebyshev多项式，定义为 \[ T_n(x) = 2^{n-1}(x-x_1)(x-x_2)...(x-x_n), x_i = \cos (2i - 1) \pi / 2n \] 这个多项式有一个等价的定义，即$T_n(x) = \cos (n \arccos x)$（考虑多项式的次数以及零点）

考虑这两个等价的定义（主要是考虑三角函数的定义），Chebyshev多项式将满足一系列的性质：

$|T_n(x)| \leq 1$

$T_n(x)$的零点恰好为$x_1, ..., x_n$

$T_n(x)$是一个单增区间和单减区间交替的函数，并且交替$n-1$次，在$\cos i \pi / n$处进行交替，每次交替都会发生在$\pm 1$上

$T_{n+1}(x) = 2x T_n(x) - T_{n-1}(x)$

其中的性质3可以用于证明开头提到的问题的结论，也被称为Chebyshev TH

（不妨取$g(x)$为使函数最小的多项式，为方便和$T_n(x)$比较，设$G(x) = 2^{n-1}g(x)$，那么$|G(x)|<1$恒成立，在$T_n(x)$的$n+1$个极值点（取$\pm 1$）上，$T_n(x) - G(x)$和$T_n(x)$有相同的正负，从而$T_n(x) - G(x)$有$n$个零点，这不可能）

Chebyshev 定理的另一种表述为：$T_n(x)$是首一多项式中$||T_n(x)||_{\infty}$最小的多项式

Chebyshev 多项式是多项式的一组正交基，在如下的内积下： \[ \langle f, g\rangle = \int_{-1}^1 f(x)g(x) \frac{d x}{\sqrt{1-x^2}} \]

Fourier Transform

省略

Part 3

Least Square Method

一般的，我们希望解决找到一个点$x^*$，$x^* \in \arg \min_x ||Ax - b||_2^2$

此时，对该函数求导，$\nabla||Ax-b||_2^2 = 2A^T A x - 2 A^T b$，一般就可以求出$x$

Linear system of equation

一般的线性方程组$Ax = b$，我们可以使用高斯消元法来计算，大约需要$O(n^3)$运算

对于线性方程组$Ax = b$，记$e_b$为$b$的相对误差，$e_x$为$x$的相对误差，我们称其条件数$cond(A) = e_x / e_b$

在矩阵$A$可逆的情况下，$e_x = A^{-1} e_b$，那么 \[ cond(A) = \max_{e_b, b \neq 0} \frac{||A^{-1} e||/||A^{-1}b||}{||e||/||b||} = \max_{e \neq 0} \frac{||A^{-1} e||}{||e||} \max_{b \neq 0} \frac{||b||}{||A^{-1}b||} = ||A^{-1}|| * ||A|| \geq 1 \] 其中$||\cdot||$为任意的向量范数和其导出的矩阵范数

Norm

一般的，称由向量范数$||\cdot||$导出的矩阵范数$||A|| := \max_{||x||=1} ||Ax||$为（$||\cdot||$导出的）算子范数

算子范数有如下的性质：

相容性：$||AB|| \leq ||A|| * ||B||$

三角不等式：$||A||-||B||\leq ||A + B|| \leq ||A|| + ||B||$

以下是几个常见的矩阵算子范数：

$||A||_1 = \max_{j} \sum_{i=1}^n |a_{i,j}|$，向量1-范数的导出范数，也称为列范数

$||A||_{\infty} = \max_{i} \sum_{j=1}^n |a_{i,j}|$，向量$\infty-$范数的导出范数，也称为行范数

$||A||_2 = \max_{||x||_2 \leq 1} \sqrt{x^T A^T A x} = \sqrt{\lambda_{max}(A^TA)}$，向量2-范数的导出范数

定义谱半径$\rho(A):= \max\{ |\lambda_1|, ..., |\lambda_n|\}$

$\rho(A) < 1$当且仅当$\lim_{k \to \infty} A^k =0$

对于实对称矩阵而言，$\rho(A) = ||A||_2$

Jacobi method & Gauss-Seide method

Jacobi method

对于线性方程组$Ax = b$

我们把$A$写为$L+D+U$的形式，其中$L, D, U$分别为下三角，对角，上三角阵

此时有$Dx = b - (L+U)x$，因此我们进行迭代：$x_{k+1} = D^{-1}(b - (L+U)x_k)$

Gauss-Seidel method

对于上式的$(L+D+U)x = b$，有$(L+D)x = b -Ux$

我们进行迭代：$x_{k+1} = D^{-1}(b - Ux_k - Lx_{k+1})$

注意到$L$是下三角阵，因此不会产生计算上的顺序问题

分析

我们有如下的定理

如果$A$是严格对角占优的，那么$A$是可逆的，并且对于任何向量$b$和初始猜测$x_0$，线性方程组$Ax = b$使用Jacobi迭代 / Gauss Seidel迭代会收敛到唯一解（即方程组的解）

为了证明上述定理，我们考虑一般性的迭代$x_{k+1} = Ax_k + b$

设$x^*$是该迭代的一个不动点，那么有$x_{k+1}-x^* = A(x_k - x^*) = A^k(x_0 -x^*)$

至少当$A^k \to 0$时，不动点迭代是收敛的，而$A^k \to 0$和$\rho(A) < 1$是等价的

因此，为了证明不动点迭代是收敛的，我们可以尝试证明$\rho(A) < 1$

也就是说，我们希望证明$\rho(D^{-1}(L+U)) < 1$以及$\rho((L+D)^{-1} U) < 1$

对$\rho(D^{-1}(L+U))$考虑，设$\lambda$为$D^{-1}(L+U)$的任意一个特征值，$v$为其特征向量

那么有$(L+U)v = \lambda Dv$，设$||v||_{\infty} = 1$，且$v_m$最大，那么根据$(L+U) v = \lambda Dv$，考虑两边的第$m$个分量，得到$|\lambda a_m| = \sum_{j \neq m} |a_{mj}| < |a_m|$，从而$|\lambda| < 1$，进而有$\rho(D^{-1}(L+U)) < 1$，另一个分析是类似的

Moore-Penrose Pseudoinverse

考虑方程$Ax = b$

当$A \in \mathbb{R}^{n \times m}(n > m)$时，并且保证$A$满秩时，我们将方程$Ax = b$改写为$A^T A x = A^T b$

此时，$A^T A$仍然满秩

那么，我们可以写出$x = (A^TA)^{-1} A^Tb$，这是一个pseudoinverse

Richardson iteration

引入

根据$Cayley-Hamilton$定理，存在$n$次多项式$p(A)$使得$p(A)A = I$

这启发我们从多项式角度来迭代地求出$A^{-1}$

方法

对于线性方程组$Ax = b$

令$x_0 = 0$，$x_{k+1} = (I - \alpha A) x_k + \alpha b$

分析

当$\rho(I - \alpha A) < 1$时收敛，由于$I, \alpha A$可交换，因此实际上等价于$|1 - \alpha \lambda_{max}| < 1$

下面我们来探讨为什么这等价于研究一个多项式

在线性方程组$Ax = b$中，实际上我们在寻找$x = A^{-1}b$，即$x = p(A)b$

寻找$p(x)$，等价于寻找$q(x) = 1-xp(x)$

假设$x_k = p_k(A)b$，那么$x_{k+1} = ((I-\alpha A)p_k(A) + \alpha)b$

相当于$p_{k+1}(x) = p_k(x)(1-\alpha x) + \alpha$，即$1-q_{k+1}(x) = (1-q_k(x))(1 - \alpha x) + \alpha x$

得到$q_{k+1}(x) = q_k(x)(1-\alpha x)$

通过归纳，可以证明我们在考虑多项式$p_k(x) = (1-(1 - \alpha x)^k)/x$

Chebyshev iteration

方法

$x_{k+1} = (I - \alpha_kA) x_k + \alpha_k b$

分析

定义$e_k = x_k - x^*$，那么$e_k = \prod(I-\alpha_i A) e_0$

因此$||e_k|| \leq ||\prod_i (I - \alpha_i A)|| *||e_0||$

我们希望关于$A$的那一项系数尽可能小，此时运用Chebyshev多项式就有不错的效果

Conjugate gradients

Krylov 子空间

记Krylov子空间为：$K_0 = \{0\}, K_i = span\{b, Ab, ..., A^{i-1} b\}$

设$x_i = \arg \min_{x \in K_i} ||x - x_*||_A^2$，其中$x_*$满足$Ax_* = b$，而$||x||_A = \langle x, x \rangle_A = x Ax^T$，是由正定矩阵$A$定义出的内积导出的范数

引理

如果$\langle x, y \rangle_A = 0$，那么我们称$x, y$关于$A$共轭

设$v_i = x_i - x_{i-1}$，那么$v_1, ..., v_n$两两共轭

注意到对$i < j$，$\nabla ||x_j - x_*||_A^2 = 2 (Ax_j - b)$与$K_j$正交（$x_j$最小），从而与$K_i$正交

类似的$A x_{j-1} - b$与$K_i$也正交，因此$Av_j = (Ax_j - b) - (Ax_{j-1} - b)$也与$K_i$正交，从而得到了我们希望的结论

$K_k = span\{v_1, ..., v_k\}$

由上面的结论就可以得出

设$r_i = b - Ax_i$，$K_k = span\{v_1,...,v_{k-1}, r_{k-1}\}$

注意到$r_{i-1} \in K_i$，但是$r_{i-1}$与$K_{i-1}$正交（根据$x_{i-1}$的定义以及梯度为0），因此$K_k = span\{v_1,...,v_{k-1}, r_{k-1}\}$

设$v_i = x_i - x_{i-1}$，那么$v_i$是关于$v_{i-1}, r_{i-1}$的函数

由于$v_i \in K_i$，因此$v_i$可以用$v_1, ..., v_{i-1}, r_{i-1}$线性表出

不妨设$v_i = c_0 r_{i-1} + \sum_{i=1}^{i-1} c_i v_{i-1}$，注意到$r_{i-1}$与$v_1,...,v_{i-1}$正交，$v_i$之间共轭

再注意到$Av_j \in K_{i-1}, \forall j \leq i-2$，因此$r_{i-1}^TAv_j = 0$。即$r_{i-1}$与$v_1,...,v_{i-2}$共轭

考虑$v_i^T r_{i-1} = c_0 ||r_{i-1}||_2^2$，得到$c_0 = v_i^T r_{i-1} / ||r_{i-1}||_2^2$

再考虑$v_i^T A v_j = 0$，我们就能相应地解出$c_1, ..., c_{i-1}$

误差分析

略

Power method

介绍

让我们考虑怎么迭代地求出特征值以及特征向量

注意到，对于可对角矩阵$A$，其存在$n$个线性无关的特征向量，设为$v_1, ..., v_n$，对应的特征值为$\lambda_1, ...$，并且$|\lambda_1| \geq |\lambda_2| \geq ...$

那么，如果$x = \sum \alpha_i v_i$，则$A^k x = \sum \lambda_i^k \alpha_i v_i = \lambda_1^k (\alpha_1 v_1 + \frac{\lambda_2^k}{\lambda_1^k} \alpha_2 v_2 + ...)$

因此，当$k$足够大时，我们就会得到一个右边是一个$\lambda_1$的近似特征向量$v_1'$

通过最小二乘法，最小化$||Av_1' - \lambda_1' v_1'||_2$来得到$\lambda_1'$，不难得到$\lambda_1' = v_1'^TAv_1' / v_1'^Tv_1'$

分析

我们对近似得到的特征值进行分析

对实对称矩阵$A$，假设近似特征向量$x$有$||x||_2 = 1$，并且$\lambda$满足$||Ax - \lambda x||_2 < \epsilon$，那么$\min |\lambda_i - \lambda| < \epsilon$

$||Ax - \lambda x||_2^2 = \sum_{i} |\alpha_i|^2 |\lambda_i - \lambda|^2 \geq \min |\lambda_i - \lambda|^2$

这个定理告诉我们，我们得到的特征值是在所给特征向量下的比较好的近似

Inverse Power method

介绍

注意到$(A-qI)^{-1}$的特征值为$\{\frac{1}{\lambda_i - q}\}$

如果我们选择适合的参数$q$，那么我们就能够得到相应的最大的$\frac{1}{\lambda_i- q}$，进而得到$\lambda_i$

QR decompesition

介绍

$Q_0 = I$，而$Q_i$满足$AQ_{i-1} = Q_i R_i(i \geq 1)$，其中$Q_iR_i$是$AQ_{i-1}$的$QR$分解

这样可以进行特征值的同时迭代