数值线性代数：Arnoldi求解特征值/特征向量

原创已于 2025-08-17 10:14:26 修改 · 3k 阅读

8 ·

CC 4.0 BY-SA版权

JiNan YouQuan Software Co. Ltd.

文章标签：

#其他

于 2023-07-12 22:18:28 首次发布

CAx 同时被 2 个专栏收录

122 篇文章

订阅专栏

HPC

27 篇文章

订阅专栏

文章介绍了数值线性代数中的主要问题，如线性方程组和特征值求解，并重点讨论了用于解决大型稀疏矩阵特征值问题的ARPACK库中的IRAM算法。该算法基于Arnoldi过程，通过隐式重启来保持正交性，以逼近矩阵的部分特征值和特征向量。IRAM通过修正的QR分解来剔除非目标特征值，适用于大规模科学计算和工程问题。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

线性方程组求解、最小二乘法、特征值/特征向量求解是(数值)线性代数的主要研究内容。

在力学、气象学、电磁学、金融等学科中，许多问题最终都归结为特征值、特征向量的求解。

ARPACK使用IRAM(Implicit Restarted Arnoldi Method)求解大规模稀疏矩阵的部分特征值与特征向量。了解或者熟悉IRAM算法，必定有助于更好地使用ARPACK中相关特征值求解函数。

本文拟就ARPACK中特征值求解的IRAM算法进行分析，希望对从事相关研究的朋友们有所帮助。

注1：限于研究水平，分析难免不当，欢迎批评指正。

零、预修

对于n阶级复矩阵 $\boldsymbol{A}\in \mathbb{C}^{n\times n}$ ，若存在n阶向量 $\boldsymbol{x}\in \mathbb{C}^{n}$ 与标量 $\lambda \in \mathbb{C}^{1}$ ，满足 $\boldsymbol{\mathbf{A}}\boldsymbol{\mathbf{x}}=\lambda \boldsymbol{x}$ ，则称 $\lambda$ 是矩阵 $\boldsymbol{A}$ 的特征值， $\boldsymbol{x}$ 是对应的特征向量。

更一般地，若对于n阶级复矩阵 $\boldsymbol{A}\in \mathbb{C}^{n\times n}$ 、 $\boldsymbol{B}\in \mathbb{C}^{n\times n}$ ，若存在n阶向量 $\boldsymbol{x}\in \mathbb{C}^{n}$ 与标量 $\lambda \in \mathbb{C}^{1}$ ，满足 $\boldsymbol{\mathbf{A}}\boldsymbol{\mathbf{x}}=\lambda\boldsymbol{B} \boldsymbol{x}$ ，则称 $\lambda$ 是矩阵 $\boldsymbol{A}$ 的广义特征值， $\boldsymbol{x}$ 是对应的广义特征向量。

若 $\boldsymbol{Q}\in \mathbb{C}^{n\times n}$ ，满足 $\boldsymbol{Q}^{H}\boldsymbol{Q}=\boldsymbol{Q}\boldsymbol{Q}^{H}=\boldsymbol{I}$ ，则称 $\boldsymbol{Q}$ 为酉矩阵。

定义Hessenberg矩阵 $\boldsymbol{H}$ ，

$\boldsymbol{H}=\begin{pmatrix} h_{1,1} & h_{1,2} & \cdots & h_{1,k}\\ h _{2,1} & h_{2,2} & \cdots &h_{2,k} \\ 0& \ddots& \ddots &\vdots \\ 0& 0 & h_{k,k-1}& h _{k,k} \end{pmatrix}$

Schur分解定理：设 $\boldsymbol{A}\in \mathbb{C}^{n\times n}$ ，则存在酉矩阵 $\boldsymbol{U}\in \mathbb{C}^{n\times n}$ ，使得

$\boldsymbol{U}^{H}\boldsymbol{A}\boldsymbol{U}=\boldsymbol{T}$

其中， $\boldsymbol{T}$ 是上三角阵。

实Schur分解：设 $\boldsymbol{A}\in \mathbb{R}^{n\times n}$ ，则存在正交矩阵 $\boldsymbol{Q}\in \mathbb{R}^{n\times n}$ ，使得

$\boldsymbol{Q}^{T}\boldsymbol{A}\boldsymbol{Q}=\begin{pmatrix} \boldsymbol{R}_{11} & \boldsymbol{R}_{12} & \cdots & \boldsymbol{R}_{1m}\\ & \boldsymbol{R}_{22}& \cdots & \boldsymbol{R}_{2m}\\ & & \ddots & \vdots \\ & & & \boldsymbol{R}_{mm} \end{pmatrix}$

