您可以自由地共享和演绎本文稿, 只需遵守开源协议[CC By 4.0]
如果这篇文章帮助到你, 可以请我喝一杯咖啡~
矩阵特征值的计算方式
幂法
设置方阵$A$的各特征值$\lambda$中仅有唯一的模长最大的值，且各特征矢$u_i$线性无关，即：
\[\vert \lambda_1\vert > \vert \lambda_2 \vert \geqslant \vert \lambda_3\vert \geqslant ... \geqslant \vert \lambda_n\vert\]设$x^{(0)}$是任意不为零的$n$维向量：
由迭代公式：
\[x^{(k+1)} = Ax^{(k)} \tag{4.1}\]可以证明，在$k$足够大时，矩阵特征值$\lambda$是向量$x^{(k+1)}$与向量$x^k$的任意两个对应分量之商$x^{(k+1)}_i/x^{k}_i$
证明：
由于$u_i$是$n$个$n$维向量且线性无关，故它是$n\times1$空间中的一组基底，故无论$x^{(0)}$的取值，它总能被表示为$x^{(0)}= \sum_{i=1}^na_iu_i$，其中，$a_i$不全为0. 即$x^{(0)}$是各个特征矢的线性组合.
则由迭代公式：
\[x^{(k+1)} = Ax^{(k)} = A^{k+1}x^{(0)}=\sum_{i=1}^nA^{k+1}a_iu_i=\sum_{i=1}^n\lambda_i^{k+1}a_iu_i\]重写上述等式：
\[x^{(k+1)}=\lambda_1^{k+1}[a_1u_1+(\frac{\lambda_2}{\lambda_1})^{k+1}a_2u_2+(\frac{\lambda_3}{\lambda_1})^{k+1}a_3u_3+...+(\frac{\lambda_n}{\lambda_1})^{k+1}a_nu_n] \tag{4.2}\]不妨假设$\lambda_1$是各个特征值中绝对值最大的一个，则在$k$充分大时，等式(4.2)可以截断为：
\[x^{(k+1)}\approxeq\lambda_1^{k+1}a_1u_1 \tag{4.3}\]此时，$x^{(k+1)}$可以看作是方阵$A$对应于特征值$\lambda_1$的一个特性矢. 同时，在$k$充分大时，
\[x^{(k)}\approxeq\lambda_1^{k}a_1u_1 \tag{4.4}\]比较式子(4.3)和(4.4)，可以知道，$x^{(k+1)}$和$x^{(k)}$线性相关，且线性系数之比为$\lambda_1$
故，可以将$x^{(k+1)}$和$x^{(k)}$的任一对应分量$x_i$作比，有：
\[\lambda_1 = x^{(k+1)}/x^{(k)}\tag{4.5}\]$\blacksquare$
但是精妙的数学逻辑应用到计算机上就会变得粗糙，由于我们需要重复计算若干多次式(4.1)，向量$x^{(k)}$中的分量可能出现溢出或舍入为0的情况. 故这时我们需要对每一次迭代出的向量$x^{(k)}$进行归一化处理：
\[\left\{ \begin{aligned} &y^{(k)} = \frac{x^{(k)}}{\Vert x^{(k)} \Vert_\infin}\\ &x^{(k+1)} = Ay^{(k)} \\ \end{aligned} \right.\]由上述原理，在实际计算中，在式子(4.5)的分量选择过程中，我们应当选择绝对值最大的分量以消除舍入误差.
有效性讨论
- 若不满足”仅有一个模取得最大值”的情况，由上述推导过程易知算法仍有效
- 特殊地，若有$\lambda_1 = -\lambda_2$时，迭代过程变为一摆动序列，则比较$k$足够大时的任意相邻三项即可得到求出$\lambda_1$和$u_1$与$u_2$
- 易知，算法的收敛速度取决于$\vert\lambda_2\vert/\vert \lambda_1\vert$
迭代加速
原点移位法
设$\lambda_i$是$A$的特征值，则$\lambda_i - \lambda_0$是$A-\lambda_0E$的特征值.
因此我们只需要找到合适的$\lambda_0$使得$\forall i \in \mathbb{N}^*, \vert \lambda_1 - \lambda_0\vert > \vert \lambda_i - \lambda_0\vert$且$max\frac{\vert \lambda_i - \lambda_0\vert}{\vert \lambda_1 - \lambda_0\vert}<\vert \frac{\lambda_2}{\lambda_1}\vert$，我们就可以用矩阵$A-\lambda_0E$来代替矩阵$A$，从而获得更快的收敛过程.
但实践中我们不可能先验地知晓$A$的所有特征值，故此该方法在实际上难以精确地执行.
Aitken方法
我们称幂法中”$k$足够大时，$x^{(k+1)}$和$x^{(k)}$成比例”为一阶收敛. 若序列$a_k$一阶收敛至$a$，即：
\[\lim_{k\to \infin}\frac{a_{k+1}-a}{a_k-a} = c \not ={0}\]作近似：在$k$充分大时有：
\[\frac{a_{k+2}-a}{a_{k+1}-a}\approx \frac{a_{k+1}-a}{a_{k}-a}\]解出$a$的表达式：
\[a' \approx a_k-\frac{(a_{k+1} - a_k)^2}{a_{k+2}-2a_{k+1}+a_k}\]每次进行幂法迭代之后，我们都依据上式计算一次$\lambda’$，若$\lambda’$与上一个$\lambda’$之间的误差满足条件，则输出结果.
函数插值法
已知平面上若干互不相同的点$(x_0, y_0), (x_1, y_1), (x_2, y_2), …, (x_n, y_n)$，其关系属于某未知对应关系$f(x)$，如何构造一函数，使得该函数可以穿过这些点并刻画$f(x)$？
这种问题被称为插值问题，上述的函数被称为插值函数.
Lagrange多项式法
Lagrange多项式法使用一个$n$次多项式来拟合上述的$(n+1)$个点，且可以得出唯一的结果.
证明：
对于点集${(x_0, y_0), (x_1, y_1), (x_2, y_2), …, (x_n, y_n)}$，找到一$n$次多项式$f(x)=a_0+a_1x+a_2x^2+…+a_nx^n$使得$\forall x_i, f(x_i) = y_i$等价于方程组：
\[\left\{ \begin{aligned} & f(x_0) = y_0\\ & f(x_1) = y_1\\ & ...\\ & f(x_n) = y_n\\ \end{aligned} \right.\tag{4.6}\]展开之，可使其变为一个系数矩阵是范德蒙矩阵的关于多项式系数向量$A$的线性方程组，且范德蒙矩阵在$x_i$互不相同时满秩.
多项式系数向量$A$存在且唯一.
$\blacksquare$
我们对Lagrange多项式（虽然还没有完全提出它的求解方式）进行误差估计：
设置有若干节点$(x_i, y_i), i\in [0, n]$是$[a, b]$上的$(n+1)$个互异节点，若待拟合函数在$(a,b)$上$(n+1)$阶可导，即$f(x)\in \mathbb{C}^{(n+1)}[a,b]$，若我们用Lagrange多项式$\phi(x)$拟合之，则其插值余项（截断误差）可以表述为：
\[R_n(x) = f(x) - \phi(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x)\tag{4.7}\]其中$\omega_{n+1}(x)=\Pi_{i=0}^n(x - x_i)$，$\xi \in (a,b)$
直球使用微分中值定理即可证明：
由定义，$R(x) = f(x) - \phi(x)$，再由Lagrange多项式的基本假设（式4.6），有$\forall i, R(x_i) = 0$.
故可以等价地将$R(x)$转写为$R(x) = K(x)\omega_{n+1}(x)$，其中$\omega_{n+1}(x)=\Pi_{i=0}^n(x - x_i)$，$K(x)$是某未知的多项式，我们需要找到它：
使用惯用的套路构造辅助函数：
\[h(t) = R(t)-K(x)\omega_{n+1}(t)\]观察到上述的函数有$(n+2)$个零点：$x, x_0, x_1, x_2, …, x_n$，而我们事先约定这些零点都在$f(x)$的解析域$(a,b)$上（否则Lagrange多项式将失去意义）
对自变量$t$反复应用罗尔中值定理$(n+1)$次，则存在某一点$\xi \in (a,b)$，使得：
\[h^{(n+1)}(\xi) = R^{n+1}(\xi)-K(x)(n+1)! = 0\]故得到$K(x) = \frac{R^{n+1}(\xi)}{(n+1)!}$，$\xi \in (a,b)$
而$R(x) = f(x) - \phi(x)$，$\phi(x)$是一个$n$次多项式，在$(n+1)$次求导后变为0，故
\[K(x) = \frac{f^{n+1}(\xi)}{(n+1)!}, \xi \in (a,b)\]$\blacksquare$
式4.7的意义在于，若记$M = \max_{(a,b)}\vert f^{(n+1)}(x) \vert$，则Lagrange插值法的截断误差$R(x) \leqslant \frac{M}{(n+1)!}\omega_{n+1}(x)$，这体现了对误差的一种控制.
Lagrange多项式的求解
最直观地，我们可以直接通过式$(4.6)$列方程组求解，但由于涉及到$n$次方运算，我们会得到一个范氏矩阵，而它往往是病态的，故考虑使用其它的方法求解.
事实上凡是要解方程求Lagrange多项式的，都会面临病态的问题，因此我们从方程形式上来考察它：
Lagrange多项式可以认为是一组线性无关的基$1, x, x^2, x^3, …, x^n$线性组合而成的，因此我们可以构造另外一组与其同构的基来解决这一问题. 我们构造这样的一组基底：
\[{l_i(x), i\in [0,n]}, s.t:l_i(x) = \frac{\prod_{j=0, j\not ={i}}^n(x-x_j)}{\prod_{j=0, j\not ={i}}^n(x_i-x_j)}\]这组基满足一个特殊的性质：对于点集中的若干离散$x$，当且仅当$x = x_i$时$l_i(x)$为$1$，否则为$0$，故我们可以构造插值函数如下：
\[f(x) = \sum_{i = 0}^n y_il_i(x)\tag{4.8}\]容易得知上述的函数满足Lagrange插值的基本假设（式4.6）
但是基于4.8的Lagrange插值法具有一明显的弊端：即一旦点集发生改变，各个基均需要进行修改，故我们希望找到另外一组基，使得即使点集发生了增减，我们也只需简单地增加或删除一项即可：
Newton插值法
课本中对Newton插值的提法太过啰嗦，故给出一种额外的提法：
我们以一种逐步构建的方式来提出Newton插值.
1 最初，我们的点集中只有$(x_0, y_0)$，因此我们给出插值式$f_0(x) = y_0$来拟合之.
2 随后我们加入了点$(x_1, y_1)$，此时插值公式需要额外满足$f(x_1)=y_1$，根据Newton插值的思想，我们通过为$f_0(x)$附加一项来达成此目的：（$s.t.$在数学符号中指”subject to”，”使得”）
\[f_1(x) = y_0 + R_1(x),\ s.t. f_1(x_0) = y_0, f_1(x_1) = y_1\]故$R_1(x)$一定满足：$R_1(x_0) = 0, R_1(x_1)=y_1-y_0$
故可以提出满足上述条件的$R_1(x)$较为简单形式：$R_1(x)=b_1(x-x_0)$，代入解出
\[b_1=\frac{y_1-y_0}{x_1 - x_0}\]此时插值公式变为
\[f_1(x) = y_0 + \frac{y_1-y_0}{x_1 - x_0}(x-x_0)\]3 继续加入点$(x_2, y_2)$，此时插值公式需要额外满足$f(x_2)=y_2$，故仍使用上述的方式：
\[f_2(x) = y_0 + \frac{y_1-y_0}{x_1 - x_0}(x-x_0) + R_2(x)\]同理，设置$R_2(x)=b(x-x_0)(x-x_1)$，解出
\[b_2 = \frac{\frac{y_2-y_0}{x_2-x_0}-\frac{y_1-y_0}{x_1-x_0}}{x_2-x_1}\]这个形式不够简洁，经过简单的数学变换之后（这种变换蕴藏着“均差(差商)的轮换对称性”，后述），可以化为
\[b_2 = \frac{\frac{y_2-y_1}{x_2-x_1}-\frac{y_1-y_0}{x_1-x_0}}{x_2-x_0}\]我们发现它有着某种对称性：我们定义$f[a, b] = \frac{f(a)-f(b)}{a-b}$，则$b_2$可以转写为：
\[b_2=\frac{f[x_2, x_1]-f_[x_1,x_0]}{x_2-x_0}\]这样的形式仍符合着某种更”高级”的对称性，故我们将$b_2$记为$f[x_0, x_1, x_2]$，表示先求”较为低级的”两个$f[x_2, x_1]$和$f[x_1, x_0]$，再将其作差，与外层的两个值之差$x_2-x_0$相除. 这样，我们递归地提出了均差(差商)的概念：
- $f[a, b] = \frac{f(a) - f(b)}{a-b}$
- $f[a, a_1, …, b_1, b] = \frac{f[a,…,b_1]-f[a_1,…,b]}{a-b}$，”$…$”可以为空且$a_1, b_1$可以相同(即本式支持三个及以上项目)
容易证明，均差具有轮换对称性，即只要参数的组合一致，参数的排列与均差的值无关.
故Newton插值法可以直观地表述为：
\[f^*(x) = f(x_0) + (x-x_0)f[x_0,x_1] + (x-x_0)(x-x_1)f[x_0,x_1,x_2] + ...\]分析地描述为：
\[f^*(x) = \sum_{i=0}^n l_i(x)\tag{4.9-1}\] \[l_i(x) = f[x_0\ to.\ x_i]\prod_{j=0}^{i-1}(x-x_j)\tag{4.9-2}\]请务必注意上述式子中连乘符号上的$(i-1)$.
Newton插值的截断误差
在非节点的任意点$(x, f(x))$上，Newton插值的截断误差可以表述为：
\[R(x) = f^*(x) - f(x)\]而$f(x)$不可直接求知，因此采用一巧妙的解法：以点$(x,f(x))$为第$(n+1)$个节点构造一新的Newton插值式$f_{(n+1)}(x)$，而此式可以精确反映$f(x)$的值. 则易知其差值：
\[R(x) = f[x_0\ to.\ x_n,x]\prod_{j=0}^{n}(x-x_0)\]即为Newton插值的截断误差. 而对于一个点系的插值式是唯一的，因此结合Lagrange插值法可以导出：
\[f[x_0\ to.\ x_n,x]\omega_{n+1}(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x)\]其中$\xi$在点集${(x_i,y_i)}$连成的开区域内部，即
\[f[x_0\ to.\ x_n,x] = \frac{f^{(n+1)}(\xi)}{(n+1)!}\]等距均差和等距节点的Newton插值
如果上述过程中的节点都是等距的，即$\forall i\in[0, n-1], x_i - x_{i+1}\equiv h$，则均差可以改写为如下的形式(证明是简单的)：
\[f[x_a, x_{a+1}, ..., x_{a+n}] = \frac{\Delta^n_{a+n, a} f}{n!h^n}\]其中$\Delta^n f$被称为$n$阶差分，根据差分与均差在内容上的一致性，我们可以递归地描述差分：
- $\Delta^1_{a+1, a}f =f(a+1) - f(a)$
- $\Delta^n_{a+n, a}f=\Delta^{n-1}{a+n-1, a}f-\Delta^{n-1}{a+n, a+1}f$
至于所谓向前差分或向后差分，我们可以通过互换下角标顺序来描述(反正阶数又不会变)，又何必创造两个不同的符号呢？数学本来是简洁明快又可爱的，为什么要弄出这么多分裂的概念来让她变得啰嗦而丑陋呢？如是，我们即可重新以差分来改写节点等距情况下的Newton插值式：
\[f^*(x) = \sum_{i=0}^n l_i(x)\tag{4.10-1}\] \[l_i(x) = \frac{\Delta^i_{i,0}f}{i!h^i} \prod_{j=0}^{i-1}(x-x_0)\tag{4.10-2}\]直观地表述为
\[f^*(x) = f(x_0)+\frac{\Delta_{1,0}^1f}{1!h^1}(x-x_0)+\frac{\Delta_{2,0}^2f}{2!h^2}(x-x_0)(x-x_1)+\frac{\Delta_{3,0}^3f}{3!h^3}(x-x_0)(x-x_1)(x-x_2)\] \[...\]而在某些情况下，$x$会使用类似于等差数列的形式表示：$x=x_0+ht$，这时我们可以使用$t$作自变量，将上述的式子再度改写：
\[f^*(x_0+ht) = \sum_{i=0}^n l_i(t)\tag{4.11-1}\] \[l_i(t) = \frac{\Delta^i_{i,0}f}{i!h^i} \prod_{j=0}^{i-1}(x_0+th-(x_0+jh))=\Delta^i_{i,0}f\frac{\prod_{j=0}^{i-1}(t-j)}{i!}\tag{4.10-2}\]直观地表述为
\[f^*(x_0+ht) = f(x_0) + t\Delta_{1,0}^1f + \frac{t(t-1)}{2!}\Delta_{2,0}^2f+\frac{t(t-1)(t-2)}{3!}\Delta_{3,0}^3f+...\]这一式子被称为向前插值公式. 代入求值时应注意自变量是$t$，若输入的是$x$则需相应的线性转化.
同理，若将$x$表示为$x=x_n+ht$，也可以得出类似的结果：
\[f^*(x_n+ht) = \sum_{i=0}^n l_i(t)\tag{4.11-1}\] \[l_i(t) = \frac{\Delta^i_{0,i}f}{i!h^i} \prod_{j=0}^{i-1}(x_0+th-(x_0+jh))=\Delta^i_{0,i}f\frac{\prod_{j=0}^{i-1}(t-j)}{i!}\tag{4.10-2}\]直观地表述为
\[f^*(x_0+ht) = f(x_0) + t\Delta_{0,1}^1f + \frac{t(t-1)}{2!}\Delta_{0,2}^2f+\frac{t(t-1)(t-2)}{3!}\Delta_{0,3}^3f+...\]向前和向后插值都是基于式子(4.10)的，只不过对于$x$的等差数列表示不同.