Loading

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

 

基本性质

定理 7.1.1 (谱分解定理)\(A\in\mathbb{R}^{n\times n}\) 是对称的,则存在正交阵 \(Q\) 使得

\[Q^TAQ = \Lambda = \mathrm{diag}(\lambda_1,\cdots,\lambda_n) \]

定理 7.1.2 (极小极大定理)\(A\in\mathbb{R}^{n\times n}\) 是对称阵,并有特征值 \(\lambda_1\ge\cdots\ge\lambda_n\) ,则有

\[\lambda_i = \max_{\mathcal{X}\in\mathcal{G}_i^n}\min_{0\neq u\in\mathcal{X}}\dfrac{u^TAu}{u^Tu} = \min_{\mathcal{X}\in\mathcal{G}_{n-i+1}^n}\max_{0\neq u\in\mathcal{X}}\dfrac{u^TAu}{u^Tu} \]

其中 \(\mathcal{G}_k^n\) 表示 \(\mathbb{R}^n\) 中所有 \(k\) 维子空间全体.

 

定理 7.1.3 (Weyl 定理)\(n\) 阶对称矩阵 \(A\)\(B\) 的特征值分别为

\[\lambda_1\ge\cdots\ge\lambda_n,\quad \mu_1\ge\cdots\ge\mu_n \]

则有

\[|\lambda_i-\mu_i|\le \|A-B\|_2,\quad i=1,\cdots,n \]

这一定理表明对称矩阵的特征值总是非常良态的.

 

定理 7.1.4\(A\)\(A+E\) 是两个 \(n\) 阶实对称矩阵,假定 \(q_1\)\(A\) 的一个单位特征向量, \(Q = [q_1,Q_2]\)\(n\) 阶正交矩阵, \(Q^TAQ\)\(Q^TEQ\) 分块如下

\[Q^TAQ = \left( \begin{matrix} \lambda & 0\\ 0 & D_2 \end{matrix} \right),\quad Q^TEQ = \left( \begin{matrix} \epsilon & e^T\\ e & E_{22} \end{matrix} \right) \]

若有

\[d = \min_{\mu\in\lambda(D_2)}|\lambda-\mu| > 0,\quad \|E\|_2\le \dfrac{1}{4}d \]

则存在 \(A+E\) 的一个单位特征向量 \(\widetilde{q}_1\) 使得

\[\sin\theta = \sqrt{1-\left|q_1^T\widetilde{q}_1\right|^2} \le \dfrac{4}{d}\|e\|_2 \le \dfrac{4}{d}\|E\|_2 \]

其中 \(\theta\)\(q_1,\ \widetilde{q}_1\) 之间所夹的锐角;

​此定理表明,特征向量的敏感性依赖于对应的特征值与其它特征值之间的分离程度.

 

定理 7.15 (SVD 分解定理)\(A\in\mathbb{R}^{m\times n}\) ,则存在正交阵 \(U\in\mathbb{R}^{m\times m}\)\(V\in\mathbb{R}^{n\times n}\) ,使得

\[U^TAV = \left( \begin{matrix} \Sigma_r & 0\\ 0 & 0 \end{matrix} \right),\quad \Sigma_r = \mathrm{diag}(\sigma_1,\cdots,\sigma_r),\ \sigma_1\ge\sigma_2\ge\cdots\ge\sigma_r>0 \]

\(\sigma_i\)\(A\) 的奇异值, \(U,\ V\) 的列向量分别称为 \(A\) 的左右奇异向量.

证明 我们考虑实对称阵 \(A^TA\) ,则有实特征值 \(\lambda_i,\ i=1,\cdots,n\) ,以及对应的标准正交特征向量 \(v_1,\cdots,v_n\) ,这是由实对称阵合同于对角元全为特征值的对角阵得到的。这样可取 \(Av_1,\cdots,Av_n\in\mathbb{R}^m\) ,我们有

