数值分析 笔记 4

矩阵特征值的计算

特征值的概念

矩阵 A\boldsymbol{A} 的特征值是指满足：

Ax=λxx≠0\boldsymbol{Ax} = \lambda\boldsymbol{x}\quad \boldsymbol{x \neq 0}

的标量 λ\lambda

关于特征值有以下结论：

  1. det(A)=∏λi(A)\mathrm{det}(\boldsymbol{A}) = \prod\lambda_{i}(\boldsymbol{A})
  2. λ(A)=λ(AT)\lambda(\boldsymbol{A}) = \lambda(\boldsymbol{A}^{T})
  3. 对于对角阵或三角阵，特征值为对角元
  4. 对于分块对角阵或对角块为方阵的分块三角阵，特征值为对角块特征值的并集
  5. 相似矩阵的特征值相等，即 λ(A)=λ(V−1AV)\lambda(\boldsymbol{A}) = \lambda(\boldsymbol{V}^{-1}\boldsymbol{AV})
  6. 对矩阵做运算可以直接对应到特征值上

定义：
代数重数为 λ\lambda 作为方程 det⁡(λI−A)=0\det(\lambda\boldsymbol{I - A}) = 0 的解的重数，几何重数代表 λ\lambda 对应的特征子空间的维数，可以证明 几何重数不大于代数重数

且不同 λ\lambda 对应的特征向量都是线性无关的，但是不一定能张成整个空间（即几何重数之和可能小于 nn ），如果几何重数和几何重数恒定相等，这种矩阵称为 非亏损阵，反之为亏损阵，非亏损阵的 nn 个特征向量一定是全空间的一组基

有如下定理：

∃Xs.t.X−1AX=Λ⇔A非亏损\exist \boldsymbol{X}\quad \text{s.t.} \boldsymbol{X}^{-1}\boldsymbol{AX} = \boldsymbol{\Lambda} \Leftrightarrow \boldsymbol{A} \text{非亏损}

其中 Λ\boldsymbol{\Lambda} 是对角阵

自然可以证明实对称矩阵一定是非亏损的，并且其特征值一定是实数，可正交对角化，即 X\boldsymbol{X} 是正交矩阵

通过重数，可以定义强有力但是我不会的 Jordan 分解：

A=XJX−1J=[J1⋱Jp]Js=[λs1λs⋱⋱1λs]\boldsymbol{A} = \boldsymbol{XJX}^{-1} \quad \boldsymbol{J} = \begin{bmatrix} \boldsymbol{J}_{1} & & \\[0.6em] & \ddots & \\[0.6em] & & \boldsymbol{J}_{p} \\ \end{bmatrix} \quad \boldsymbol{J}_{s} = \begin{bmatrix} \lambda_{s} & 1 & & \\[0.6em] & \lambda_{s} & \ddots & \\[0.6em] & & \ddots & 1 \\[0.6em] & & & \lambda_{s} \\ \end{bmatrix}

其中每个 λ\lambda 对应几何重数个 Jordan 块，阶数之和为代数重数

特征值的分布

很多时候我们关心的是矩阵特征值的极值，例如谱半径是矩阵特征值绝对值的最大值，2-范数、2-条件数也与特征值有关

定义矩阵 A∈Cn×n\boldsymbol{A} \in \mathbb{C}^{n\times n} 的 Gerschgorin 圆盘为：圆心为 aiia_{ii}，半径 ri=(∑j=1n∣aij∣)−∣aii∣r_{i} = \Big(\sum\limits_{j = 1}^{n}|a_{ij}|\Big) - |a_{ii}|，则有：

A\boldsymbol{A} 的特征值一定在 nn 个圆盘的并集上

如果有 mm 个圆盘组成了独立的连通支，则这 mm 个圆盘中恰好有 mm 个特征值

幂法与反幂法

幂法的作用是为了求主特征值与主特征向量，即模最大的特征值与其对应的向量

