跳转至正文
GZ Cloudhome Logo

用于快速求解大型稀疏矩阵特征值的算法:Arnoldi 分解

发布于:2023 年 7 月 7 日 at 22:28
更新于:2023 年 7 月 11 日 at 19:08

最近在做科研的时候需要快速地解一个大型稀疏矩阵最小的几个特征值和对应的特征向量,经过文献调研发现 Arnoldi iteration 方法可以应用到这上面来。我花了一些时间了解 Arnoldi iteration 以及与其相关的 Lanczos iteration、Krylov 子空间,在这里整理一下。

从幂法讲起

我们知道,求解一个 n×nn \times n 矩阵 A\mathbf{A} 的最大特征值和对应的特征向量可以使用幂法(power method)。具体地,我们迭代地计算并归一化序列

x,Ax,A2x,\boldsymbol{x}, \mathbf{A} \boldsymbol{x}, \mathbf{A}^2 \boldsymbol{x}, \dots

在第 kk 步迭代时,我们可以得到 x(k)=Akx/Akx\boldsymbol{x}^{(k)} = \mathbf{A}^k \boldsymbol{x} / \| \mathbf{A}^k \boldsymbol{x} \|,随着 kk 的增大,我们得到的 x(k)\boldsymbol{x}^{(k)} 也将收敛于最大特征值对应的特征向量。

听起来很好,对吧?但是仔细一想,我们进行了多次迭代,最终只用到了最后的结果来得到输出,前面的所有迭代是否浪费了?如果利用前面的迭代,我们对特征向量的估计能否更加准确,或者说具有输出多个特征值对应的特征向量的近似值的能力?

想要回答这个问题,就需要来研究序列 x,Ax,A2x,\boldsymbol{x}, \mathbf{A} \boldsymbol{x}, \mathbf{A}^2 \boldsymbol{x}, \dots。事实上,这个序列在数学上有其对应的名字,这个序列的前 rr 个向量张成 Krylov 子空间(假设 r<=mr <= m

Kr(x)=Kr(x,A)span{x,Ax,A2x,,Ar1x}=R(Kr(x))Fn.\mathcal{K}^r(\boldsymbol{x}) = \mathcal{K}^r(\boldsymbol{x}, \mathbf{A}) \triangleq \mathrm{span} \Big\{ \boldsymbol{x}, \mathbf{A} \boldsymbol{x}, \mathbf{A}^2 \boldsymbol{x}, \cdots, \mathbf{A}^{r-1} \boldsymbol{x} \Big\} = \mathcal{R}(\mathcal{K}^r(\boldsymbol{x})) \subset \mathbb{F}^n.

下面我们看如何基于这个子空间来近似求解矩阵 A\mathbf{A} 的特征值和特征向量。

Rayleigh-Ritz 法求特征值的近似

我们可以利用 Gram-Schmidt 正交化方法提取上述子空间的一组单位正交基 VFn×r\mathbf{V} \in \mathbb{F}^{n \times r},然后即可得到 A\mathbf{A} 在基 V\mathbf{V} 下的投影矩阵

T=VAVFj×j,\mathbf{T} = \mathbf{V}^* \mathbf{A} \mathbf{V} \in \mathbb{F}^{j \times j},

T\mathbf{T} 做特征值分解,可以得到

TY=YΛT,\mathbf{T} \mathbf{Y} = \mathbf{Y} \mathbf{\Lambda_T},

其中 Y\mathbf{Y} 的各列 y1,,yr\boldsymbol{y}_1, \dots, \boldsymbol{y}_rT\mathbf{T} 的特征向量,对应 ΛT\mathbf{\Lambda_T} 的对角线元素(特征值)。

能否用 T\mathbf{T} 的特征值来近似代替 A\mathbf{A} 的特征值?幸运的是,答案是肯定的,这样的方法称作 Rayleigh-Ritz 方法。设 T\mathbf{T} 的特征值和特征向量分别为 μi\mu_iyi\boldsymbol{y}_i1ir1 \le i \le r),那么我们可以用 λ~i=μi\tilde{\lambda}_i = \mu_i 来近似表示 A\mathbf{A} 的第 ii 个特征值,用 x~i=Vyi\tilde{\boldsymbol{x}}_i = \mathbf{V} \boldsymbol{y}_i 来表示该特征值对应的特征向量。

