基本性質
定理 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\) ,按照變號數的規則,此時我們將其替換為一個小正數.
到這里關於數值代數的內容暫時告一段落,本文只整理了課程教材《數值線性代數》中的部分重要內容,有一些較難的部分在課程中沒有展開教學,因此這部分的筆記只能等以后再溫習課程時添加了。