非对称特征值问题的计算方法

目录

• 幂法
  • 多项式的根
• 反幂法
  • 半次迭代法
  • 复特征值迭代
• QR 方法
  • 基本迭代与收敛性
  • 实 Schur 标准型
  • 上 Hessenberg 化
  • 带原点位移的 QR 迭代
  • 双重步位移的 QR 迭代
  • 隐式 QR 算法

 

首先简单介绍几个预备知识：

 

多项式 \(p_A(\lambda) = \det(\lambda I-A)\) 称为 \(A\) 的特征多项式，记 \(\lambda(A)\) 为 \(A\) 的特征值全体，通常称之为 \(A\) 的谱集；假定 \(p_A(\lambda)\) 有分解

\[p_A(\lambda) = (\lambda-\lambda_1)^{n_1}\cdots(\lambda-\lambda_p)^{n_p} \]

其中 \(n_1+n_2+\cdots+n_p=n,\ \lambda_i\neq \lambda_j\ (i\neq j)\) ，则称 \(n_i\) 为 \(\lambda_i\) 的代数重数，称数

\[m_i = n -\mathrm{rank}(\lambda_i I-A) \]

为 \(\lambda_i\) 的几何重数；若 \(n_i=1\) ，则称 \(\lambda_i\) 是 \(A\) 的一个单特征值，否则称为重特征值；若 \(n_i=m_i\) ，则称其是半单特征值。

显然单特征值必是半单特征值，事实上，此时 \(\lambda_i\) 的特征子空间是 \(1\) 维，故 \(\mathrm{rank}(\lambda_i I-A)=n-1,\ m_i = 1\) ； 如果所有特征值都是半单的，则称 \(A\) 是非亏损的； \(A\) 非亏损当且仅当 \(A\) 可相似对角化.

 

定理 6.1.1 （Jordan 分解定理） 设 \(A\in\mathbb{C}^{n\times n}\) 有 \(r\) 个互不相同的特征值 \(\lambda_1,\cdots, \lambda_r\) ，其重数分别为 \(n(\lambda_1),\cdots, n(\lambda_r)\) ，则存在非奇异阵 \(P\in\mathbb{C}^{n\times n}\) ，使得

\[J = P^{-1}AP = \left( \begin{aligned} J(\lambda_1)\\ & J(\lambda_2)\\ & & \ddots\\ & & & J(\lambda_r) \end{aligned} \right) \]

其中 \(J(\lambda_i)\) 是若当块，矩阵 \(J\) 称为 \(A\) 的若当标准型.

 

定理 6.1.2 （Schur 分解定理） 设 \(A\in\mathbb{C}^{n\times n}\) ，则存在酉矩阵（列向量是 \(\mathbb{C}^n\) 中的标准正交基） \(U\in\mathbb{C}^{n\times n}\) ，使得

\[U^*AU = T \]

其中 \(T\) 是上三角阵.

 

幂法

我们先假定 \(A\in\mathbb{C}^{n\times n}\) 可对角化，则有

\[A = X\Lambda X^{-1} \]

其中 \(\Lambda\) 是对角阵， \(X\) 非奇异，再假定

\[|\lambda_1| > |\lambda_2|\ge \cdots\ge |\lambda_n| \]

那么任取向量 \(u_0\in\mathbb{C}^{n\times n}\) ，由于 \(X\) 的列向量构成一组基，从而

\[u_0 = \alpha_1x_1+\alpha_2x_2+\cdots+\alpha_nx_n \]

那么有

\[A^ku_0 = \sum_{j=1}^n\alpha_j A^kx_j = \sum_{j=1}^n\alpha_j\lambda_j^kx_j = \lambda_1^k\left[\alpha_1x_1+\sum_{j=2}^n\alpha_j\left(\dfrac{\lambda_j}{\lambda_1}\right)^kx_j\right] \]

这样得到

\[\lim_{k\to\infty}\dfrac{A^ku_0}{\lambda_1^k} = \alpha_1x_1 \]

这意味着当 \(\alpha_1\neq 0\) 且 \(k\) 足够大