如果投影向量 V\mathbf{V} 恰巧是由 A\mathbf{A}rr 个特征向量 x1,x2,,xr\boldsymbol{x}_1, \boldsymbol{x}_2, \dots, \boldsymbol{x}_r 张成的空间的单位正交基,那么可以证明求得的 λ~i\tilde{\lambda}_ix~i\tilde{\boldsymbol{x}}_iA\mathbf{A} 的特征值和特征向量完全相同。此时存在单位正交矩阵 PFr×r\mathbf{P} \in \mathbb{F}^{r \times r},使得

V=XP,\mathbf{V} = \mathbf{X} \mathbf{P},

因此

T=VAV=PXAXP=PΛAP,\begin{align*} \mathbf{T} = \mathbf{V}^* \mathbf{A} \mathbf{V} &= \mathbf{P}^* \mathbf{X}^* \mathbf{AXP} \\ &= \mathbf{P}^* \mathbf{\Lambda_A} \mathbf{P}, \end{align*}

因此 T\mathbf{T}A\mathbf{A} 具有相同的特征值,即 μi\mu_i 也是 A\mathbf{A} 的特征值。进一步,

yiVAVyi=yiVAVyiVyi2=μi,\boldsymbol{y}_i^* \mathbf{V}^* \mathbf{A} \mathbf{V} \boldsymbol{y}_i = \frac{\boldsymbol{y}_i^* \mathbf{V}^* \mathbf{A} \mathbf{V} \boldsymbol{y}_i}{\| \mathbf{V} \boldsymbol{y}_i \|^2} = \mu_i,

因此 x~i=Vyi\tilde{\boldsymbol{x}}_i = \mathbf{V} \boldsymbol{y}_i 也是 A\mathbf{A} 的特征向量。

由此可以预见,近似特征值和特征向量的 Rayleigh-Ritz 方法的精度取决于 V\mathbf{V} 的精度,即 V\mathbf{V} 是否包含了我们想要的特征值和特征向量。幸运的是,Krylov 子空间近似为 A\mathbf{A} 的最大的几个特征值所张成,这让我们能够以较高的精度来近似。

Arnoldi / Lanczos 分解

Arnoldi 分解Lanczos 分解 都是计算 Krylov 子空间的单位正交基 V\mathbf{V} 的高效算法,尤其适用于大型稀疏的 A\mathbf{A}。其中,Lanczos 方法是 Arnoldi 方法的一个特例,当 A\mathbf{A} 是埃尔米特矩阵(Hermitian)时 Arnoldi 方法变为 Lanczos 方法,更加高效。

下面是 Arnoldi 分解的步骤:

  1. 选择一个随机的向量 v1\boldsymbol{v}_1,要求 v1=1\| \boldsymbol{v}_1 \| = 1
  2. 对于所有的 j=1,,r1j = 1, \dots, r-1
    1. wj+1=Avj\boldsymbol{w}_{j+1}' = \mathbf{A} \boldsymbol{v}_j
    2. 对于所有的 k=1,,jk = 1, \dots, j,令 Tk,j=vkwj+1T_{k,j} = \boldsymbol{v}_k^* \boldsymbol{w}_{j+1}'
    3. wj+1=wj+1k=1jTk,jvk\boldsymbol{w}_{j+1} = \boldsymbol{w}_{j+1}' - \sum_{k=1}^j T_{k,j} \boldsymbol{v}_k
    4. Tj+1,j=wj+1T_{j+1, j} = \| \boldsymbol{w}_{j+1} \|
    5. Tj+1,j0T_{j+1, j} \neq 0,令 vj+1=wj+1/Tj+1,j\boldsymbol{v}_{j+1} = \boldsymbol{w}_{j+1} / T_{j+1, j};否则,选择任意一个模长为 1 且对称于 v1,,vj\boldsymbol{v}_1, \dots, \boldsymbol{v}_j 的向量作为 vj+1\boldsymbol{v}_{j+1}

