您可以自由地共享和演绎本文稿, 只需遵守开源协议[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$的等差数列表示不同.