\[u_k = \dfrac{A^ku_0}{\lambda_1^k} \]

是一个很好的近似特征向量.

 

观察上述表达式，我们不能事先知道 \(A\) 的最大特征值，并且 \(A\) 的幂计算量太大。注意到我们只需要 \(u_k\) 的方向，因此不必每次用 \(\lambda_1\) 进行约化，但是为了防止溢出，约化是必要的；因此我们有迭代格式

\[\begin{aligned} y_k &= Au_{k-1}\\ \mu_k &= \zeta_j^{(k)},\quad \|y_k\|_{\infty} = \left|\zeta_j^{(k)}\right|\\ u_k &= y_k/\mu_k \end{aligned} \]

其中 \(u_0\) 是初始迭代向量，通常要求 \(\|u_0\|_{\infty} = 1\) .

 

定理 6.2.1 设 \(A\in\mathbb{C}^{n\times n}\) 有 \(p\) 个互不相同的特征值满足 \(|\lambda_1|>|\lambda_2|\ge\cdots\ge|\lambda_p|\) ，且模最大特征值 \(\lambda_1\) 是半单的，如果初始向量 \(u_0\) 在 \(\lambda_1\) 的特征子空间上的投影不为零，则幂法产生的向量序列 \(\{u_k\}\) 收敛到 \(\lambda_1\) 的一个特征向量 \(x_1\) ，且数值序列 \(\{\mu_k\}\) 收敛到 \(\lambda_1\) .

证明 假设有分解

\[A = X\mathrm{diag}(J_1,\cdots,J_p)X^{-1} \]

其中 \(J_i\) 是若当块，且 \(n_1+n_2+\cdots+n_p = n\) ，则 \(\lambda_1\) 半单意味着 \(J_1 = \lambda_1I_{n1}\) ，令 \(y = X^{-1}u_0\) ，并作分块

\[y = \left(y_1^T,y_2^T,\cdots,y_p^T\right),\quad X = (X_1,X_2,\cdots,X_p) \]

则我们有

\[\begin{aligned} A^ku_0 &= X\mathrm{diag}(J_1,\cdots,J_p)X^{-1}u_0\\ &= X_1J_1^ky_1 + X_2J_2^ky_2+\cdots +X_pJ_p^ky_p\\ &= \lambda_1^kX_1y_1 + X_2J_2^ky_2+\cdots +X_pJ_p^ky_p\\ &= \lambda_1^k\left[X_1y_1+X_2\left(\dfrac{J_2}{\lambda_1}\right)^ky_2+\cdots+X_p\left(\dfrac{J_p}{\lambda_1}\right)^ky_p\right] \end{aligned} \]

注意到右边 \(\lambda_1^{-1}J_i\) 的谱半径小于 \(1\) ，因此有

\[\lim_{k\to\infty}\dfrac{A^ku_0}{\lambda_1^k} = X_1y_1 \]

根据假设 \(u_0\) 在 \(\lambda_1\) 的特征子空间上的投影不为零，则 \(X_1y_1\neq 0\) ；

​注意到迭代格式中 \(\|u_k\|_{\infty}=1\) ，且

\[u_k = \dfrac{Au_{k-1}}{\mu_k} = \dfrac{A^ku_0}{\mu_k\mu_{k-1}\cdots\mu_1} \]

每一步 \(\mu_k\) 都是模最大分量，则 \(\zeta_k = \mu_k\mu_{k-1}\cdots\mu_1\) 是 \(A^ku_0\) 的一个模最大分量，从而 \(\xi_k/\lambda_1^k\) 是 \(A^ku_0/\lambda_1^k\) 的模最大分量，上面我们知道 \(A^ku_0/\lambda_1^k\) 的极限存在，因此

\[\zeta = \lim_{k\to\infty}\dfrac{\zeta_k}{\lambda_1^k} \]

存在，从而

\[\lim_{k\to\infty}u_k = \lim_{k\to\infty}\dfrac{A^ku_0}{\zeta_k} = \lim_{k\to\infty}\left(\dfrac{A^ku_0}{\lambda_1^k}\bigg/\dfrac{\zeta_k}{\lambda_1^k}\right) = \dfrac{X_1y_1}{\zeta} = x_1 \]