由此方法得到的各个 vi\boldsymbol{v}_i 构成单位正交基,且 T=VAV\mathbf{T} = \mathbf{V}^* \mathbf{AV}。通过 Arnoldi 分解得到的正交基 V\mathbf{V} 可以让进一步得到的矩阵 T\mathbf{T} 为上 Hessenberg 型1矩阵。当矩阵 A\mathbf{A} 为 Hermitian 时,T\mathbf{T} 进一步简化为三对角矩阵。Hessenberg 矩阵的特征值可以高效地求得,如使用 QR 分解;再加上 T\mathbf{T} 的大小仅为 r×rr \times r,当 rnr \ll n 时这种算法将会相当高效。

在实践中,我们往往发现通过 Arnoldi 方法得到的 Ritz 特征值能够接近 A\mathbf{A} 的部分特征值,一般来说特别是最大的几个特征值们。下图是一个摘自维基百科的 Arnoldi 方法求一个 400×400400 \times 400 大小的矩阵的特征值的动图示例,红点为 Ritz 特征值,黑点为实际特征值—— illustration_arnoldi_iteration

Implicitly Restarted Arnoldi / Lanczos Method

前面提到,Arnoldi 分解可以用来高效地近似部分矩阵的特征值和特征向量。但是,普通的 Arnoldi 方法可能会遇到数值问题,且精度不够。Implicitly Restarted Arnoldi Method(IRAM)算法可以看作是对其的改进(当输入矩阵为 Hermitian 时,变为 IRLM)。

简单来讲,IRAM 通过多次重启 Arnoldi 分解来让求得的特征值近似每次重启都变得更好。每次重启中,IRAM 都会尝试去掉特征值中不满足要求的部分。这样,经过多次重启之后,求得的特征值和特征向量将会收敛到一个稳定的值。

著名的 ARPACK 包即实现了 IRAM / IRLM 算法,在 Python 和 Matlab 中都能方便地调用。在 Python 中,可以通过 scipy.sparse.linalg.eigs / scipy.sparse.linalg.eigsh 来调用。下面以 eigsh 为例,讲解一下部分调用的参数:

复杂度分析

IRLM 的时间复杂度和 nn 呈线性关系。下图是我做的不同 nn 下对 A\mathbf{A} 运行 scipy.sparse.linalg.eigsh 的实验。横轴是 nn,纵轴是时间,单位秒。在不同的 nn 下,A\mathbf{A} 每行的非零元素均保持为 d50d \approx 50,且每次只求 r=100r=100 个最大特征值(即 scipy.sparse.lilnalg.eigshwhich 参数为 LM)。 arnoldi-time-vs-n 可以看到,当 nn 超过 50000 之后,nn 和时间之间的关系确实是线性关系。为了看得更明显,我还基于 nn 从 50000 到 200000 部分的数据做了一个单变量线性拟合(如橙色线段所示)。拟合得到的斜率为 0.00158,说明每增加一个样本,求特征向量的用时将会增加 1.58 ms。

此外,我还用一个例子对比了 numpy.linalg.eighscipy.sparse.linalg.eigsh 的时间差别,该例子中,矩阵 A\mathbf{A}10000×1000010000 \times 10000 方阵,平均每行 50 个非零元素。可以看到,做全分解、且未利用稀疏特性的 numpy.linalg.eigh 花费的时间要多出很多。

函数运行时间
numpy.linalg.eigh2min 19s
scipy.sparse.linalg.eigshk=6which='LM'26s

附录

Krylov 子空间

我们定义矩阵

Km(x)=Km(x,A)=[xAxA2xAm1x]\mathbf{K}^m (\boldsymbol{x}) = \mathbf{K}^m (\boldsymbol{x}, \mathbf{A}) = \big[ \begin{matrix} \boldsymbol{x} & \mathbf{A} \boldsymbol{x} & \mathbf{A}^2 \boldsymbol{x} & \cdots & \mathbf{A}^{m-1} \boldsymbol{x} \end{matrix} \big]

为 Krylov 矩阵。其实这正是迭代到第 mm 步时,我们把各个迭代得到的向量作为矩阵的各列得到的。同时,Krylov 矩阵的各列张成 Krylov (子)空间

