对称特征值问题的计算方法 您所在的位置:网站首页 实对称矩阵的特征值求解技巧有哪些 对称特征值问题的计算方法

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

2023-12-17 19:56| 来源: 网络整理| 查看: 265

目录基本性质对称 QR 方法三对角化隐式对称 QR 迭代Jacobi 方法经典 Jacobi 方法过关 Jacobi 方法二分法

 

基本性质

定理 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)


【本文地址】

公司简介

联系我们

今日新闻

    推荐新闻

    专题文章
      CopyRight 2018-2019 实验室设备网 版权所有