显然 \(x_1\) 是属于 \(\lambda_1\) 的特征向量，再由 \(Au_{k-1} = \mu_ku_k\) 和 \(x_1\) 有一个模最大分量为 \(1\) ，则有 \(\mu_k\to\lambda_1\) .

​幂法的收敛速度由 \(|\lambda_2|/|\lambda_1|\) 决定。

 

多项式的根

应用幂法求取多项式

\[P(z) = z^n+\alpha_1z^{n-1} + \cdots + \alpha_n \]

模最大根，我们利用 Laplace 展开式

\[\begin{aligned} P(z) &= \left| \begin{matrix} 1 & \alpha_1 & \cdots & \alpha_{n-1} & \alpha_n\\ -1 & z & \cdots & 0 & 0\\ \vdots & \vdots & & \vdots & \vdots\\ 0 & 0 & \cdots & z & 0\\ 0 & 0 & \cdots & -1 & z \end{matrix} \right| = \left| \begin{matrix} 1 & 0 & \cdots & 0 & 0\\ -1 & z+\alpha_1 & \cdots & z+\alpha_{n-1} & z+\alpha_n\\ \vdots & \vdots & & \vdots & \vdots\\ 0 & 0 & \cdots & z & 0\\ 0 & 0 & \cdots & -1 & z \end{matrix} \right| \\ &= \left| \begin{matrix} z+\alpha_1 & z+\alpha_2 & \cdots & z+\alpha_{n-1} & z+\alpha_n\\ -1 & z & \cdots & 0 & 0\\ \vdots & \vdots & & \vdots & \vdots\\ 0 & 0 & \cdots & z & 0\\ 0 & 0 & \cdots & -1 & z \end{matrix} \right| \end{aligned} \]

因此我们得到 \(P(z) = \det(zI + A)\) ，其中

\[A = \left( \begin{matrix} \alpha_1 & \alpha_2 & \cdots & \alpha_{n-1} & \alpha_n\\ -1 & 0 & \cdots & 0 & 0\\ \vdots & \vdots & & \vdots & \vdots\\ 0 & 0 & \cdots & 0 & 0\\ 0 & 0 & \cdots & -1 & 0 \end{matrix} \right) \]

从而只需对 \(A\) 应用幂法求取最大模特征值，则该特征值的相反数就是 \(P(z)\) 的根.

 

反幂法

​反幂法又称为反迭代法，也就是应用幂法于 \(A^{-1}\) 上求解 \(A\) 的模最小特征值和对应的特征向量，因此有迭代格式

\[\begin{aligned} &Ay_k= z_{k-1}\\ &\mu_k = \zeta_j^{(k)},\quad \|y_k\|_{\infty} = \left|\zeta_j^{(k)}\right|\\ &z_k = y_k/\mu_k \end{aligned} \]

由之前的讨论知道，若 \(A\) 有特征值 \(|\lambda_n|<|\lambda_{n-1}|\le\cdots\le|\lambda_1|\) ，则 \(\{z_k\}\) 收敛到 \(A\) 对应于 \(\lambda_n\) 的一个特征向量，而 \(\{\mu_k\}\) 收敛于 \(\lambda_n^{-1}\) .

​在实际应用中，主要应用反幂法求解特征向量，在求得特征值 \(\lambda_i\) 的近似值 \(\widetilde{\lambda}_i\) 后，应用反幂法于 \(A-\widetilde{\lambda}_iI\) 上，则有

\[(A-\mu I)v_k = z_{k-1},\quad z_k = v_k/\|v_k\|_{\infty} \]

只需要每次求解线性方程组即可。

​如果假设特征值有

\[0 <|\lambda_1-\mu|<|\lambda_2-\mu|\le|\lambda_3-\mu|\le\cdots\le|\lambda_n-\mu| \]