Km(x)=Km(x,A)span{x,Ax,A2x,,Am1x}=R(Km(x))Fn.\mathcal{K}^m(\boldsymbol{x}) = \mathcal{K}^m(\boldsymbol{x}, \mathbf{A}) \triangleq \mathrm{span} \Big\{ \boldsymbol{x}, \mathbf{A} \boldsymbol{x}, \mathbf{A}^2 \boldsymbol{x}, \cdots, \mathbf{A}^{m-1} \boldsymbol{x} \Big\} = \mathcal{R}(\mathcal{K}^m(\boldsymbol{x})) \subset \mathbb{F}^n.

下面我们来探讨一下 Krylov 子空间的维度。首先 Krylov 子空间 Km(x)\mathcal{K}^m(\boldsymbol{x}) 的维度不可能超过 nn,因为它各列向量的长度为 nn

dim(Km(x))n.\dim(\mathcal{K}^m(\boldsymbol{x})) \le n.

这也就意味着矩阵 Kn+1(x)\mathbf{K}^{n+1}(\boldsymbol{x}) 的各列是线性相关的。从另外一个角度来说,如果 x\boldsymbol{x}A\mathbf{A} 的一个特征向量,那么 Ax=λx\mathbf{A} \boldsymbol{x} = \lambda \boldsymbol{x},这意味着 Km(x)=K1(x)  m>=1\mathcal{K}^m(\boldsymbol{x}) = \mathcal{K}^1(\boldsymbol{x}) ~~ \forall m >= 1。初始的 x\boldsymbol{x} 必然由张满 Fn\mathbb{F}^n 的各个特征向量线性组合而成,因此 x\boldsymbol{x} 成分中有几个特征向量,Km(x)\mathbf{K}^m(\boldsymbol{x}) 各列便(这里的简单分析不考虑广义特征向量情形)最多张成几维空间。所以,可以预见,对于任意 xRn\boldsymbol{x} \in \mathbb{R}^n 存在一个最小的 mm1mn1 \le m \le n,使得

K1(x)K2(x)Km(x)=Km+1(x)\mathcal{K}^1(\boldsymbol{x}) \subsetneqq \mathcal{K}^2(\boldsymbol{x}) \subsetneqq \dots \subsetneqq \mathcal{K}^m(\boldsymbol{x}) = \mathcal{K}^{m+1}(\boldsymbol{x})

成立。对于这个 mmKm+1(x)Fn×(m+1)\mathbf{K}^{m+1}(\boldsymbol{x}) \in \mathbb{F}^{n \times (m+1)} 开始有线性相关的列。即存在 aFm+1\boldsymbol{a} \in \mathbb{F}^{m+1},使得

Km+1(x)a=0\mathbf{K}^{m+1}(\boldsymbol{x}) \boldsymbol{a} = \boldsymbol{0}

成立,即

p(A)x=amAmx++a1Ax+a0x=0.p(\mathbf{A}) \boldsymbol{x} = a_m \mathbf{A}^m \boldsymbol{x} + \dots + a_1 \mathbf{A} \boldsymbol{x} + a_0 \boldsymbol{x} = \boldsymbol{0}.

其中

p(λ)=amλm++a1λ+a0.p(\lambda) = a_m \lambda^m + \dots + a_1 \lambda + a_0.

由于从 mm 开始矩阵的各列才开始线性相关,所以

am0.a_m \neq 0.

我们称 p(λ)p(\lambda)A\mathbf{A} 相对于 x\boldsymbol{x} 的最小多项式(我之前博客中关于矩阵的最小多项式的介绍见这里)。

Galerkin 条件和 Ritz 对

设有 Fn\mathbb{F}^n 的某 rr 维子空间 V\mathcal{V},一个向量 vV\boldsymbol{v} \in \mathcal{V} 满足 Galerkin 条件

w,Avλv=0,  wV\langle \boldsymbol{w}, \mathbf{A} \boldsymbol{v} - \lambda \boldsymbol{v} \rangle = 0, ~~ \forall \boldsymbol{w} \in \mathcal{V}