其中， $\boldsymbol{R}_{ii}$ 是实数，或者是一个具有一对复共轭特征值的2阶仿真。

为了叙述方便，本文就标准特征值问题进行叙述，而不讨论广义特征值问题。

Gram-Schmidt正交化：

一、大型稀疏矩阵特征值算法

1.1 Projection Method

投影法在较小的线性空间内寻找满足精度要求的近似解，也即在某个线性空间内寻找真解的投影。这其实就是投影法得名的原因。

对于广义特征值问题 $\boldsymbol{A}\boldsymbol{x}=\lambda \boldsymbol{B}\boldsymbol{x}$ ， $\boldsymbol{A}\in \mathbb{C}^{n\times n}$ ， $\boldsymbol{B}\in \mathbb{C}^{n\times n}$ ，

首先，构造列满秩矩阵 $\boldsymbol{\mathcal{K}}\in \mathbb{C}^{n\times m}$ 与 $\boldsymbol{\mathcal{L}}\in \mathbb{C}^{n\times m}$ ，其中 $m\leq n$ ；然后，设真解 $\boldsymbol{x}$ 在线性空间 $\boldsymbol{\mathcal{K}}$ 的投影为 $\boldsymbol{\tilde{x}}$ ，即 $\boldsymbol{x}=\boldsymbol{\mathcal{K}}\boldsymbol{\tilde{x}}$ ，令其满足Petrov-Galerkin条件，即 $\forall \boldsymbol{y}\in \boldsymbol{\mathcal{L}}$ ，均有 $\boldsymbol{\mathcal{L}}^{T}\left (\boldsymbol{A}-\lambda \boldsymbol{B}\right )\boldsymbol{\mathcal{K}}\boldsymbol{\tilde{x}} =\boldsymbol{0}$ 。 $\boldsymbol{\mathcal{K}}$ 称为搜索空间， $\boldsymbol{\mathcal{L}}$ 称为约束空间。

若 $\boldsymbol{\mathcal{L}}=\boldsymbol{\mathcal{K}}$ 时，则有 $\boldsymbol{\mathcal{K}}^{T}\left (\boldsymbol{A}-\lambda \boldsymbol{B}\right )\boldsymbol{\mathcal{K}}\boldsymbol{\tilde{x}}=\boldsymbol{0}$ 。

Ref. from SLEPc Users Manual

Most eigensolvers provided by slepc perform a Rayleigh-Ritz projec tion for extracting the spectral approximations, that is, they project the problem onto a low dimensional subspace that is built appropriately.

Starting from this general idea, eigensolvers differ from each other in which subspace is used, how it is built and other technicalities aimed at improving convergence, reducing storage re quirements, etc.

1.2 Power Method

1.3 Lancos Method

1.4 Arnoldi Method

Arnoldi在1951年发表的文章中提出了一种获取大型稀疏矩阵的部分特征值、特征向量的方法。在这种方法中，通过Ritz\Galerkin余量法，将矩阵 $\boldsymbol{u}$ 约化为低阶Hessenberg矩阵，进而求解该Hessenberg矩阵的特征值与特征向量，用以近似原矩阵 $\boldsymbol{u}$ 部分特征值、特征向量，即

$\left [ \kappa _{i} \right ]\left ( \lambda \boldsymbol{I}-\boldsymbol{u} \right )\left [ k_{i} \right ]\boldsymbol{c}=0$

上式也可以写作

$\boldsymbol{A}\boldsymbol{V}=\boldsymbol{V}\boldsymbol{H}+\boldsymbol{f}\boldsymbol{e}_{k}^{T}$

其中， $\boldsymbol{V}=\left ( \boldsymbol{v}_{1}, \boldsymbol{v}_{2},\cdots ,\boldsymbol{v}_{k} \right )\in \mathbb{C}^{n\times k }$ ， $\boldsymbol{V}^{H}\boldsymbol{V}=\mathbf{I}$ ， $\boldsymbol{H}=h_{i,j}\in \mathbb{C}^{k\times k }$ ， $\boldsymbol{f}\in \mathbb{C}^{n }$ ， $\boldsymbol{V}^{H} \boldsymbol{f}=\boldsymbol{0}$ ， $e_{k}^{T}=\left ( 0,0,\cdots,0,1 \right )\in \mathbb{C}^{k}$ 。