则收敛速度取决于 \(|\lambda_1-\mu|/|\lambda_2-\mu|\) ，则 \(\mu\) 越靠近 \(\lambda_1\) ，收敛越快；但是此时 \(A-\mu I\) 会离一个奇异矩阵非常近，每次要解一个非常病态的线性方程组，好在其病态性不影响收敛速度，并且常常只需迭代一次就可以得到相当好的近似特征向量。事实上，求解该线性方程组的误差主要对解在特征子空间上的长度有影响，但我们只关心其方向。

 

半次迭代法

我们可以随机地选取向量，进行规范化后作为初始向量，但更好的是采用半次迭代法：首先进行分解

\[P(A-\mu I) = LU \]

那么第一次迭代为

\[P^TLUv_1 = z_0 \]

只需令 \(z_0 = P^TLe\) ，其中 \(e\) 是分量全为 \(1\) 的向量，这样就只需要求解三角方程组

\[Uv_1 = e \]

这样一来，不需要明确给出初始向量 \(z_0\) ，而完成这一次迭代就只需求解 “半次” 线性方程组.

 

复特征值迭代

​应用反幂法于拟上三角阵 \(A\in\mathbb{R}^{n\times n}\) ，具有形式

\[A = \left( \begin{matrix} A_{11} & A_{12} & \cdots & A_{1k}\\ & A_{12} & \cdots & A_{2k}\\ & & \ddots & \vdots\\ & & & A_{kk} \end{matrix} \right) \]

其中 \(A_{ii}\) 是实数或二阶矩阵，考虑一种求取 \(A\) 全部特征向量的数值方法；

​首先，对于其中每一个 \(A_{ii}\) ，如果它是实数，直接应用反幂法即可；否则，我们假设其有特征值 \(a+bi\) ，则反幂法有格式

\[\begin{aligned} &(A-aI-ibI)y_{k+1} = u_k\\ &u_{k+1} = y_{k+1}/\mu_{k+1} \end{aligned} \]

可以将虚实部分开迭代，设 \(y_{k+1} = x_{k+1}+iz_{k+1},\ u_k=p_k+iq_k\) ，那么

\[\begin{aligned} &(A-aI)x_{k+1}+bz_{k+1} = p_k\\ &-bx_{k+1}+(A-aI)z_{k+1} = q_k \end{aligned} \]

分别消元得到

\[\begin{aligned} &\left[(A-aI)^2+b^2I\right]x_{k+1} = (A-aI)p_k-bq_k\\ &\left[b^2I+(A-aI)^2\right]z_{k+1} = bp_k +(A-aI)q_k \end{aligned} \]

只需要分别解出两式即可.

 

QR 方法

QR 方法是利用正交相似变换将给定矩阵逐步约化为上三角阵或拟上三角阵的迭代方法，其基本收敛速度是二次的.

 

基本迭代与收敛性

​给定 \(A_0 = A\in\mathbb{C}^{n\times n}\) ， QR 算法的基本迭代格式为

\[\begin{aligned} &A_{m-1} = Q_mR_m,\\ &A_m = R_mQ_m \end{aligned} \quad m = 1,2,\cdots \]

其中 \(Q_m\) 为酉矩阵， \(R_m\) 为上三角阵。为了分析方便，我们先假定 \(R_m\) 对角元非负，则

\[A_m = Q_m^*A_{m-1}Q_m \]

即矩阵序列 \(\{A_m\}\) 中每个矩阵都与 \(A\) 相似。反复应用迭代格式有

\[A_m = \widetilde{Q}_m^*A\widetilde{Q}_m,\quad \widetilde{Q}_m = Q_1Q_2\cdots Q_m \]

再代入 \(A_m = R_mQ_m\) 有

\[\widetilde{Q}_mQ_{m+1}R_{m+1} = A\widetilde{Q}_m \]

从而有

\[\widetilde{Q}_mQ_{m+1}R_{m+1}R_m\cdots R_1 = A\widetilde{Q}_mR_m\cdots R_1 \]

得到递推关系

\[\widetilde{Q}_{m+1}\widetilde{R}_{m+1} = A\widetilde{Q}_m\widetilde{R}_{m},\quad \widetilde{R}_{m} = R_{m}\cdots R_1 \]