幂法的定义为：取几乎任意的非零向量 v0\boldsymbol{v}_{0}，通过以下方式计算向量序列：

vk=Avk−1\boldsymbol{v}_{k} = \boldsymbol{Av}_{k - 1}

设主特征值为 λ1\lambda_{1}，其对应的特征向量为 x1\boldsymbol{x}_{1}，则有，当 A\boldsymbol{A} 有唯一主特征时：

lim⁡k→∞vk=λ1lim⁡k→∞vk−1=Cx1\lim\limits_{k\to\infty}\boldsymbol{v}_{k} = \lambda_{1}\lim\limits_{k\to\infty}\boldsymbol{v}_{k - 1} = C\boldsymbol{x}_{1}

我们针对非亏损阵 A\boldsymbol{A} 进行证明：

设 A\boldsymbol{A} 线性无关的特征向量为 x^1,⋯ ,x^n\hat{\boldsymbol{x}}_{1}, \cdots, \hat{\boldsymbol{x}}_{n}，v0=∑i=1nαix^i\boldsymbol{v}_{0} = \sum\limits_{i=1}^{n}\alpha_{i}\hat{\boldsymbol{x}}_{i}，其中 α1≠0\alpha_{1} \neq 0，这也是 几乎随机 的原因，则有：

vk=Akv0=∑i=1nαiλikx^i=λ1k[∑i=1sαix^i+∑i=s+1nαi(λiλ1)kx^i]\begin{align*}\boldsymbol{v}_{k} &= \boldsymbol{A}^{k}\boldsymbol{v}_{0} \\&= \sum\limits_{i=1}^{n}\alpha_{i}\lambda_{i}^{k}\hat{\boldsymbol{x}}_{i} \\&= \lambda_{1}^{k}\Bigg[\sum\limits_{i=1}^{s}\alpha_{i}\hat{\boldsymbol{x}}_{i} + \sum\limits_{i=s+1}^{n}\alpha_{i}(\frac{\lambda_{i}}{\lambda_{1}})^{k}\hat{\boldsymbol{x}}_{i}\Bigg]\end{align*}

其中，ss 为 λ1\lambda_{1} 的代数重数

由于 ∣λi/λ1∣<1|\lambda_i / \lambda_1 | < 1，因此当 kk 充分大的时候：

lim⁡k→∞vk≈λ1k(∑i=1sαix^i)\lim\limits_{k\to\infty}\boldsymbol{v}_{k} \approx \lambda_{1}^{k}\Big(\sum\limits_{i=1}^{s}\alpha_{i}\hat{\boldsymbol{x}}_{i}\Big)

证毕

规格化

幂法的主要问题在于，当 kk 很大的时候，可能在收敛之前就已经发生了溢出现象

解决溢出的一种方式为规格化，即在每一步计算完成之后除以向量的模最大的元素 max⁡‾(vk)\overline{\max}(\boldsymbol{v}_{k})，即：

uk=vkmax⁡‾(vk)vk+1=Auk\begin{align*} \boldsymbol{u}_{k} &= \frac{\boldsymbol{v}_{k}}{\overline{\max}(\boldsymbol{v}_{k})} \\ \boldsymbol{v}_{k + 1} &= \boldsymbol{Au}_{k}\\ \end{align*}

这样的话有：

lim⁡k→∞uk=x1max⁡‾(x1)lim⁡k→∞∣∣vk∣∣∞=λ1\begin{align*} \lim\limits_{k\to\infty}\boldsymbol{u}_{k} &= \frac{\boldsymbol{x}_{1}}{\overline{\max}(\boldsymbol{x}_{1})} \\ \lim\limits_{k\to\infty}||\boldsymbol{v}_{k}||_{\infty} &= \lambda_{1}\\ \end{align*}

其中 x1\boldsymbol{x}_{1} 是某个主特征向量