\[v_i^Tv_j = \left\{ \begin{aligned} 1\quad i=j\\ 0\quad i\neq j \end{aligned} \right.\quad\quad (Av_i)^T(Av_j) = v_iA^TAv_j = \lambda_jv_i^Tv_j = \left\{ \begin{aligned} \lambda_j\quad i=j\\ 0\quad i\neq j \end{aligned} \right. \]

也就是说,只要 \(\lambda_j\neq 0\) ,就有 \(Av_i,\ Av_j\) 是正交向量;

​假设特征值按照模从大到小排列,仍记为 \(\lambda_1,\cdots,\lambda_n\) ,令 \(\mathrm{rank}(A^TA) = \mathrm{rank}(A) = r\) ,那么就有

\[Av_1,\ Av_2,\ \cdots,\ Av_r\in\mathbb{R}^m,\quad Av_{r+1} = \cdots = Av_n = 0 \]

\(r\) 项是 \(\mathbb{R}^m\) 中的正交向量,剩余项对应于 \(\lambda_i=0\) 的特征值,因此为零;

​注意到 \(0\le\|Av_i\|_2^2 = (Av_i)^T(Av_i) = \lambda_i\) ,重新进行标准化得到

\[u_1,\ u_2,\ \cdots,\ u_r,\quad u_i = \dfrac{Av_i}{\sqrt{\lambda_i}},\ i=1,\cdots,r \]

然后扩张为 \(\mathbb{R}^m\) 中的基 \(u_1,\cdots,u_m\) ,从而我们有

\[A(v_1,\cdots,v_n) = (\sqrt{\lambda_1}u_1,\cdots,\sqrt{\lambda_r}u_r,0,\cdots,0) = (u_1,\cdots,u_m)\cdot\mathrm{diag}(\sqrt{\lambda_1},\cdots,\sqrt{\lambda_r},0,\cdots,0) \]

\(V = (v_1,\cdots,v_n),\ U = (u_1,\cdots,u_m)\) ,则有分解

\[U^TAV = \mathrm{diag}(\sqrt{\lambda_1},\cdots,\sqrt{\lambda_r},0,\cdots,0) \]

此定理说明 \(A\) 作为线性运算满足

\[\mathcal{A}:\mathbb{R}^n\to\mathbb{R}^m,\quad (v_1,\cdots,v_n) \mapsto (\sqrt{\lambda_1}u_1,\cdots,\sqrt{\lambda_r}u_r,0,\cdots,0) \]

它将一组正交基映到另一组部分正交基.

 

推论 7.1.1\(A,B\in\mathbb{R}^{m\times n}\) ,假定有奇异值

\[\sigma_1\ge\cdots\ge\sigma_n,\quad \tau_1\ge\cdots\ge\tau_n \]

则有

\[|\sigma_i-\tau_i| \le \|A-B\|_2,\quad i=1,\cdots,n \]

这表明奇异值也是良态的.

 

对称 QR 方法

我们将 QR 方法应用于对称矩阵,并充分利用其对称性。

 

三对角化

​若 \(A\)\(n\) 阶实对称矩阵,并假定有上 Hessenberg 分解

\[Q^TAQ = T \]

\(T\) 必为对称三对角阵。事实上,设 \(Q=(\xi_1,\cdots,\xi_n)\) ,则有

\[(Q^TAQ)_{ij} = \xi_i^TA\xi_j = \left(\xi_i^TA\xi_j\right)^T = \xi_j^TA\xi_i = (Q^TAQ)_{ji} \]

因此 \(T\) 是对称阵,由上 Hessenberg 分解的形式即证;

​考虑简化计算每一步 \(H_kA_{k-1}H_k\) ,设

\[H_k = I -\beta vv^T,\quad v\in\mathbb{R}^{n-k} \]

这里 \(H_k\) 是作用于 \(A\)\(n-k\) 行列的元素。由对称性得到

\[H_kA_{k-1}H_k = A_{k-1} - v\omega^T - \omega v^T \]

其中

\[\omega = u - \dfrac{1}{2}\beta\left(v^Tu\right)v,\quad u = \beta A_{k-1}v \]

从而简化了计算过程.

 

隐式对称 QR 迭代

​由于 \(A\) 的特征值均为实数,因此没有必要再使用双重步位移迭代,有迭代格式

\[\begin{aligned} &T_k -\mu_kI = Q_kR_k\\ &T_{k+1} = R_kQ_k + \mu_kI \end{aligned} \]

在每一步位移时,最简单的选取 \(\mu_k\) 的方法是 \(T_k(n,n)\) ,但更好的是取 \(T_k\) 右下角的二阶矩阵

\[\left( \begin{matrix} \alpha_{n-1} & \beta_{n-1}\\ \beta_{n-1} & \alpha_n \end{matrix} \right) \]

的两个特征值中靠近 \(\alpha_n\) 的一个,即为

\[\mu_k = \alpha_n + \delta - \mathrm{sgn}(\delta)\sqrt{\delta^2+\beta_{n-1}^2},\quad \delta = (\alpha_{n-1}-\alpha_n)/2 \]

这就是 Wilkinson 位移.

​我们考虑实现一次 QR 迭代,虽然可以直接用 Givens 变换进行分解,但更好的是用隐含的方式实现变换。观察一步迭代

\[T - \mu I = QR,\quad \widetilde{T} = RQ + \mu I \]

容易验证 \(\widetilde{T} = Q^TTQ\) ,事实上

\[Q^TTQ = Q^T(\mu I +QR)Q = \mu I + RQ = \widetilde{T} \]

类似于非对称矩阵的隐式迭代,注意到在 QR 分解过程中, \(G_1^Te_1 = Qe_1\) ,只需找到一个第一列为 \(e_1\) 的正交阵 \(\widetilde{Q}\) ,则

\[G_1^T\widetilde{Q}e_1 = G_1^Te_1 = Qe_1,\quad \widetilde{T} = \left(G_1^T\widetilde{Q}\right)^TT\left(G_1^T\widetilde{Q}\right) = \widetilde{Q}^TG_1TG_1^T\widetilde{Q} \]

这意味着 \(G_1^T\widetilde{Q}\)\(Q\) 变换在差正负号的意义下相同,而

\[G_1TG_1^T = \left( \begin{matrix} \times & \times & + & 0 & \cdots \\ \times & \times & \times & 0 & \cdots \\ + & \times & \times & \times & \cdots \\ 0 & 0 & \times & \times & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{matrix} \right) \]

之后每一步 Givens 变换 \(G_k\) 都有 \(G_k^Te_1 = e_1\) ,因此最终

\[\widetilde{Q} = G_2^T\cdots G_{n-1}^T \]

就得到满足的正交阵.

 

Jacobi 方法

​Jacobi 方法是计算实对称矩阵全部特征值和特征向量最古老的方法之一,它具有编程简单、并行效率高的特点,近年来又重新受到人们的重视。

 

经典 Jacobi 方法

​设 \(A=[\alpha_{ij}]\)\(n\times n\) 实对称矩阵, Jacobi 方法的目标是将非对角范数

\[E(A) = \left(\|A\|_F^2-\sum_{i=1}^n\alpha_{ii}^2\right)^{\frac{1}{2}} = \left(\sum_{i=1}^n\sum_{j=1;j\neq i}^n\alpha_{ij}^2 \right)^{\frac{1}{2}} \]

逐步约化为零。我们通过平面旋转来进行约化 \((p,q)\) 平面

\[\left( \begin{matrix} \beta_{pp} & \beta_{pq}\\ \beta_{qp} & \beta_{qq} \end{matrix} \right) = \left( \begin{matrix} c & s\\ -s & c \end{matrix} \right)^T \left( \begin{matrix} \alpha_{pp} & \alpha_{pq}\\ \alpha_{qp} & \alpha_{qq} \end{matrix} \right) \left( \begin{matrix} c & s\\ -s & c \end{matrix} \right) \]

其中 \(\beta_{pq} = \beta_{qp} = 0\) ,则有

\[0 = \beta_{pq} = \beta_{qp} = \left(c^2-s^2\right)\alpha_{pq} + cs(\alpha_{pp}-\alpha_{qq}) \]

如果 \(\alpha_{pq} = 0\) ,只需取 \(c=1,\ s=0\) 即可;否则考虑

\[\dfrac{s^2}{c^2}\alpha_{pq} + \dfrac{s}{c}(\alpha_{qq}-\alpha_{pp}) - \alpha_{pq} = 0 \]

则可以进行变量代换

\[t = \dfrac{s}{c} = \tan\theta,\quad \tau = \dfrac{\alpha_{qq}-\alpha_{pp}}{2\alpha_{pq}} \]

则得到二次方程

\[t^2+2\tau t - 1 = 0 \]

这样就有两种可选择的根,我们选择绝对值较小的根

\[t = \dfrac{\mathrm{sgn}(\tau)}{|\tau|+\sqrt{1+\tau^2}} \]

这样选择保证了旋转角 \(|\theta|\le \pi/4\) ,这对 Jacobi 方法的收敛性至关重要;接着就可以确定

\[c = \dfrac{1}{\sqrt{1+t^2}},\quad s = tc \]

​考虑旋转平面的选取,由于 Frobenius 范数在正交变换下不变,因此 \(\|B\|_F=\|A\|_F\) ,应用于 \((p,q)\) 平面的二阶矩阵有

\[\alpha_{pp}^2+\alpha_{qq}^2+2\alpha_{pq}^2 = \beta_{pp}^2+\beta_{qq}^2+2\beta_{pq}^2 = \beta_{pp}^2+\beta_{qq}^2 \]

这样就有

\[E(B)^2 = \|B\|_F^2 - \sum_{i=1}^n\beta_{ii}^2 = \|A\|_F^2 - \sum_{i=1}^n\alpha_{ii}^2 + \alpha_{pp}^2 + \alpha_{qq}^2 - \beta_{pp}^2 - \beta_{qq}^2 = E(A)^2 - 2\alpha_{pq}^2 \]

我们的目的就是让 \(E(B)\) 尽可能小,因此应当选择模最大的元素 \(\alpha_{pq}\) ,然后选取此 \((p,q)\) 平面;

​最终我们得到迭代格式

\[A_k = J_k^TA_{k-1}J_k,\quad k=1,2,\cdots \]

其中 \(A_0 = A\)\(J_k\) 是由上述方式确定的平面旋转变换.

 

定理 7.3.1 存在 \(A\) 的特征值的一个排列 \(\lambda_1,\lambda_2,\cdots,\lambda_n\) ,使得

\[\lim_{k\to\infty}A_k = \mathrm{diag}(\lambda_1,\lambda_2,\cdots,\lambda_n) \]

证明 我们已经知道

\[E(A_{k})^2 = E(A_{k-1})^2-2\left(\alpha_{pq}^{(k-1)}\right)^2 \]

这里 \(\alpha_{pq}^{(k-1)}\)\(A_{k-1}\) 的非对角元中绝对值最大者,又注意到

\[E(A_{k-1})^2 = \|A_{k-1}\|_F^2 - \sum_{i=1}^n\alpha_{ii}^2\le n^2\alpha_{pq}^2- \sum_{i=1}^n\alpha_{ii}^2 \le n(n-1)\left(\alpha_{pq}^{(k-1)}\right)^2 \]

将其代入就得到

\[E(A_k)^2 \le \left(1-\dfrac{1}{N}\right)E(A_{k-1})^2,\quad N = \dfrac{1}{2}n(n-1) \]

这样就容易验证

\[\lim_{k\to\infty}E(A_k)^2 \le \lim_{k\to\infty}\left(1-\dfrac{1}{N}\right)^kE(A)^2 = 0 \]

即证;关于特征值的排列证明省略.

 

过关 Jacobi 方法

​每次旋转平面都只需要 \(O(n)\) 的操作,但是选择旋转平面却需要 \(O(n^2)\) 的比较操作,这是得不偿失的。因此可以按照一定顺序依次消去每个对角元,而不进行选择,这就得到循环 Jacobi 方法;

​实际计算中运用最多的是过关 Jacobi 方法:首先确定一个关值,对绝对值超过关值的非对角元进行旋转操作,直到所有非对角元小于关值,此时减小关值,重复操作直到达到足够小的关值。常用关值为

\[\delta_0 = E(A),\quad \delta_k = \dfrac{\delta_{k-1}}{\sigma},\quad k=1,2,\cdots \]

其中 \(\sigma\ge n\) 是一个固定的常数.

 

​Jacobi 方法的优势在于计算特征向量非常方便,对于 \(k\) 步变换

\[A_k = J_k^TJ_{k-1}^T\cdots J_1^TAJ_1J_2\cdots J_k \]

我们有

\[AQ_k = Q_kA_k,\quad Q_k = J_1J_2\cdots J_k \]

由于 \(A_k\) 对角元已经非常小,就可以把它们看做是近似特征值,从而 \(Q_k\) 的列向量可以看做对应的近似特征向量.

 

二分法

​我们考虑求实对称三角阵的指定特征值的二分法,将二分法与三对角化技巧相结合,就可得到求任意实对称阵指定特征值和对应特征向量的数值方法;

​设

\[T = \left( \begin{matrix} \alpha_1 & \beta_2\\ \beta_2 & \alpha_2 & \beta_3\\ & \ddots & \ddots & \ddots\\ & & \ddots & \ddots & \beta_n\\ & & & \beta_n & \alpha_n \end{matrix} \right) \]

我们可以假定 \(\beta_i\neq 0\) ,即 \(T\) 不可约,否则可以进行分块处理;

​记 \(p_i(\lambda)\)\(T-\lambda I\)\(i\) 阶顺序主子式,则容易验证如下递推关系

\[p_0(\lambda) = 1,\quad p_1(\lambda) = \alpha_1 - \lambda\\ p_i(\lambda) = (\alpha_i-\lambda)p_{i-1}(\lambda) - \beta_i^2p_{i-2}(\lambda)\\ i = 2,\cdots, n \]

由实对称,我们知道 \(p_i(\lambda)\) 的根都是实根,并且有

定理 7.4.1 上述多项式 \(p_i(\lambda)\) 满足

  • \(\exist M>0,\quad \mathrm{s.t.}\quad \forall\lambda>M,\ p_i(-\lambda)>0\) ,并且 \(p_i(\lambda)\) 的符号为 \((-1)^i\)
  • 相邻两个多项式没有公共根
  • \(p_i(\mu)=0\) ,则 \(p_{i-1}(\mu)p_{i+1}(\mu)<0\)
  • \(p_i(\lambda)\) 的根全是单重的,且严格分隔 \(p_{i+1}(\lambda)\) 的根

证明 由顺序主子式的首项 \((-\lambda)^i\) 容易发现 (1) 成立;我们假设有公共根 \(\mu\) 使得 \(p_{i-1}(\mu) = p_i(\mu)= 0\) ,则

\[0 = p_{i}(\mu) = (\alpha_i-\lambda)p_{i-1}(\mu) - \beta_i^2p_{i-2}(\mu) = -\beta_{i}^2p_{i-2}(\mu) \]

由于 \(\beta_i\neq 0\) ,从而 \(p_{i-2}(\mu) = 0\) ,以此类推有 \(p_0(\mu) = 0\) ,矛盾;

​由上面的关系,若 \(p_i(\mu)=0\) ,则

\[p_{i+1}(\mu) = (\alpha_{i+1}-\lambda)p_{i}(\mu) - \beta_{i+1}^2p_{i-1}(\mu) = - \beta_{i+1}^2p_{i-1}(\mu) \]

因此它们异号;

​最后,我们用数学归纳法证明根的分隔。当 \(i=1\) ,显然 \(\alpha_1\)\(p_1(\lambda)\) 的根,而 \(p_2(\alpha_1) = -\beta_2^2<0\) ,之前得到 \(p_2(\lambda)\) 首项为 \(\lambda^2\) ,因此它在 \(\alpha_1\) 两边有两个根,这样就证明了 \(i=2\) 的情况成立;

​假设对 \(i=k\) 成立,设 \(p_{k-1}(\lambda),\ p_{k}(\lambda)\) 有严格分隔的根

\[\nu_1<\nu_2<\cdots<\nu_{k-1},\quad \mu_1<\mu_2<\cdots<\mu_k \]

满足

\[\mu_1<\nu_1<\mu_2<\nu_2<\cdots<\nu_{k-1}<\mu_k \]

根据递推公式有

\[p_{k+1}(\mu_j) = -\beta_{k+1}^2p_{k-1}(\mu_j),\quad j=1,\cdots,k \]

根据 \(p_{k-1}(\nu_j) = 0\) ,并且由全是单重根,它在根处变号,可以得到

\[(-1)^{j-1}p_{k-1}(\mu_j) > 0,\quad j=1,\cdots,k \]

于是我们有

\[(-1)^jp_{k+1}(\mu_j) >0,\quad j=1,\cdots,k \]

这意味着 \(p_{k+1}(\lambda)\) 在区间

\[(-\infty,\mu_1),(\mu_1,\mu_2),\cdots,(\mu_{k-1},\mu_k),(\mu_k,+\infty) \]

上都有根,根据这 \(n+1\) 个区间,每个区间上有一个根,即证.

 

​给定实数 \(\mu\) ,定义 \(s_k(\mu)\) 是数列 \(p_0(\mu),\cdots,p_k(\mu)\) 的变号数,并规定:若 \(p_i(\mu)=0\) ,则 \(p_i(\mu)\)\(p_{i-1}(\mu)\) 同号,例如

\[T = \left( \begin{matrix} 1 & 1 & 0\\ 1 & 1 & 1\\ 0 & 1 & 1 \end{matrix} \right) \]

则有

\[\begin{aligned} &p_0(\lambda) = 1,\\ &p_1(\lambda) = 1-\lambda,\\ &p_2(\lambda) = (1-\lambda)^2-1,\\ &p_3(\lambda) = (1-\lambda)^3-2(1-\lambda) \end{aligned} \]

于是当 \(\mu=1\) ,就有

\[p_0(1) = 1,\quad p_1(1) = 0,\quad p_2(1) = -1,\quad p_3(1) = 0 \]

这样一来,总共有 \(p_1, p_2\) 之间的一次变号, \(s_3(1) = 1\) .

 

定理 7.4.2\(T\) 为不可约对称三角阵时, \(s_k(\mu)\) 恰好是 \(p_k(\mu)\) 在区间 \((-\infty,\mu)\) 内根的个数.

 

推论 7.4.1\(T\) 为不可约对称三角阵,则 \(s_n(\mu)\) 恰好是该矩阵在区间 \((-\infty,\mu)\) 内特征值的个数.

 

​于是我们得到利用二分法求 \(T\) 的指定特征值的过程:设有特征值

\[\lambda_1<\lambda_2<\cdots<\lambda_n \]

则必有

\[|\lambda_i|\le \rho(T) \le \|T\|_{\infty} \]

现在假定要求 \(\lambda_m\) ,则可以取区间 \((l_0,u_0) = (-\|T\|_{\infty},\|T\|_{\infty})\) ,然后每次取中点 \(r\) ,计算 \(s_n(r)\) 判断在 \((-\infty,r)\) 上特征值的个数,选取含有 \(\lambda_m\) 的部分,直到最终区间长度达到要求;

​由于高阶多项式计算容易发生溢出,因此定义

\[q_i(\lambda) = \dfrac{p_i(\lambda)}{p_{i-1}(\lambda)},\quad i=1,\cdots,n \]

利用之前的递推关系有

\[q_1(\lambda) = p_1(\lambda) = \alpha_1-\lambda,\quad q_i(\lambda) = \alpha_i-\lambda-\dfrac{\beta_i^2}{q_{i-1}(\lambda)},\quad i=2,\cdots, n \]

容易发现 \(s_n(\mu)\) 恰为 \(q_1(\mu),\cdots,q_n(\mu)\) 中负数的个数;

​注意到 \(q_i(\lambda)\) 可能为 \(0\) ,按照变号数的规则,此时我们将其替换为一个小正数.

 

到这里关于数值代数的内容暂时告一段落,本文只整理了课程教材《数值线性代数》中的部分重要内容,有一些较难的部分在课程中没有展开教学,因此这部分的笔记只能等以后再温习课程时添加了。

posted @ 2022-02-19 20:16  Bluemultipl  阅读(1297)  评论(0编辑  收藏  举报