由此即得

\[A^m = \widetilde{Q}_{m}\widetilde{R}_{m} \]

​利用这一基本关系式，我们导出 QR 算法和幂法的关系。记 \(\widetilde{R}_{m}\) 元素为 \(\gamma_{ij}\) ， \(\widetilde{Q}_{m}\) 第一列为 \(q_1^{(m)}\) ，则有

\[A^me_1 = \gamma_{11}q_1^{(m)} \]

那么 \(q_1^{(m)}\) 可看做以 \(e_1\) 为初始向量的幂法得到的向量。若 \(A\) 的模最大特征值 \(\lambda_1\) 与其它特征值分离，则 \(q_1^{(m)}\) 将收敛到对应的特征向量.

 

定理 6.4.1 设 \(A\) 的 \(n\) 个特征值满足 \(|\lambda_1|>|\lambda_2|>\cdots>|\lambda_n|>0\) ，并设 \(n\) 阶方阵 \(Y\) 的第 \(i\) 行是 \(\lambda_i\) 的特征向量，如果 \(Y\) 有 LU 分解，则上述迭代格式得到的矩阵序列 \(\{A_m\}\) 对角线以下元素趋于 \(0\) ，并且对角元 \(\alpha_{ii}^{(m)}\) 趋于 \(\lambda_i\) .

 

实 Schur 标准型

​我们希望有只涉及实数运算的 QR 迭代，但由于复共轭特征值的存在，我们不能期望迭代格式得到的矩阵序列逼近上三角阵，事实上我们有

定理 6.4.2 （实 Schur 分解） 设 \(A\in\mathbb{R}^{n\times n}\) ，则存在正交阵 \(Q\) 使得

\[Q^TAQ = \left( \begin{matrix} R_{11} & R_{12} & \cdots & R_{1m}\\ & R_{12} & \cdots & R_{2m}\\ & & \ddots & \vdots\\ & & & R_{mm} \end{matrix} \right) \]

其中 \(R_{ii}\) 是一个实数或具有复共轭特征值的 \(2\) 阶方阵.

 

​我们自然想到迭代格式

\[\begin{aligned} &A_{m-1} = Q_mR_m,\\ &A_m = R_mQ_m \end{aligned} \quad m = 1,2,\cdots \]

其中 \(Q_m\) 是正交阵， \(R_m\) 是上三角阵，并且此迭代格式收敛于实 Schur 标准型；

​然而，此迭代格式每次迭代运算量太大，收敛太慢，因此并不实用，我们需要进行优化.

 

上 Hessenberg 化

​实际计算时，总是将 \(A\) 约化为准上三角阵，再对约化后的矩阵进行迭代；

​设 \(A\in\mathbb{R}^{n\times n}\) ，我们希望得到一个非奇异矩阵 \(Q\) 使得

\[\widetilde{A} = QAQ^{-1} \]

具有某种特殊形式，我们考虑 Householder 变换 \(H_1\) 使 \(A\) 的第一列有尽可能多的零元素，然而为了保证相似变换，还需要右乘 \(H_1\) ，为了保证 \(H_1A\) 第一列的零元素不被破坏掉，选取

\[H_1 = \left( \begin{matrix} 1 & 0\\ 0 & \widetilde{H}_1 \end{matrix} \right) \]

这样一来，就有

\[H_1AH_1 = \left( \begin{matrix} \alpha_{11} & a_2^T\widetilde{H}_1\\ \widetilde{H}_1a_1 & \widetilde{H}_1A_{22}\widetilde{H}_1 \end{matrix} \right) \]

只需要选取合适的变换使得

\[\widetilde{H}_1a_1 = pe_1 \]

这样第一列就有 \(n-2\) 个零元素，以此类推得到上 Hessenberg 矩阵