幂法加速

  1. 原点位移：
    幂法的收敛速度取决于次大特征值与最大值的模长之比，而特征值的运算是线性的，因此考虑 B=A−sI\boldsymbol{B} = \boldsymbol{A} - s\boldsymbol{I}，则可以选取合适的 ss，使得 λ1−s\lambda_{1} - s 仍然是 B\boldsymbol{B} 的主特征值，且：

    ∣λ2(B)λ1−s∣\Big|\frac{\lambda_{2}(\boldsymbol{B})}{\lambda_{1} - s}\Big|

    尽可能小，其中 λ2(B)\lambda_{2}(\boldsymbol{B}) 代表模次大的特征值

  2. Relay 商加速：
    实对称矩阵 A\boldsymbol{A} 对应向量 x≠0\boldsymbol{x} \neq \boldsymbol{0} 的 Relay 商定义为：

    R(x)=⟨Ax,x⟩⟨x,x⟩=xTAxxTxR(\boldsymbol{x}) = \frac{\lang\boldsymbol{Ax}, \boldsymbol{x}\rang}{\lang\boldsymbol{x}, \boldsymbol{x}\rang} = \frac{\boldsymbol{x}^{T}\boldsymbol{Ax}}{\boldsymbol{x}^{T}\boldsymbol{x}}

    可以证明：

    λn≤R(x)≤λ1\lambda_{n} \leq R(\boldsymbol{x}) \leq \lambda_{1}

    其中 λ1,λn\lambda_{1}, \lambda_{n} 为最大和最小特征值

    可以证明规格化幂法满足：

    R(uk)=λ1+O((λ2λ1)2k)max⁡‾(uk)=λ1+O((λ2λ1)k−1)\begin{align*} R(\boldsymbol{u}_{k}) &= \lambda_{1} + O\Big(\Big(\frac{\lambda_{2}}{\lambda_{1}}\Big)^{2k}\Big) \\ \overline{\max}(\boldsymbol{u}_{k}) &= \lambda_{1} + O\Big(\Big(\frac{\lambda_{2}}{\lambda_{1}}\Big)^{k - 1}\Big) \end{align*}

    因此 Relay 商可以比最大值更快收敛

反幂法

即通过求 A−1\boldsymbol{A}^{-1} 模最小的特征值，取倒数后得到 A\boldsymbol{A} 的主特征值

QR 分解

Householder 变换

对于 w∈Rn\boldsymbol{w} \in \mathbb{R}^{n}, wTw=1\boldsymbol{w}^T\boldsymbol{w} = 1，称 H(w)=I−2wwT\boldsymbol{H(w)} = \boldsymbol{I} - 2\boldsymbol{ww}^{T} 为 Householder 矩阵，Hx\boldsymbol{Hx} 被称为 Householder 变换

定义 v=wwTx=(wTx)w\boldsymbol{v} = \boldsymbol{ww}^{T}\boldsymbol{x} = (\boldsymbol{w}^{T}\boldsymbol{x})\boldsymbol{w} ，几何上来看 v\boldsymbol{v} 为 x\boldsymbol{x} 在 w\boldsymbol{w} 方向上的投影

显然 H\boldsymbol{H} 是正交矩阵，因此 Householder 变换不会改变向量的二范数，并且可以证明以下定理：

对于 x≠y∈Rn\boldsymbol{x} \neq \boldsymbol{y} \in \mathbb{R}^{n}, 且 ∣∣x∣∣2=∣∣y∣∣2||\boldsymbol{x}||_{2} = ||\boldsymbol{y}||_{2}，则存在 Householder 矩阵 H\boldsymbol{H} 使得 Hx=y\boldsymbol{Hx = y}

取 y=−σe1\boldsymbol{y} = -\sigma\boldsymbol{e}_{1}，其中 σ=sgn(x1)∣∣x∣∣2\sigma = \text{sgn}(x_1)||\boldsymbol{x}||_{2}，则构造 H\boldsymbol{H} 时，满足 v=x+σe1,w=v/∣∣v∣∣2\boldsymbol{v = x} + \sigma\boldsymbol{e_{1}}, \boldsymbol{w = v} / \boldsymbol{||v||_{2}}