对于 $k+1$ , 设

$\boldsymbol{A}\boldsymbol{\tilde{V}}=\boldsymbol{\tilde{V}}\boldsymbol{\tilde{H}}+\boldsymbol{\tilde{f}}\boldsymbol{e}_{k+1}^{T}$

其中 $\boldsymbol{\tilde{V}}=\left ( \boldsymbol{V},\boldsymbol{v}_{k+1} \right )\in \mathbb{C}^{n\times \left (k+1 \right ) }$ ， $\mathbf{\tilde{H}}=\begin{pmatrix} \boldsymbol{H} & \mathbf{h}_{2}\\ \mathbf{h}_{1}& h_{3} \end{pmatrix}\in \mathbb{C}^{\left (k+1 \right )\times \left (k+1 \right ) }$ ， $\boldsymbol{\tilde{f}}\in \mathbb{C}^{n }$ ， $\boldsymbol{\tilde{V}}^{H}\boldsymbol{\tilde{f}}=0$ ， $e_{k+1}^{T}=\left ( 0,0,\cdots,0,1 \right )\in \mathbb{C}^{k+1}$ ，则有

$\left\{\begin{matrix} \boldsymbol{A}\boldsymbol{\tilde{V}}=\left ( \boldsymbol{A}\boldsymbol{V},\boldsymbol{A}\boldsymbol{v}_{k+1} \right )\\ \boldsymbol{\tilde{V}}\boldsymbol{\tilde{H}}=\left ( \boldsymbol{V},\boldsymbol{v}_{k+1} \right ) \begin{pmatrix} \boldsymbol{H} & \boldsymbol{h}_{2}\\ \boldsymbol{h}_{1}& h_{3} \end{pmatrix} \end{matrix}\right.\Rightarrow \left\{\begin{matrix} \boldsymbol{A}\boldsymbol{V}=\boldsymbol{V}\boldsymbol{H}+\boldsymbol{v}_{k+1}\boldsymbol{h}_{1}\\ \boldsymbol{A}\boldsymbol{v}_{k+1}=\boldsymbol{V}\boldsymbol{h}_{2}+\boldsymbol{v}_{k+1}h_{3}+\boldsymbol{\tilde{f}} \end{matrix}\right.$

若取 $\boldsymbol{v}_{k+1}=\frac{\boldsymbol{f}}{\left \| \boldsymbol{f} \right \|}$ ，则有

$\left\{\begin{matrix} \boldsymbol{h}_{1}=\left \| \mathbf{f} \right \|\boldsymbol{e}_{k}^{T} \\ \boldsymbol{h}_{2}=\frac{\boldsymbol{V}^{H}\boldsymbol{A}\boldsymbol{f}}{\left \| \boldsymbol{f} \right \|}\\ h_{3}=\frac{\boldsymbol{f}^{H}\boldsymbol{A}\boldsymbol{f}}{\left \| \boldsymbol{f} \right \|^{2}}\end{matrix}\right.$ ， $\left\{\begin{matrix} \boldsymbol{\tilde{V}}=\left ( \boldsymbol{V}, \frac{\boldsymbol{f}}{ \left \| \boldsymbol{f} \right \| }\right )\\ \boldsymbol{\tilde{H}}=\begin{pmatrix} \boldsymbol{H} &\frac{\boldsymbol{V}^{H}\boldsymbol{A}\boldsymbol{f}}{\left \| \boldsymbol{f} \right \|} \\ \left \| \mathbf{f} \right \|\boldsymbol{e}_{k}^{T} &\frac{\boldsymbol{f}^{H}\boldsymbol{A}\boldsymbol{f}}{\left \| \boldsymbol{f} \right \|^{2}} \end{pmatrix}\\ \boldsymbol{\tilde{f}}= \frac{\boldsymbol{A}\boldsymbol{f}}{\left \| \boldsymbol{f} \right \|}-\boldsymbol{V}\frac{\boldsymbol{V}^{H}\boldsymbol{A}\boldsymbol{f}}{\left \| \boldsymbol{f} \right \|}-\frac{\boldsymbol{f}^{H}\boldsymbol{A}\boldsymbol{f}}{\left \| \boldsymbol{f} \right \|^{3}}\boldsymbol{f}\end{matrix}\right.$

另外可以看出 $\boldsymbol{v}_{k+1}=\frac{\boldsymbol{A}\boldsymbol{v}_{k}-\boldsymbol{V}\boldsymbol{h}\left ( :,k \right )}{\left \| \boldsymbol{f} \right \|}$ ，即 $\boldsymbol{V}=\left ( \boldsymbol{v}_{1}, \boldsymbol{v}_{2},\cdots ,\boldsymbol{v}_{k} \right )$ 是Krylov子空间 $\mathcal{K}\left ( \boldsymbol{A},\mathbf{v}_{1},k \right )=\left ( \mathbf{v}_{1} ,\boldsymbol{A}\mathbf{v}_{1},\boldsymbol{A}^{2}\mathbf{v}_{1},\cdots ,\boldsymbol{A}^{k-1}\mathbf{v}_{1}\right )$ 的一组标准正交基。

在实际数值计算过程中，由于计算精度等问题，通常很难保证 $\boldsymbol{V}^{H}\boldsymbol{V}=\boldsymbol{I}$ ，也就是说在矩阵 $\boldsymbol{V}=\left ( \boldsymbol{v}_{1}, \boldsymbol{v}_{2},\cdots ,\boldsymbol{v}_{k} \right )$ 逐列计算过程中，存在正交性丢失的问题。

针对此问题，Sorensen在1992的论文中提出了一种称之IRAM(Implicit Restarted Arnoldi Method)的算法。Sorensen认为误差 $\boldsymbol{f}$ 实际上是初始向量 $\boldsymbol{v}_{1}$ 的函数，可以通过不断选取新的初始向量 $\boldsymbol{v}_{i}^{new}=\boldsymbol{A}\boldsymbol{v}_{k+1}$ ，迭代Arnoldi过程 $p$ 次来强迫 $\boldsymbol{f}\rightarrow \mathbf{0}$ ，并且给出了收敛性证明。

$\boldsymbol{A}\boldsymbol{V}_{k+p}=\boldsymbol{V}_{k+p}\boldsymbol{H}_{k+p}+\boldsymbol{f}_{k+p}\boldsymbol{e}_{k+p}^{T}$

对于得到的 $\boldsymbol{H}_{k+p}$ ，需要施加 $p$ 次带偏移的QR分解，以剔除不想要的特征值/特征向量。这实际上就是ARPACK使用的IRAM(Implicit Restarted Arnoldi Method)算法。

ARPACK中，采用了统一的特征值表述，

$\mathbf{OP}\cdot \boldsymbol{x}=w \boldsymbol{x}$

对于实特征值问题提供了5种计算模式，而对复特征值问题提供了3种计算模式，

Real	Mode	Formulation	Eigenvalue	Description	Applications
	Mode 1	$\boldsymbol{A}\boldsymbol{x}=\lambda \boldsymbol{x}$ $\mathbf{OP} = \boldsymbol{A}$ $\boldsymbol{B} = \boldsymbol{I}$ $w=\lambda$	$\lambda =w$	$\boldsymbol{A}$ : symmetric
	Mode 2	$\boldsymbol{A}\boldsymbol{x}=\lambda \boldsymbol{M}\boldsymbol{x}$ $\mathbf{OP} = \boldsymbol{M}^{-1}\boldsymbol{A}$ $\boldsymbol{B} = \boldsymbol{M}$ $w=\lambda$	$\lambda =w$	$\boldsymbol{A}$ : symmetric $\boldsymbol{M}$ : symmetric positive definite
	Mode 3 Shift-and-Invert mode	$\boldsymbol{K}\boldsymbol{x}=\lambda \boldsymbol{M}\boldsymbol{x}$ $\mathbf{OP} = \left ( \boldsymbol{K}-\sigma \boldsymbol{M}\right )^{-1}\boldsymbol{M}$ $\boldsymbol{B}=\boldsymbol{M}$ $w=\frac{1}{\lambda -\sigma }$	$\lambda=\frac{1}{w }+\sigma$ $w\neq 0$	$\boldsymbol{K}$ : symmetric $\boldsymbol{M}$ :symmetric positive semi-definite	Mode Analysis Mode Solver
	Mode 4 Buckling mode	$\boldsymbol{K}\boldsymbol{x}=\lambda \boldsymbol{K}_{G}\boldsymbol{x}$ $\mathbf{OP}=\left ( \boldsymbol{K}-\sigma \boldsymbol{K}_{G}\right )^{-1}\boldsymbol{K}$ $\boldsymbol{B}=\boldsymbol{K}$ $w=\frac{\lambda }{\lambda-\sigma }=1+\frac{\sigma }{\lambda-\sigma }$ $\sigma \neq 0$	$\lambda =\frac{\sigma w}{w-1}$ $w\neq 1$	$\boldsymbol{K}$ : symmetric positive semi-definite $\boldsymbol{K}_{G}$ : symmetric indefinite	Buckling Analysis
	Mode 5 Cayley transformed mode	$\boldsymbol{A}\boldsymbol{x}=\lambda \boldsymbol{M}\boldsymbol{x}$ $\mathbf{OP}=\left ( \boldsymbol{A}-\sigma \boldsymbol{M}\right )^{-1}\left ( \boldsymbol{A}+ \sigma \boldsymbol{M}\right )$ $\boldsymbol{B}=\boldsymbol{M}$ $w =\frac{ \lambda+\sigma}{\lambda-\sigma}==1+\frac{2\sigma }{\lambda-\sigma }$	$\lambda =\sigma\frac{ w+1}{w-1}$ $w\neq 1$	$\boldsymbol{A}$ : symmetric $\boldsymbol{M}$ :symmetric positive semi-definite
Complex	Mode 1	$\boldsymbol{A}\boldsymbol{x}=\lambda \boldsymbol{x}$ $\mathbf{OP} = \boldsymbol{A}$ $\boldsymbol{B} = \boldsymbol{I}$ $w=\lambda$	$\lambda =w$
	Mode 2	$\boldsymbol{A}\boldsymbol{x}=\lambda \boldsymbol{M}\boldsymbol{x}$ $\mathbf{OP}= \boldsymbol{M}^{-1}\boldsymbol{A}$ $\boldsymbol{B} = \boldsymbol{M}$ $w=\lambda$	$\lambda =w$	$\boldsymbol{M}$ :hermitian positive definite
	Mode 3 shift-and-invert mode	$\boldsymbol{A}\boldsymbol{x}=\lambda \boldsymbol{M}\boldsymbol{x}$ $\mathbf{OP} = \left ( \boldsymbol{A}-\sigma \boldsymbol{M}\right )^{-1}\boldsymbol{M}$ $\boldsymbol{B}=\boldsymbol{M}$ $w=\frac{1}{\lambda -\sigma }$	$\lambda=\frac{1}{w }+\sigma$ $w\neq 0$	$\boldsymbol{M}$ :hermitian semi-definite	Mode Analysis Mode Solver

1.4.1 Mode 3

对于ARPACK中的Mode 3, 即有

$\boldsymbol{K}\boldsymbol{x}=\lambda \boldsymbol{M}\boldsymbol{x}\Rightarrow \left ( \boldsymbol{K}-\sigma \boldsymbol{M} \right )=\left (\lambda -\sigma \right )\boldsymbol{M}\boldsymbol{x}$

整理可得，

$\left ( \boldsymbol{K}-\sigma \boldsymbol{M} \right )^{-1}\boldsymbol{M}\boldsymbol{x}=\frac{1}{\lambda -\sigma }\boldsymbol{x}$

若针对上式得到的特征值为 $w$ ，即 $w=\frac{1}{\lambda -\sigma }$ ，则有 $\lambda=\frac{1}{w }+\sigma$ 。

1.4.2 Mode 4

对于ARPACK中的Mode 4, 即有

$\boldsymbol{K}\boldsymbol{x}=\lambda \boldsymbol{K}_{G}\boldsymbol{x}\Rightarrow \boldsymbol{K}_{G}\boldsymbol{x}=\frac{1}{\lambda }\boldsymbol{K}\boldsymbol{x}\Rightarrow \left (\boldsymbol{K}_{G} -\frac{1}{\sigma } \boldsymbol{K}\right )\boldsymbol{x}=\left ( \frac{1}{\lambda } -\frac{1}{\sigma }\right )\boldsymbol{K}\boldsymbol{x}$

整理上式，则有

$\left ( \sigma \boldsymbol{K}_{G}-\boldsymbol{K} \right )\boldsymbol{x}=\frac{\sigma -\lambda }{\lambda }\boldsymbol{K}\boldsymbol{x}\Rightarrow \left ( \boldsymbol{K}-\sigma \boldsymbol{K}_{G} \right )^{-1}\boldsymbol{K}\boldsymbol{x}=\frac{\lambda }{\lambda-\sigma }\boldsymbol{x}$

若针对上式得到的特征值为 $w$ ，即 $w=\frac{\lambda }{\lambda-\sigma }$ ，则有 $\lambda =\frac{\sigma w}{w-1}$ 。

以上推导，可由ARPACK/SRC/dseupd.f代码可以得到，

         if (type .eq. 'SHIFTI') then
            do 40 k=1, ncv
               workl(ihd+k-1) = one / workl(ihd+k-1) + sigma
  40        continue
         else if (type .eq. 'BUCKLE') then
            do 50 k=1, ncv
               workl(ihd+k-1) = sigma * workl(ihd+k-1) /
     &                          (workl(ihd+k-1) - one)
  50        continue
         else if (type .eq. 'CAYLEY') then
            do 60 k=1, ncv
               workl(ihd+k-1) = sigma * (workl(ihd+k-1) + one) /
     &                          (workl(ihd+k-1) - one)
  60        continue
         end if

1.4.3 Mode 4*

~~可以看到，Mode 4需要满足 $\sigma \neq 0$ ，实际应用不太方便，可以考虑Mode 3来计算Mode 4问题，~~

~~$\boldsymbol{K}\boldsymbol{x}=\lambda \boldsymbol{K}_{G}\boldsymbol{x}\Rightarrow \boldsymbol{K}_{G}\boldsymbol{x}=\frac{1}{\lambda }\boldsymbol{K}\boldsymbol{x}$~~

$\left (\boldsymbol{K}_{G} -\sigma \boldsymbol{K}\right )\boldsymbol{x}=\left ( \frac{1}{\lambda } -\sigma\right )\boldsymbol{K}\boldsymbol{x}\Rightarrow\left (\boldsymbol{K}_{G} -\sigma \boldsymbol{K}\right )^{-1}\boldsymbol{K}=\frac{1}{\frac{1}{\lambda } -\sigma}\boldsymbol{x}$

$\left (\boldsymbol{K}_{G} -\sigma \boldsymbol{K}\right )^{-1}\boldsymbol{K}=\frac{1}{\frac{1}{\lambda } -\sigma}\boldsymbol{x}\Rightarrow\left (\boldsymbol{K}_{G} -\sigma \boldsymbol{K}\right )^{-1}\boldsymbol{K}=\frac{\lambda}{1 -\sigma\lambda}\boldsymbol{x}$

~~若针对上式得到的特征值为 $w$ ，即 $w=\frac{\lambda}{1 -\sigma\lambda}$ ，则有 $\lambda =\frac{w}{1+\sigma w }$ 。~~

~~更近一步，若 $\sigma=0$ ，则有 $\lambda =w$ 。~~

参考书籍

Golub G H , Loan C F V .Matrix Computations.Johns Hopkins University Press,1996.

Ford W .Numerical Linear Algebra with Applications using MATLAB. 2014.

徐树方. 数值线性代数(第二版). 北京大学出版社, 2010.

参考文献

Golub,G.H.a Eigenvalue Computation in the 20thcentury.J.Comput. Appl.Math.,123(1-2):35-65.

Lanczos C .An Iteration Method for the solution of the Eigenvalue Problem of Linear Differential and Integral Operators.Journal of research of the National Bureau of Standards, 1950, 45(4).

Arnoldi W E .The principle of minimized iterations in the solution of the matrix eigenvalue problem.Quarterly of Applied Mathematics, 1951, 9(1).

Lehoucq R B , Sorensen D C , Yang C .ARPACK Users' Guide: Solution of Large Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods.

D,C,Sorensen.Implicit Application of Polynomial Filters in a k-Step Arnoldi Method.Siam Journal on Matrix Analysis & Applications, 1992.

Burrus C S , Cox S J , Radke R J .A Matlab Implementation of the Implicitly Restarted Arnoldi Method for Solving Large-Scale Eigenvalue Problems[J]. 2000.

R.B. Lehoucq, Analysis and Implementation of an Implicitly Restarted Arnoldi Iteration, Rice University Technical Report TR95-13, Department of Computational and Applied Mathematics.

Hernandez V. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems. ACM Transactions on Mathematical Software, 2005, 31(3):351-362.

Jose E. Roman. SLEPc Users Manual Scalable Library for Eigenvalue Problem Computations. 2025.

Zhang Hong. CS 595 - Advanced Scientific Computing.