\[H = \left( \begin{matrix} h_{11} & h_{12} & h_{13} & \cdots & h_{1,n-1} & h_{1n}\\ h_{21} & h_{22} & h_{23} & \cdots & h_{2,n-1} & h_{2n}\\ & h_{23} & h_{33} & \cdots & h_{3,n-1} & h_{3n}\\ & & \ddots & \ddots & \vdots & \vdots\\ & & & \ddots & \vdots & \vdots\\ & & & & h_{n,n-1} & h_{nn} \end{matrix} \right) \]

​现在令

\[Q_0 = H_1H_2\cdots H_{n-2} \]

则有

\[Q_0^TAQ_0 = H \]

该分解式就称为 \(A\) 的上 Hessenberg 分解.

 

定理 6.4.3 设 \(A\in\mathbb{R}^{n\times n}\) 有如下两个上 Hessenberg 分解

\[U^TAU = H,\quad V^TAV = G \]

若 \(U,\ V\) 第一列相同，且 \(H\) 的次对角元 \(h_{i+1,i} \neq 0\) ，则存在对角元均为 \(1\) 或 \(-1\) 的对角阵 \(D\) 使得

\[U = VD,\quad H = DGD \]

一个上 Hessenberg 矩阵 \(H\) ，如果次对角元 \(h_{i+1,i} \neq 0\) ，则称它不可约。上述定理表明，如果 \(Q^TAQ = H\) 不可约，则 \(Q,\ H\) 在相差一个正负号的意义下，由 \(Q\) 的第一列确定.

 

​我们尝试通过平面旋转实现对 \(H\) 的 QR 分解，事实上，我们只需要对对角块应用 Givens 变换即可，得到

\[P_{n-1,n}P_{n-2,n-1}\cdots P_{12}H = R \]

其中 \(P_{ij}\) 是对 \((i,j)\) 平面的旋转变换， \(R\) 是上三角阵，令 \(Q^T = P_{n-1,n}P_{n-2,n-1}\cdots P_{12}\) ，则 \(H = QR\) ，完成分解；接着进行的迭代

\[\widetilde{H} = RQ = RP_{12}^TP_{23}^T\cdots P_{n-1,n}^T \]

容易发现 \(\widetilde{H}\) 仍为上 Hessenberg 矩阵.

 

带原点位移的 QR 迭代

​由于基本 QR 迭代是线性收敛的，其收敛速度取决于特征值之间的分离程度。为了加速收敛，可考虑进行原点位移。

​设第 \(m\) 步迭代位移为 \(\mu_m\) ，则带原点位移的 QR 迭代格式为

\[\begin{aligned} &H_{m}-\mu_m I = Q_mR_m,\\ &H_{m+1} = R_mQ_m + \mu_mI \end{aligned} \quad m = 0,1,\cdots \]

其中 \(H_0\) 是给定的上 Hessenberg 矩阵.

​由于 \(H_m\) 是上 Hessenberg 阵，因此如果它收敛，那么其最后一行的两个非负元素中 \(h_{n,n-1}\) 很小， \(h_{nn}\) 接近 \(H\) 的特征值，这样我们选取 \(\mu_m = h_{nn}^{(m)}\) ，可以证明若 \(h_{n,n-1}^{(m)}=\epsilon\) ，则 \(h_{n,n-1}^{(m+1)}=O(\epsilon^2)\) ，事实上只需要考虑最后一步 Givens 变换作用于

\[\widetilde{H}_{m2} = \left( \begin{matrix} \alpha & \beta\\ \epsilon & 0 \end{matrix} \right) \]

注意右下角 \(h_{nn}^{(m)}\) 已被位移消除，则变换结果得到

\[h_{n,n-1}^{(m+1)} = -s^2\beta = -\dfrac{\beta}{\sigma^2}\epsilon^2 \]

从而加速为二次收敛.

 

双重步位移的 QR 迭代

​上面的带原点位移的 QR 迭代有严重缺陷：若 \(A\) 具有复共轭特征值，则实位移一般不能起到加速作用，于是我们考虑将两步带原点位移的 QR 迭代合为一步，以避免复数运算。

​设 \(A\in\mathbb{R}^{n\times n}\) ，考虑迭代格式