这实际上和矩阵消元的过程是一样的，而且由于 H\boldsymbol{H} 是正交的，因此可以连续左乘 householder 矩阵，以达到正交三角化的效果

矩阵的 QR 分解

可以证明，任意实矩阵都一定存在 QR 分解，即 A=QR\boldsymbol{A = QR}，其中 Q\boldsymbol{Q} 为正交阵，如果 A\boldsymbol{A} 非奇异且 R\boldsymbol{R} 对角元全正，则分解唯一

下面利用 Householder 变换构造 QR 分解

令：

Hk⋯H2H1A=R\boldsymbol{H_{k}\cdots\boldsymbol{H}_{2}\boldsymbol{H}_{1}A = R}

则

A=H1H2⋯HkR=QRA = \boldsymbol{H_{1}\boldsymbol{H}_{2}\cdots\boldsymbol{H}_{k}R = QR}

QR 分解的直接构造是简单的，大致思路为：

  1. 构造 H1\boldsymbol{H}_{1} 给第一列消元
  2. 考虑右下角的 (n−1)×(n−1)(n - 1) \times (n - 1) 阶子阵，构造H2′\boldsymbol{H}_{2}' 用于消除子阵第一列，对应原矩阵中的第二列，之后即有：

    H2=[I10T0H2′]\boldsymbol{H_{2}} = \begin{bmatrix} \boldsymbol{I}_{1} & \boldsymbol{0}^{T} \\ \boldsymbol{0} & \boldsymbol{H}_{2}' \end{bmatrix}

  3. 以此类推即可

可以使用算法构造每一个 H\boldsymbol{H} 对应的 v\boldsymbol{v}，需要的乘法次数为 O(mn2)O(mn^{2})，可减少为 nm2−n3/3nm^{2} - n^{3} / 3，如果生成 Q\boldsymbol{Q} 则需要额外 O(mn2)O(mn^{2})

Givens 旋转变换

正交变换的另一种方式，通过旋转矩阵来做到

二维平面上的旋转矩阵为：

G=[cos⁡θsin⁡θ−sin⁡θcos⁡θ]\boldsymbol{G} = \begin{bmatrix} \cos \theta & \sin\theta \\ -\sin\theta & \cos\theta \end{bmatrix}

这代表顺时针旋转 θ\theta 角

如果想旋转 nn 维 空间中的两个元素xi,xjx_i, x_j，则只需要将这个矩阵嵌入 I\boldsymbol{I} 的 {i,j}×{i,j}\{i, j\}\times \{i, j\} 这四个位置即可

如果取：

G=1α[x1x2−x2x1],x12+x22=α2\boldsymbol{G} = \frac{1}{\alpha}\begin{bmatrix} x_1 & x_2 \\ -x_2 & x_1 \end{bmatrix},\quad x_1^2 + x_2^2 = \alpha^2

则

G[x1x2]=[α0]\boldsymbol{G}\begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} \alpha \\ 0 \end{bmatrix}

这也就实现了消元，并且嵌入到 nn 维之后仍然每次只影响两行，可以通过不断旋转第 (1,i)(1, i) 两个元素来给向量消元，之后按照 QR 分解同样的算法可以构造另一种 QR

QR 分解的迭代算法

仍然是矩阵求特征值的问题，考虑我们已知 A\boldsymbol{A} 的一个特征向量 x\boldsymbol{x}，并且通过正交变换将其消元：Hx=σe1\boldsymbol{Hx} = \sigma\boldsymbol{e}_{1}，则：