时,我们说向量 v\boldsymbol{v}λ\lambda 构成 Ritz 对,其中 λ\lambda 为 Ritz 值,v\boldsymbol{v} 为 Ritz 向量。

V\mathcal{V} 的某组单位正交基构成矩阵 WFn×r\mathbf{W} \in \mathbb{F}^{n \times r},定义

GWAWFr×r,\mathbf{G} \triangleq \mathbf{W}^* \mathbf{A} \mathbf{W} \in \mathbb{F}^{r \times r},

G\mathbf{G} 的特征值和特征向量构成对 si,θi,i=1,,r\boldsymbol{s}_i, \theta_i, i = 1, \dots, r,我们有

Gsi=θisi.\mathbf{G} \boldsymbol{s}_i = \theta_i \boldsymbol{s}_i.

viWsi\boldsymbol{v}_i \triangleq \mathbf{W} \boldsymbol{s}_i,我们可以看到

Aviθivi=AWsiθiWsi.\mathbf{A} \boldsymbol{v}_i - \theta_i \boldsymbol{v}_i = \mathbf{A} \mathbf{W} \boldsymbol{s}_i - \theta_i \mathbf{W} \boldsymbol{s}_i.

进一步,

Gsi=θisi    W(Aviθivi)=0,\mathbf{G} \boldsymbol{s}_i = \theta_i \boldsymbol{s}_i ~~ \Rightarrow ~~ \mathbf{W}^* \big( \mathbf{A} \boldsymbol{v}_i - \theta_i \boldsymbol{v}_i \big) = \boldsymbol{0},

Range(W)=V\mathrm{Range}(\mathbf{W}) = \mathcal{V},所以

w,Aviθivi=0,  wV.\langle \boldsymbol{w}, \mathbf{A} \boldsymbol{v}_i - \theta_i \boldsymbol{v}_i \rangle = 0, ~~ \forall \boldsymbol{w} \in \mathcal{V}.

所以我们得出结论:当且仅当 s,θ\boldsymbol{s}, \thetaG\mathbf{G} 的一对特征向量和特征值时,vWs,λθ\boldsymbol{v} \triangleq \mathbf{W} \boldsymbol{s}, \lambda \triangleq \thetaA\mathbf{A} 的一对 Ritz 向量和 Ritz 值,即 Ritz 对。

Householder 变换

Householder 变换可以理解成(实数域下,复数域下几何性质有所不同)给定一个超平面,把输入向量沿着这个超平面做一个对称,得到的向量作为输出。

具体地,给定一个超平面的单位法向量 v\boldsymbol{v},线性变换

xx2v(vx)\boldsymbol{x} \mapsto \boldsymbol{x} - 2 \boldsymbol{v} (\boldsymbol{v}^* \boldsymbol{x})

即为 Householder 变换。此线性变换对应矩阵

P=I2vv.\mathbf{P} = \mathbf{I} - 2 \boldsymbol{v v}^*.

容易看出,矩阵 P\mathbf{P} 是酉(unitary)矩阵且 Hermitian 的,满足下列性质:

  1. P=P1\mathbf{P}^* = \mathbf{P}^{-1}
  2. P=P\mathbf{P}^* = \mathbf{P}
  3. P1=P\mathbf{P}^{-1} = \mathbf{P},由 1 和 2 导出。

Householder 变换可以把一个矩阵转换为 Hessenberg 型。见 维基百科相关词条

Schur 分解

设矩阵 ACn×n\mathbf{A} \in \mathbb{C}^{n \times n},存在一个酉矩阵 Q\mathbf{Q} 和一个上三角矩阵 R\mathbf{R},使得

AQ=QR.\mathbf{AQ} = \mathbf{QR}.

其中 R\mathbf{R} 的对角元为 A\mathbf{A} 的特征值。

Footnotes

  1. 即黑森贝格矩阵。它与三角阵很相似。一个上黑森伯格矩阵 H\mathbf{H} 的次对角元以下的所有元素都为 0,即 Hij=0,i>j+1H_{ij} = 0, i > j+1。一个下黑森伯格矩阵的次对角元以上的所有元素都为 0。