\[\begin{aligned} &H_1 = Q_0^TAQ_0,\\ &H_k-\mu_kI=Q_kR_k,\\ &H_{k+1} = R_kQ_k+\mu_kI \end{aligned} \quad k = 1,2,\cdots \]

我们假定迭代格式中的上 Hessenberg 阵都不可约，否则可以进行分块迭代；

​由于实矩阵可能有复特征值，因此对于有复共轭特征值 \(\mu_1,\ \mu_2\) 的右下角矩阵，不能期望 \(h_{nn}^{(k)}\) 收敛于特征值，这时应当选取某个特征值位移，但这样会涉及复数运算。为了避免这一问题，可以同时进行两次位移

\[\begin{aligned} &H - \mu_1I = U_1R_1,\quad H_1 = R_1U_1 + \mu_1I,\\ &H_1 - \mu_2I = U_2R_2,\quad H_2 = R_2U_2 + \mu_2I \end{aligned} \]

记 \(H = H_k\) ，则有

\[M = QR,\quad H_2 = Q^*HQ \]

其中

\[M = (H-\mu_1I)(H-\mu_2I),\quad Q = U_1U_2,\quad R = R_2R_1 \]

实际上可以验证

\[Q^*HQ = U_2^*U_1^*(\mu_1I+U_1R_1)U_1U_2 = U_2^*(\mu_1I+R_1U_1)U_2 = U_2^*(\mu_2I+U_2R_2)U_2 = H_2 \]

由上面的式子得到

\[M = H^2 -sH+tI,\quad s = \mu_1+\mu_2 = h_{mm}^{(k)} + h_{nn}^{(k)},\ t = \mu_1\mu_2 = \det G_k\in\mathbb{R} \]

因此 \(M\) 是实矩阵，这样一来我们自然想到计算两步迭代的步骤：

1. 计算 \(M = H^2-sH+tI\)
2. 计算 \(M = QR\) 分解
3. 计算 \(H_2 = Q^THQ\)

但第一步的计算量就已经很大，幸运的是我们有

定理 6.4.4 若 \(H\) 是不可约上 Hessenberg 矩阵，且 \(\mu_1,\ \mu_2\) 均非 \(H\) 的特征值，则 \(H_2\) 也是不可约上 Hessenberg 矩阵.

 

​根据之前的定理，只需计算一个正交阵 \(\widetilde{Q}\) 使得 \(\widetilde{Q}^TH\widetilde{Q} = \widetilde{H}_2\) ，只要 \(\widetilde{Q}\) 和 \(Q\) 第一列相同，那么 \(\widetilde{H}_2\) 和 \(H_2\) 在相差正负号的意义上是相同的。我们通过另外的途径实现 \(H\) 到 \(H_2\) 的变换；

​首先，由于 \(M = QR\) ， \(R\) 为上三角阵，因此 \(Q,\ M\) 第一列共线，实际上由于 \(Q\) 每个列向量都是标准化的，因此 \(Q\) 第一列就是 \(M\) 第一列单位化得到的。计算 \(M = H^2 -sH + tI\) 得到

\[Me_1 = (\xi_1,\xi_2,\xi_3,0,\cdots,0)^T \]

其中

\[\begin{aligned} &\xi_1 = \left(h_{11}^{(k)}\right)^2 + h_{12}^{(k)}h_{21}^{(k)} - sh_{11}^{(k)} + t\\ &\xi_2 = h_{21}^{(k)}\left(h_{11}^{(k)}+h_{22}^{(k)}-s\right)\\ &\xi_2 = h_{21}^{(k)}h_{32}^{(k)} \end{aligned} \]

​其次，考虑作用于 \(Me_1\) 的 Householder 变换 \(P_0Me_1 = \alpha e_1\) ，同样有 \(P_0\) 第一列与 \(Me_1\) 共线，我们直接把 \(P_0\) 的第一列作为 \(Q\) 的第一列，即 \(P_0e_1 = Qe_1\) ，而 \(P_0\) 可由上面的 \(Me_1\) 得到。

​现在， 我们寻找一个正交阵 \(\widetilde{Q}\) ，使其第一列为 \(e_1\) ，那么