HAHTe1=HA1σx=1σHλx=λσσe1=λe1\boldsymbol{HAH}^{T}\boldsymbol{e}_{1} = \boldsymbol{HA}\frac{1}{\sigma}\boldsymbol{x} = \frac{1}{\sigma}\boldsymbol{H}\lambda\boldsymbol{x} = \frac{\lambda}{\sigma}\sigma\boldsymbol{e}_{1} = \lambda\boldsymbol{e}_{1}

根据单位向量的特性，可以得到：

HAHT=[λrT0A1]\boldsymbol{HAH}^{T} = \begin{bmatrix} \lambda & \boldsymbol{r}^{T} \\ \boldsymbol{0} & \boldsymbol{A}_{1} \end{bmatrix}

于是我们只需要求 (n−1)×(n−1)(n - 1) \times (n - 1) 维矩阵的特征值即可，这种技术被称为 收缩技术

迭代算法

首先介绍实 Schur 分解：

对于任意实方阵 A\boldsymbol{A}，存在正交阵 Q\boldsymbol{Q}，使得 QTAQ=S\boldsymbol{Q}^{T}\boldsymbol{AQ} = \boldsymbol{S} 为实 Schur 型，或称为拟上三角阵

实 Schur 型是指对角块为 1 阶或 2 阶 的分块上三角阵，并且满足 2 阶分块阵的特征值是 A\boldsymbol{A} 的两个共轭复特征值

则利用 QR 分解求特征值的一种迭代算法为：

Ak=QkRkAk+1=RkQk\begin{align*} \boldsymbol{A}_{k} &= \boldsymbol{Q}_{k}\boldsymbol{R}_{k} \\ \boldsymbol{A}_{k + 1} &= \boldsymbol{R}_{k}\boldsymbol{Q}_{k} \\ \end{align*}

可以证明：若矩阵 A\boldsymbol{A} 的等模的特征值为实重特征值或复共轭特征值, 则QR迭代所得矩阵序列 {Ak}\{\boldsymbol{A}_{k}\} 基本收敛于拟上三角阵，而拟上三角阵是方便求特征值的

对于实对称阵，其迭代的极限也为对角阵

实用的 QR 算法

QR 算法的不足之处：

  1. 计算量很大（化简为上 Hessenberg 阵）
  2. 手链可能很慢（带原点位移的 QR）

上 Hessenberg 阵是指仅有上三角、主对角线旁边的副对角线有元素的矩阵，如下：

[1234567898765]\begin{bmatrix} 1& 2& 3& 4 \\ 5& 6& 7& 8 \\ & 9& 8& 7 \\ & & 6& 5 \end{bmatrix}

对于上 Hessenberg 型矩阵，可以应用 Givens 旋转的 QR 分解，每步的计算量为 O(n2)O(n^{2})，而普通矩阵 QR 分解的计算量为 O(n3)O(n^{3})，并且可以证明迭代过程中仍然是保持 Hessenberg 的

对于任意矩阵，都可以用 Householder 变换将其正交相似为上 Hessenberg 矩阵，思路是将矩阵二分块为 1×11 \times 1 与 (n−1)×(n−1)(n - 1)\times(n - 1)，这样只用将第一列第 2∼n2\sim n 行进行消元，即得到了 hessenberg 的形式

具体来说：

用Householder变正交相似约化为上Hessenberg矩阵
用Householder变正交相似约化为上Hessenberg矩阵

而实对称阵的带原点位移迭代是指：

Ak=QkRk−skIAk+1=RkQk+skI\begin{align*} \boldsymbol{A}_{k} &= \boldsymbol{Q}_{k}\boldsymbol{R}_{k} - s_k\boldsymbol{I} \\ \boldsymbol{A}_{k + 1} &= \boldsymbol{R}_{k}\boldsymbol{Q}_{k} + s_k\boldsymbol{I} \\ \end{align*}

如果是非对称阵，则可能有复特征值，因此需要做二维平面上的平移

许可协议

本文采用 署名-非商业性使用-相同方式共享 4.0 国际 许可协议，转载请注明出处。