\[P_0\widetilde{Q}e_1 = P_0e_1 = Qe_1 \]

这意味着 \(P_0\widetilde{Q}\) 和 \(Q\) 第一列相同，那么有

\[H_2 = Q^THQ,\quad \widetilde{H}_2 = (P_0\widetilde{Q})^TH(P_0\widetilde{Q}) = \widetilde{Q}^T(P_0^THP_0)\widetilde{Q} = \widetilde{Q}^T(P_0HP_0)\widetilde{Q} \]

根据之前的定理， \(H_2\) 和 \(\widetilde{H}_2\) 在差正负号的意义下相同，我们就直接取 \(H_2 = \widetilde{H}_2\) 即可；

​最后，考虑求取 \(\widetilde{Q}\) ，由 \(P_0\) 的形式，它只改变 \(H\) 的前三行和前三列，即

\[P_0HP_0 = \left( \begin{matrix} \times & \times & \times & \times & \cdots & \times & \times\\ \times & \times & \times & \times & \cdots & \times & \times\\ + & \times & \times & \times & \cdots & \times & \times\\ + & + & \times & \times & \cdots & \times & \times\\ & & & \times & \cdots & \times & \times\\ & & & & \ddots & \vdots & \vdots\\ & & & & & \times & \times \end{matrix} \right) \]

接着对第一列，进行 Householder 变换 \(P_1\) ，为了不破坏形式，我们只变换第二行到第四行，即 \(P_1 = \mathrm{diag}(1, \widetilde{P}_1, I_{n-4})\) ，则有

\[P_1P_0HP_0P_1 = \left( \begin{matrix} \times & \times & \times & \times & \cdots & \times & \times\\ \times & \times & \times & \times & \cdots & \times & \times\\ & \times & \times & \times & \cdots & \times & \times\\ & + & \times & \times & \cdots & \times & \times\\ & + & + & \times & \cdots & \times & \times\\ & & & & \ddots & \vdots & \vdots\\ & & & & & \times & \times \end{matrix} \right) \]

以此类推，每一步 \(P_k\) 变换三行三列，直到最后一步为

\[P_{n-2} = \mathrm{diag}(I_{n-2},\widetilde{P}_{n-2}) \]

这里 \(P_{n-2}\) 是二阶 Householder 变换。于是

\[\widetilde{Q} = P_1P_2\cdots P_{n-2} \]

观察得到 \(\widetilde{Q}e_1 = e_1\) ，因此得到满足条件的正交阵.

 

隐式 QR 算法

​我们给出上述方法收敛的判定准则：当

\[|h_{i+1,i}|\le (|h_{ii}|+|h_{i+1,i+1}|)\epsilon_u \]

时，就将 \(h_{i+1,i}\) 看做零；接着给出如下隐式 QR 算法，它计算 \(n\) 阶实矩阵 \(A\) 的实 Schur 分解 \(Q^TAQ = T\) ，其中 \(T\) 是拟上三角阵，即对角块为实数或二阶方阵的上三角阵.

 

算法 6.4.3 （隐式 QR）

• 计算 \(A\) 的上 Hessenberg 分解： \(H = U_0^TAU_0,\ Q=U_0\)

• 收敛性判定

  • 将满足 \(|h_{i+1,i}|\le (|h_{ii}|+|h_{i+1,i+1}|)\epsilon_u\) 的 \(h_{i+1,i}\) 归零
  • 确定最大的非负整数 \(m\) 和最小的非负整数 \(l\) 使得
  \[H = \left( \begin{matrix} H_{11} & H_{12} & H_{13}\\ 0 & H_{22} & H_{23}\\ 0 & 0 & H_{33} \end{matrix} \right) \begin{aligned} &\tiny{l}\\ &\tiny{n-l-m}\\ &\tiny{m} \end{aligned} \]

  其中 \(H_{33}\) 是约化完成的拟上三角阵，而 \(H_{22}\) 是不可约上 Hessenberg 矩阵

  • 如果 \(m=n\) 退出，否则继续

• 进行 QR 迭代，然后返回收敛性判定