多项式特征值问题的数值方法

求解大型复杂结构特征值问题的Lanczos分布式并行算法研究
IV
图表清单
图 1.1自由度弹簧-质点系统 (3)
图 3.1单质量弹簧系统 (22)
图 3.2 n自由度阻尼的弹簧-质点系统 (25)
图 3.3 GJD和QJD方法的残差范数随迭代次数增加的变化 (26)
图 3.4 GJD和QJD方法的残差范数随迭代次数增加的变化 (27)
表 2.1 Case 1的结果和比较12
表 2.2 Case 2的结果和比较 (12)
表 2.3 Case 3的结果和比较 (13)
表 2.4 Case 4的结果和比较 (14)
表 3.1 CPUtime比较 (26)
表 3.2 CPUtime比较 (28)
表 3.3特征值比较 (28)
承诺书
本人声明所呈交的硕士学位论文是本人在导师指导下进行的研究工作及取得的研究成果。

除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得南京航空航天大学或其他教育机构的学位或证书而使用过的材料。

本人授权南京航空航天大学可以将学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。

（保密的学位论文在解密后适用本承诺书）
作者签名：
日 期：
南京航空航天大学硕士学位论文
1
第一章 绪论
1.1多项式特征值问题的来源
我们考虑m 次的矩阵多项式（或λ-矩阵）
110(,),m m m m P A A A A λλλ−−=+++L （1.1） 其中,0:n n
k A k m ×∈=£。

多项式特征值问题（PEP ）是要找到一个特征值λ和相应的非
零特征向量x 满足
(,)0.P A x λ=
m=1的情况对应于广义特征值问题（GEP ）
A x
B x
λ=
并且如果0A I =，我们得到标准特征值问题（SEP ）
.A x x λ= （1.2）
另一个重要的情况是当m=2时，这时就是二次特征值问题（QEP ）[1]。

假定(,)P A λ是一个如（1.1）所示的m 次的n n ×阶矩阵多项式
定义1.1 当(,)P A λ不恒等于0时(,)P A λ是正则的。

我们规定整篇文章中都有此假定，即(,)P A λ是正则的。

PEP 问题是要找出标量λ和非零向量x 和y 满足
(,)0P A x λ=，*(,)0y P A λ=
λ称作特征值，x 和y 是对应的右和左特征向量。

等价地，特征值是下面特征多项式的根
det((,))0P A λ=
假定d 是这个标量多项式det((,))P A λ的次数。

det((,))P A λ的d 个根称为有限特征值。

如果d
m n <，那么我们说(,)P A λ有mn d −个无限特征值。

假定
m A 是奇异的。

我们考
虑关于(,)P A μ的PEP 问题，其中1μλ=。

我们看到0μ=是这个PEP 问题的一个特征值，并
且特征向量是生成null （A ）的向量。

关于(,1)m
P A λλ的PEP 问题的零特征值与关于(,)P A λ的
PEP 问题的无限特征值相对应。

求解大型复杂结构特征值问题的Lanczos 分布式并行算法研究
2
定义1.2 假定λ是(,)P A λ的一个特征值。

λ的代数重数（由p 表示）是λ作为特征多项式det((,))P A λ的一个根的重数。

λ的几何重数（由q 表示）是张成null (,)P A λ的线性无关向量的数目。

假定p 是λ的代数重数，
而q 是λ的几何重数。

如果p=1，则λ是简单特征值。

当1p >时，λ是一个重特征值。

当1p >且q p =时，λ是半单的。

否则，λ就是一个亏损特征值。

假定1
01(,,,)()
m m n A A A A +=∈ΜK £。

我们定义齐次矩阵多项式(,,)P A αβ如下
(,,)m
k m k
k k P A A αβαβ
−==
∑
即(,,)P A αβ是m 次齐次矩阵多项式，其中2
(,)αβ∈£。

我们假设是正则的，即对于所有(,)(0,0)αβ≠，det((,,))P A αβ不恒为零。

该齐次多项式特征值问题（PEP ）是要找到一对标量(,)(0,0)αβ≠和非零向量,n
x y ∈£满足
*(,,)0,(,,)0P A x y P A αβαβ==
向量x ，y 称为对应于特征值(,)αβ的右和左特征向量。

因此，一个特征值现在是
det((,,))0P A αβ=的解中任何一条通过2£上原点的线。

对于0β≠，我们定义λαβ=。

从而我们将非齐次矩阵多项式(,)P A λ通过以下式子与
齐次矩阵多项式建立起联系
(,)(,,)m P A P A βλαβ=
我们看到，求解齐次PEP 问题等价于求解非齐次PEP 问题。

例如，对于齐次广义特征值问题，众所周知，我们有形式
()0A B x βα−=
PEPs 问题的重要性体现在求解科学和工程问题中它所起到的各种作用。

我们简要列出一些例子。

QEPs 问题以及更普遍的PEPs 问题在各种问题中有着广泛的应用。

有着自然产生PEPs 问题的许多例子。

一些物理现象可由一个具有如下矩阵系数的二阶常微分方程（ODE ）建模
(),M z D z Kz f t ++=g g g
（1.3） (0),z a = （1.4） (0).z b =g
（1.5）
南京航空航天大学硕士学位论文
3
齐次方程的解的形式为t e u λ，其中u 是一个常向量。

这引出了QEP 问题
2()0.M D K u λλ++= （1.6）
图 1.1自由度弹簧-质点系统
一个著名的例子是阻尼弹簧-质点系统。

在图1.1中，我们考虑2自由度的弹簧-质点阻尼系统。

在一些假设下，该系统的动态由一个形式为（1.3）-（1.5）的ODE 控制。

在这种情况下，
1122(,,,)z x y x y =表示质量1m 和2m 的坐标，1122(,,,)M m m m m =是质量矩阵，12462357(,,,)D diag d d d d d d d d =++++是阻尼矩阵，12462357(,,,)K diag k k k k k k k k =++++是
刚度矩阵，其中0,0(17)i i d k i >>≤≤。

QEPs 问题出现在结构力学、控制理论、流体力学以及其他更具体的应用，我们可以参考Tissue 和Meerbergen 的概括论述[1]。

关于更高阶的PEPs 的有趣的实际例子在[2]中有给出。

1.2 多项式特征值问题数值方法概述
GEPs 问题的数值求解的标准方法是：将相关矩阵约简为某个揭示特征值的更简单的形式，例如一个对(,)A B 的广义Schur 形式。

不幸的是，这些规范形式不能推广到次数大于1的λ−矩阵。

这是PEPs 问题的数值求解的一个主要困难。

对于求解PEPs 问题的数值方法
[37]
−，有两种主要的划分。

1. 可划分为在原始形式上解决问题的方法，以及把它线性化为一个具有更大维数的
GEP 问题，然后应用GEP 问题的方法进行求解的方法。

2. 可划分为对于稠密的、小到中型问题的方法，以及对于大型稀疏问题的迭代方法。

求解大型复杂结构特征值问题的Lanczos 分布式并行算法研究
4
求解二次特征值问题（QEP ）的一般方法
[812]
−是将其线性化为广义特征值问题（GEP ），
QEP 问题转化成GEP 问题后，可以借助GEP 问题的数值方法求解QEP 问题，但是这样一来，所要求解的问题的维数就扩大为原始问题的两倍，存储量和计算量明显增加，并且QEP 问题中矩阵的很多良好性质比如对称性、正定性、和稀疏性等可能不会在GEP 问题中得到体现，这是线性化求解PEP 问题的一个缺点。

对于大型特征值问题，由于受到计算机存储量和计算速度等条件的限制，使得对中小型特征值问题行之有效的方法，对大型特征值问题变得无能为力，这就急需我们去探讨对于大型特征值问题有效的数值方法。

有时人们所关心的仅仅是一小部分特征值，这就不必求出大型系统的所有特征值。

而投影方法就是这样一种只求一部分特征值的方法，例如多项式特征值问题的广义投影算法
[1318]
−，这类方法构建一个子空间k κ，它的一组正交基为列组成的矩阵为k V ，将
PEP 问题投影到子空间k κ上，然后去求解小规模的投影问题*
(,)0k k V P A V z λ=。

本论文研究的多项式特征值问题的shift-and-invert Arnoldi 方法和Jacobi-Davidson 方法都属于子空间投影方法。

1.3本文的主要研究内容
本课题主要研究多项式特征值问题的一些子空间投影方法，包括shift-and-invert Arnoldi 方法和Jacobi-Davidson 方法。

后面各章内容如下：
第二章研究了多项式特征值问题的shift-and-invert Arnoldi 方法，它基于Arnoldi 算法的直接应用，并结合了shift-and-invert 技术，是一种有效的快速收敛的直接投影方法。

重点研究了CEP 问题的shift-and-invert Arnoldi 方法，给出了CEP 问题的迭代shift-and-invert Arnoldi 过程的算法和数值试验。

第三章研究了多项式特征值问题的Jacobi-Davidson 方法，以二次特征值问题为例进行研究，包括线性化Jacobi-Davidson 方法和直接Jacobi-Davidson 方法。

给出了实际问题的数值算例来比较这两种方法。

第四章给出总结与展望。

南京航空航天大学硕士学位论文
5
第二章 求解多项式特征值问题的shift-and-invert Arnoldi 方法
2.1预备知识
2.1.1 shift-and-invert Arnoldi 方法简介
shift-and-invert Arnoldi 方法
[19]
是一种子空间投影方法，这种方法基于对单一矩阵1
A B −使
用标准Arnoldi 算法产生的Krylov 子空间。

本章节我们重点研究三次特征值问题（CEP ）的shift-and-invert Arnoldi 方法，提出一个迭代的shift-and-invert Arnoldi 算法，也讨论了不精确的shift-and-invert 的情况，并给出数值算例验证了这个新方法。

2.1.2三次特征值问题（CEP ）
大型三次特征值问题（CEP ）
32()()0L x A B C D x λλλλ=+++= （2.1.1）
（其中,,,n n
A B C D ×∈¡
且n
x ∈¡
）出现于许多科学工程应用中，例如，三维（3D ）Schrodinger
方程描述的，具有非抛物线形能带结构的半导体量子点模型可以由有限差分法产生一个阶数高达211400的CEP 问题；见[4,20]。

目前已经开发出一些方法来求解大型CEP 问题。

一种直接方法是使用所谓的“线性化”，将CEP 问题转化为一个等价的线性特征值问题，然后使用标准Krylov 子空间投影方法（如Arnoldi 算法）[1]。

然而，线性化使得问题的阶数增至三倍，即线性化问题作用的空间的维数是原始空间维数的三倍，这些对于大型问题，会极大增加计算和存储成本。

这使得我们对直接投影方法产生兴趣，在直接投影方法中，构造一个子空间，并且将投影直接应用到CEP 问题()L λ。

特别地，Ye 考虑了QEP 的情况
[19]。

这里，我们考虑CEP 的情况。

构
造出某个子空间的正交基k Q ，并且我们通过例如正交投影1
()()T
k k k L Q A L Q λλ−=来近似()L λ的极端特征值。

这种类型的方法有残量迭代方法[21]
，Jacobi-Davidson 方法
[622]
，和Krylov-type
子空间方法
[18]。

尽管这些方法使用一个类似的投影过程，但是它们所构造的应用于投影的子空
间是不同的。

我们注意到，残量迭代和Jacobi-Davidson 方法需要使用shift-and-invert 矩阵1
()L σ−。

下一章我们将会对Jacobi-Davidson 方法进行研究。

求解大型复杂结构特征值问题的Lanczos 分布式并行算法研究
6
我们仅仅考虑一种直接投影方法，它是使用1
A B −产生的Krylov 子空间的。

Krylov 子空间的基可以通过标准Arnoldi 或者Lanczos 算法构造（见[7,23,24]），但是这样的一个方法可能显得过于简单，这是因为此Krylov 子空间没有使用矩阵C 和D 来构造，而除非C 和D 是相对小的矩阵或者B 的小扰动，这时才可以期望不会产生影响。

然而，本文的一个主要发现是，当这个直接投影与shift-and-invert 相结合时，C 矩阵的影响在shift-and-invert 转化后将会降低，并且完全由1
A B −构造的Krylov 子空间也会足以选取好的近似。

我们将会提出一个迭代的shift-and-invert Arnoldi 算法作为基础部分。

由于这个算法将会需要反置一个shifted 矩阵，我们也会讨论shift-and-invert 方程组通过例如预处理迭代方法来不精确的求解的情况。

我们将会给出数值算例来验证这个算法。

2.2直接Arnoldi 过程
广义特征值问题（GEP ）()0A B x λ+=为Krylov 子空间
()
1
1111111(,){,,,}k k A B q span q A Bq A B q −−−−Κ=L
构造了一组正交基12,,k q q q L。

它同时产生了
Hessenberg 矩阵()1
,[]T
k k
k
i j H Q A B Q
h −==作
为1
A B −在Krylov 子空间上的投影，其中12[,,]k k Q q q q =L 。

我们有
()1
1,1T
k
k k k k k k A B Q
Q H h q e −++=+ （2.2.1）
其中n
k e ∈¡
表示I 的第k 列。

如果0(,)(1)u u θ=是k H −的一个特征对，那么称为Ritz 值和
Ritz 向量的0(,)u θ和k x Q u =可以用作1
A B −−的一个近似特征对。

这个近似的残差由
()01,1k k k k A B x h u Aq θ+++=
给出，其中12[,,]T k u u u u =L ，1,k k h +衡量Krylov 子空间的1
A B −-不变性，它随着k 增加而变小。

此外，对应于极端特征值的特征向量的最后分量k u 也变小；例如见[7]。

这导致0(,)x θ的收敛。

我们提议也对QEP 问题使用上述Arnoldi 过程。

特别地，对于CEP 32()L A B C D λλλλ=+++，我们对1A B −应用Arnoldi 算法，产生1
A B −的投影k H 。

然后我们在这个Krylov 子空间上建立1A C −和1A D −的投影
11,T T k k k k k k G Q A CQ F Q A DQ −−==.
注意到，k G 和k F 通常是全矩阵。

然后我们通过这个投影CEP 问题来近似特征值问题()L λ
南京航空航天大学硕士学位论文
7
32()k k k k L I H G F λλλλ=+++.
如果0θ是()k L λ的一个特征值，12[,,]T
k u u u u =L 是一个单位特征向量，即有
()320
000k k k I H G F u θ
θθ+++=
那么使用0(,)x θ（其中k x Q u =）作为()L λ的一个近似特征对，它们再一次被称为Ritz 值和Ritz 向量。

Ritz 对的近似质量通过检查近似的残差来确定，这些残差表示为
()320003211100032111000321100001,132********()()()()k k k k k k T
k k k k k k k k k k k k k k k k k k k k k k r A B C D x A I A B A C A D Q u A Q A BQ A CQ A DQ u
A Q Q H A CQ A DQ u h Aq e u
A Q Q H Q G Q F Q G Q F A CQ A DQ u h θθθθθθθθθθθθθθθθθθθ−−−−−−−−++−−=+++=+++=+++=++++=+++−−+++1,1320000001,101,10001,1001,1()()()()()k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k
Au q A Q Q H Q G Q F u CQ DQ AQ G AQ F u h Au q h Au q CQ DQ AQ G AQ F u h Au q CQ AQ G u DQ AQ F u h Au q u θθθθθθθθθθθθ++++++++++=+++++−−+=++−−=+−+−=+Δ （2.2.2）
其中0()()k k k k k k k CQ AQ G DQ AQ F θΔ=−+−。

和Arnoldi 算法一样，期望第一项随着k 增加而减小。

另一方面，除非在某些特殊的情况，例如C 和D 非常小的情况，第二项中的k Δ不必小。

因此，当第一项变得足够小时，
1100k k k k k k k r CQ AQ G DQ AQ F A A C A A D θθ−−≈−+−≤+
注意到，k Δ是空间k spanQ 的1A C −-不变性和1A D −-不变性的一个量度。

因此，上述的一个直接应用将会有一个k u Δ水平的停滞残量。

这是因为构造的Krylov 子空间没有保持来自1A C −和1A D −的信息，并且投影k G 和k F 可能没有捕获到足够的信息。

这些困难可以通过结合一个shift-and-invert 策略降低k Δ项的影响来克服。

2.3 CEP 问题的迭代的shift-and-invert Arnoldi 过程
众所周知，shift-and-invert 转换是加速标准特征值问题的Krylov 子空间方法收敛的一项非常有效的技术。

它通常也用作二次特征值问题的加速。

它也可以直接用于以直接投影方法构造一个子空间
[21,22]。

这里我们研究一种迭代过程，它结合了shift-and-invert 技术与第2节中Arnoldi
算法的直接应用，去得出对于CEP 问题的一种有效的算法。

求解大型复杂结构特征值问题的Lanczos 分布式并行算法研究
8
考虑计算接近某个σ（或离σ的距离不超过η）的()L λ的的特征值，0λσ=是一个初始近似特征值，0x 是一个近似特征向量。

使用shift-and-invert 转换()'01
1μλλλ==−，即
''0,1λλλμλ=+=得到
'3'2'000'3'22'32000000'3'22'0000()()()()(3)(32)(3)(32)()
L A B C D
A A
B A B
C A B C
D A A B A B C L λλλλλλλλλλλλλλλλλλλλλλλ=++++++=+++++++++=++++++ 然后我们寻找这个shift-and-invert 问题
$µµµµ332()()L L A B
C D λμλμμμ==+++ 的最大值（绝对值最大），其中µµµµ2
0000
(),(32),3,A L B A B C C A B D A λλλλ==++=+=。

现在，应用第2节的直接Arnoldi 过程到$()L μ，并且通过$3
2
()k k k k
L I H G F μμμμ=+++来近似$()L μ，其中k H 是Arnoldi 算法产生的µµ1A B −的投影，并且k G 和k F 是µµ1A C −和µµ1
A D
−的投影，即
µµµµ11
,T
T k k
k k k k G Q A CQ F Q A DQ −−==
定理2.1 如果0μ是$()k L μ的最大值
（绝对值最大），12[,,]T
k u u u u =L 是对应的初始特征向量，设定1k x Q u =，1001λλμ=+，我们有
µ323311111001,110()()()k k k k k A B C D x h u Aq u λλλλλμλλ+++++≤−+−Δ （2.3.1）
其中µµµµ0()()k k k k k
k k CQ AQ G DQ AQ F μΔ=−+−。

证明：对于$()L
μ,由（3）我们有 $µµµµµ0101,10()[()()]k k k k k k k k
k k L x h u Aq CQ AQ G DQ AQ F u μμμ++=+−+− 现在，使用$31100()()()L L
λλλμ=−，我们得到 $µ3331110011001,110()()()()()k k k k k L x L x h u Aq u λλλμλλμλλ++=−=−+−Δ 其中µµµµ0()()k k k k k
k k CQ AQ G DQ AQ F μΔ=−+−。

取范数，得到本定理。

□
残差的首项随着Arnoldi 算法的进行不断减小。

第二项包含k Δ，并表现出B ，C 对于三次特征值问题的近似总精度的影响。

然而，通过我们的特定公式，第二项具有更高阶
3
10λλ−。

因此通过更新shift 0λ，消除了B 、C 和D 在近似上的影响。

特别地，我们将执行Arnoldi 迭代直到第一项减到比第二项低的水平（当更多的迭代不会导致残差的减少时），并且我们以新的近似特征值作为shift 重新开始。

这导致如下迭代Arnoldi 过程，它在每一步迭代都有更新的shift ，
找到λ满足λση−≤。

算法2.1 CEP 问题的迭代shift-and-invert Arnoldi 过程 1． 输入：,ση，初始x ；最大Arnoldi 步数m ；设定0λσ=； 2． For 1,2,,l =L 直到收敛，do
3． µµµµ200001
(),(32),3,,A L B A B C C A B D A q x x λλλλ==++=+== 4． For 1,2,,j m =L ，do
5． $µµ1
j
q A Bq −=； 6． For 1,2,,i j =L ，do
•$*,i j i h q q =
•$$,i j i q q h q =−；
•µµµµ11
**,,();()i j i j j i j i g q A Cq g q A Cq −−==； •µµµµ11
**,,();()i j i j j i j i
f q A Dq
g q A Dq −−==； 7． End
for 8． $1,j j h q
+=； 9． If 1,0j j h +>，$11,;j j j
q q h ++=；else BREAK ； 10. 计算32j j j I H G F μμμ+++的最大特征对0(,)u μ； 11.
110λλμ−=+；
12. µ311001,1()j j j j h u Aq γλλμ++=−；
13. µµµµ32100()(()())j k j j
j j CQ AQ G DQ AQ F u γλλμ=−−+−； 14. If
120.1γγ≤，BREAK ；
15. End for 16. k x Q u =； 17.
101,λσηλλ−<=；
18. End for
如[19]中所描述的，有如下remarks 。

10Remark 2.1 在第5行和第6行，我们需要µA 的一个分解。

如果这个问题的规模太大以至于得到一个分解不实际，我们需要近似求解这个系统。

有趣的是，这里只需要低精度的解，这是因为残量的阶数将最低是
10λλ−，甚至具有不精确的矩阵向量乘积。

我们的数值试验（见2.4
节）表明，具有一个大小为110−的残量阀值1δ（即µ$µµ1j j
Aq Bq Bq δ−≤）的$q 的求解足以满足收敛。

对于第6行，由于较大误差可以被这里的第二项3
10λλ−抵消，一个残量阀值21
δ=就足够了。

Remark 2.2 1γ和2γ是残量的两个部分，通过存储µ
j AQ ，并且j BQ 、j CQ 和j DQ 在迭代中可用，它们可以在不引起额外矩阵向量乘积的基础上计算出来。

可能也会用它们的界限来替代它们。

Remark 2.3 当进一步降低1γ的更多的Arnoldi 步骤不一定降低残差时，第14行的条件是待定的。

使用的阀值0.1可以放大。

我们也可以使用如下条件：
3210.1()A B C D x γλλλ≤+++，
的确，当µA 被不精确的反置时，（4）不再有效，应该使用这一条件。

Remark 2.4 在这个算法中，当近似特征值落入区域
λση−≤时，我们只更新了shift （17行）。

因此，如果η接近一个非常小的值，这个shift 的更新将会延时，这会降低收敛的速度。

另一方面， 如果η具有一个较大的值，一个更新的shift 可能不会靠近最接近μ的特征值。

然而，这个困难没有在我们的数值试验（见第4节）中出现。

设计一个比第17行在理论上更可靠的标准是很有趣的，但是这显得很难[25]。

2.4 数值试验
在这一节，我们给出一些数值例子来说明这些方法的影响以及它们在实际中的应用。

由于实施起来有很大的CPUtime ，这里我们选取一些不是太大的矩阵作为例子，它们来源于真实的工程问题。

我们希望将来通过同样的方法去研究更大的问题。

这里，所有的计算是在以下环境中完成的：MATLAB 7.1 on PC （Intel （R ） Pentium （R ）4，CPU Processor 3.4GHz 3.39GHz ，
Memory 0.99GB ）。

我们使用算法3.1，固定m=10内部Arnoldi 步骤和外迭代20个循环去计算靠近不同σ值（相对不同的η值）的特征值。

所谓的精确解是通过MATLAB polyeig 计算出来的。

表中的“\”表示在给定参数条件下不收敛。

算例1 结构工程中的动态分析
我们使用结构工程矩阵，它们来源于矩阵市场
[26]
中对于CEPs （三次特征值问题）的
Harwell-Boeing collection 。

这些矩阵全部变现结构工程中的动态分析。

32105,3(1,3,1),,M I M tridiag M bcsstm X M bcsstk X ==−==
矩阵对10(,)M M 的数据分别从（bcsstm01，bcsstk01）、（bcsstm04，bcsstk04）、（bcsstm06，bcsstk06）和（bcsstm08，bcsstk08）中提取。

这些矩阵的大小分别为4848×、132132×、420420×和
10741074×。

citss15中描述的详细情况如下所示：
•Case 1：BCS 结构工程矩阵（特征值矩阵）
来自Harwell-Boeing Collection 的集合BCSSTRUC1中的small test problem 。

bcsstk01：4848×的实对称正矩阵，具有400个非零元素，并且condest=1.5976e+006； bcsstk01：4848×的实对称正半定矩阵，具有24个非零元素，并且condest=inf （该矩阵是奇异的）；
•Case 2：BCS 结构工程矩阵（特征值矩阵）
来自Harwell-Boeing Collection 的集合BCSSTRUC1中的oil rig-not condensed 。

bcsstk04：132132×的实对称正定矩阵，具有3648个非零元素，并且condest=5.6094e +006； bcsstk04：132132×的实对称正半定矩阵，具有66个非零元素，并且condest=inf （该矩阵是奇异的）；
•Case 3：BCS 结构工程矩阵（特征值矩阵）
来自Harwell-Boeing Collection 的集合BCSSTRUC1中的medium test problem-lumped
mass 。

bcsstk06：420420×的实对称正定矩阵，具有7860个非零元素，并且condest=1.2162e +007； bcsstk06：420420×的实对称正半定矩阵，
具有420个非零元素，并且condest=3.4571e +006； •Case 4：BCS 结构工程矩阵（特征值矩阵）
bcsstk06：10741074×的实对称正定矩阵，具有12960个非零元素，并且condest=4.7262e +007；
bcsstk06：10741074×的实对称正半定矩阵，具有12960个非零元素，并且condest=8.2662e +006；
下面的表格给出结果和比较。

12
表 2.1 Case 1的结果和比较
σ
η
计算解 精确解
0.01 -845.4541661789810
0.1 -845.4541661789436 -850 （extreme ）
1.5 -845.4541661783189 -845.454154586488 0.01 -8.43956181917732
0.1 -8.43956180700153 -10 （close 0）
1.5 -8.43956156325571 -8.439561208444
0.01 -114.9561661904075
0.1 -114.9561661342683 -110 （middle ）
1.5 -114.9561661952685
-114.956166445518
表 2.2 Case 2的结果和比较
σ
η CPUtime
计算解 精确解 CPUtime
0.01 1.53 -125.35586039474340.1 1.89 -125.3558603964378-130 （extreme ）
1.5 1.78 -125.3558603731391-125.355860451909 1.40
0.01
1.96
0.06349573884947 -0.58002804445412i 0.1 1.92
0.06349198419845 -0.58002311446271i
0 （close 0）
1.5\ \ 0.06349132173989 -0.58002034582751i
1.64
0.01
1.58 -9.98385917585581 0.1 1.72
-9.89162951742579 -0.00000041223767i
-10 （middle ）
1.5
\ \
-9.983859176103 1.30
表 2.3 Case 3的结果和比较
σ η CPUtime
计算解 精确解 CPUtime
0.01 22.83 -887.2983551993787
0.1 22.87 -887.2983551984910 -900 （extreme ）
1.5 2
2.95 -887.2983551991607
-887.298354460802 96.89 0.01 26.50
1.54141194092581 +3.83872071437039i
0.1 26.54 1.54140653501693 -3.83861091046029i
0 （close 0）
1.5 26.67
\
1.54139177089082
+3.83864678368469i
96.89
0.01 \
\
0.1 23.98 -59.99020405206311 -60 （middle ）
1.5 \
\ -59.990204064020 96.89
14
表 2.4 Case 4的结果和比较
σ
η CPUtime
计算解 精确解 CPUtime
-1250 +2200i
0.1
2333.49
1227.258001278532 +2174.470996486237i 1227.257505828701 +2174.4712080i72931i
1683.50
0.01
607.20
-2456.349471128609+0.000000002723704i 0.1607.99 -2456.327698667066-0.000000170918531i -2500 （extreme ）
1.5608.03
-2456.349260741057-0.000000013174813i -2456.34942855535 1683.50
0.01
612.05
-6.10047457382738 -0.00000000000058i 0.1613.30 -6.10047389513217 -0.00000015348459i 0 （close 0）
1.5611.00
-6.10047433056008 -0.00000000006028i
-6.10024290418462 1683.50
0.01634.36 -1500.1474198940660.1582.78 -1500.147419894066
-1500 （middle ）
1.5
\ \
-1500.14544150022 1683.50
在这个例子中，通过数值试验我们可以做出如下remarks ：
Remark 4.1 所谓的精确解是通过MATLAB polyeig 计算出来的。

同样，我们也可以使用
MATLAB eigs ，但是通常eigs 可能会不收敛。

Remark 4.2 当解收敛时，对于不同的σ值，CPUtime 变化很小。

然而，当选取不合适的σ时，解可能是不收敛的。

Remark 4.3 当σ是复数时，算法的CPUtime 会很长。

Remark 4.4 当矩阵的阶数n 比较小时，CPUtime 几乎是相等的。

然而，当n 比较大时，这个新算法的CPUtime 会小得多。

2.5 结论
我们的结论是：基于标准Arnoldi算法的算法3.1 提供了对于CEP问题的一种快速收敛的方法。

从实用的角度来看，即使当以低精度近似求解shift-and-invert 方程组时，收敛性质也可以保持，这是非常有趣的。

此外，本文的观点可以推广到多项式特征值问题（PEP）。

而且，本文的观点可以对一般的非线性特征值问题（NEP）提供一些见解。

许多问题有待研究。

例如，对于所有类型的大型矩阵（特别当A的条件数很大时）的有效方法有哪些？如何选取这个任选参数η？这些需要进一步的研究。

近似shift-and-invert的详细研究和分析是将来研究的主题。

16第三章 求解多项式特征值问题的Jacobi-Davidson 方法
3.1预备知识
3.1.1 Jacobi-Davidson 方法简介
这一章我们研究另一种子空间投影方法——Jacobi-Davidson 方法
[20,22,27]。

在这种方法中，
子空间k κ由一种Jacobi-Davidson 方法构建，该法利用校正方程的解扩充子空间，使Ritz 对很快收敛于特征对。

求解多项式特征值问题的方法有线性化方法和直接投影方法等，本章将以二次特征值问题（QEP ）为代表，研究QEP 问题的线性化Jacobi-Davidson 方法和直接Jacobi-Davidson 方法。

3.1.2线性化简介
为了数值求解次数大于1的PEPs 问题，我们将这些问题转化为更大维数（mn ）的GEPs 问题。

这与将一个高阶常微分方程（ODE ）转化为一个1阶的ODE 的想法相同。

这些转化被认为是线性化，并且对于一个给定的PEP 问题，它们不是唯一的。

定义 3.1
[28]
我们说这个对（A,B ）是次数为m 的n n ×阶矩阵多项式(,)P A λ的一个线性化，
如果存在维数为mn 且行列式不为零的两个矩阵()E λ和()F λ，使得
(1)(,)
()()()0n m P A E A B F I λλλλ−⎛⎞
−=⎜
⎟⎝
⎠
这样就将一个多项式特征值问题（PEP ）转化为一个广义特征值问题（GEP ），然后可以
借助GEP 问题的数值方法求解相应的PEP 问题。

常用的多项式线性化方法有相伴线性化和对称线性化等。

假定
000000
n n A I A I ⎛⎞⎜⎟⎜⎟⎜
⎟=⎜⎟⎜⎟⎜⎟⎝
⎠L O M M O L
L
， 12
00000m n n
A A A I
B I −−−⎛⎞
⎜⎟⎜⎟⎜⎟=⎜⎟⎜⎟⎜⎟⎝
⎠
L L O M L
，
则矩阵束A B λ−是(,)P A λ的一个线性化，称为相伴线性化。

另外，当
0011
200
2
11
000m m m A A A A A A A A A A A −−−⎛⎞⎜⎟⎜⎟
⎜⎟=⎜⎟⎜⎟⎜⎟⎝⎠L L M M
M L
， 0
0120100000000
m m A A A B A A A A −⎛
⎞
⎜⎟
⎜
⎟⎜⎟=⎜⎟
⎜
⎟⎜⎟−⎝
⎠
L M M M M L L
L
时，如果m A 是非奇异的，则A B λ−是一个对称线性化
[28]。

在这两种线性化之中，相伴线性化是最常用的。

假定(,)P A λ是一个矩阵多项式，并且A B λ−是它的一个线性化。

λ作为PEP 问题的一个特征值的条件数应该比λ作为通过线性化得到的GEP 问题的一个特征值的条件数小。

对这一结果没有证明。

启发式的论点是对于PEP 问题的这类可能的扰动比对于相应的GEP 问题的这类可能的扰动小。

一个有趣的问题是如何找到“最佳”线性化，即一个λ作为GEP 问题的一个特征值的条件数最小的线性化。

然而，将PEP 问题线性化后，存储量和计算量明显增加，并且原始PEP 问题中矩阵的很多良好性质比如对称性、正定性和稀疏性等可能不会在GEP 问题中得到体现，这是线性化的一个不足之处，因此寻找一种直接将PEP 问题投影到子空间求解的数值方法又是非常必要的。

下面几节我们将以二次特征值问题为代表，给出其线性化Jacobi-Davidson 方法和直接Jacobi-Davidson 方法的算法，并给出数值算例比较结果。

3.2 线性化Jacobi-Davidson 方法
二次特征值问题（QEP ）的求解主要有一些线性化方法和投影方法[1,4,7,21,27]。

这一节我们研
究QEP 问题的线性化Jacobi-Davidson 方法。

QEP 问题的两类常用的线性化是
18000N N K C M λ⎛⎞⎛⎞−⎜⎟⎜⎟−−⎝⎠⎝⎠ 和 000
K
C
M
N N
λ−⎛⎞⎛⎞−⎜⎟⎜⎟⎝⎠⎝⎠。

其中N 可以是任意非奇异n×n 阶矩阵。

这两种线性化之间的选择可能要依赖于M 和K 的非奇异性。

一般地，N 选取为单位矩阵或单位矩阵的倍数，例如M I 或K I
[1]。

关于Jacobi-Davidson 方
法，通常方式是利用线性化将QEP 问题转化为等价的GEP 问题（2）或（3），然后对GEP 问题采取Jacobi-Davidson 方法。

我们将这种方法简称为“GJD ”。

Sleijpen 和van der Vorst
[29,30]
介绍的对于线性特征值问题
Ax x λ= （3.1）
的Jacobi-Davidson 方法是一个迭代投影方法。

如果V是一个具有正交列的矩阵，且(,)k k v σ是投影问题T
V AVv v σ=的一个特征对，那么对应的逼近（3.1）的一个特征对的Ritz 对(,)k k u σ，
k k u Vv =由如下过程改善。

矩阵V 通过k u 的一个正交校正t 扩充，并且V 由[V , t]代替。

最理想的正交校正t
[31]
是求解方程
()()k k A u t u t λ+=+，k t u ⊥ . （3.2） 由于k t u ⊥，算子A 可被限制到正交于k u 的子空间，产生了()()T
T
k k k k I u u A I u u −−，并且（3.2）可以重写为
()()()()T T
k k k k k k I u u A I I u u t A I u λσ−−−=−−
这里我们假设k u 标准化1k u =。

通过k σ近似未知λ，最终我们得出对于更新k t u ⊥的Jacobi-Davidson 校正方程：
()()()T T
k k k k k k I u u A I I u u t r σ−−−=− （3.3） 其中:()k k k r A I u σ=−表示(,)k k u σ的残量。

对于不同类型特征值问题的Jacobi-Davidson 方法的执行细节可以在[31]中找到。

特别地，数值试验表明
[32]
校正方程只需要近似求解出。

通常，一个预处理Krylov 子空间方法的仅仅少数几
步足以得到一个子空间V 的好的扩充t 。

由上述分析，我们很容易将Jacobi-Davidson 方法推广到广义特征值问题
Ax Bx λ=，
校正方程变为
()12H H k k k k k k k H H k k k k k t Bu u u u I A B I r t u Bu u u θ⎛⎞⎛⎞⎛⎞−−−=−⎜⎟⎜⎟⎜⎟⎝⎠⎝⎠⎝
⎠
其中:()k k k r A I u θ=−表示(,)k k u θ的残量。

我们记：i V 代表n 阶方阵的投影子空间的基组成的矩阵，MGS 表示修正的Gram-Schmidt 正交化过程。

通过线性化QEP 问题2
()0M C K x λ
λ++=，我们得到GEP 问题()0A B u λ−=，其形
式为（2）或（3）。

然后我们将GEP 问题投影到以矩阵2i V （其形式见下面的算法）的列为基的子空间上，从而得到线性化Jacobi-Davidson （GJD ）算法：
算法 3.1
1: 选取 1n
v ∈¡
，11[]V v =，精度ε, 正整数 m;
2: 对 1,2,,i m =L 执行以下步骤 2.1: 20;0i i i V V V ⎛⎞
=⎜
⎟⎝⎠
2.2: 计算 22,;Ai i Bi i W AV W BV == 2.3: 计算 22,;H
H
Ai i Ai Bi i Bi H V W H V W ==
2.4: 计算Ai i i Bi i H z H z θ=的目标特征对 (,)i i z θ； 2.5: 计算Ritz 向量2i i i u V z =； 2.6: 计算残量()i Ai i Bi i r W W z θ=−； 2.7: 收敛性检查，如果2
i
r ε<，则停止;
2.8: 如果不收敛，则求解校正方程
()12H H i i i i i i i H H i i i i i t Bu u u u I A B I r t u Bu u u θ⎛⎞⎛⎞⎛⎞−−−=−⎜⎟⎜⎟⎜⎟⎝⎠⎝⎠⎝
⎠；
2.9: 如果dim()1i V m ≤−， 则11(,)i i i V MGS V t +=，否则转 3; 3: 1[]m V u =, 转2，重新开始。

算法3.1将n 阶的QEP 问题转变成了2n 阶的GEP 问题，增加了存储量。

3.3直接Jacobi-Davidson 方法
为了扩展Jacobi-Davidson 方法到非线性特征值问题
()0T x λ=，
20我们使用校正方程
()()()T T k k k k
k k T T k k k k
p u u u I T I t r u p u u σ−−=−，k t u ⊥， （3.4）
其中'
:()k k k p T u σ=且:()k k k r T u σ=，(,)k k u σ是()0T x λ=的一个特征对的当前近似，由投影问题获得
()0T V T Vv λ=，u Vv =.
显然方程（3.4）等价于
()k k k T t p r σα−=−，
其中α需选择使得k t u ⊥。

求解t 我们得到
11'
()()()k k k k k k k t u T p u T T u ασασσ−−=−+=−+， （3.5）
并且由于k u V ∈，因此{,}{,}span V t span V t =%，其中
1'()()k k k t T T u ασσ−=%.
因此，由V 的列张成的空间在逆迭代得到的方向上进行扩展，如果k σ被选为k u 的Rayleigh 函数，该逆迭代会有三次收敛。

因而，如果方程（3.4）被精确地求解，可以得到三次收敛。

数值试验表明，即使是（3.4）的一个具有一些预处理GMRES 步骤的一般的近似求解也会导致快速收敛，它几乎是三次的。

现在我们采用将QEP 问题直接投影到子空间中的Jacobi-Davidson 方法，然后得到QEP 问
题的直接Jacobi-Davidson 方法的算法如下。

其中%(,)i i x θ是该QEP 问题的一个Ritz 对，残差
%2()i i i i r M C K x θθ=++%。

算法 3.2
1: 选取1n
v ∈¡
，11[]V v =，精度ε, 正整数 m;
2: 对 1,2,,i m =L 执行以下步骤
2.1: 计算,,;Ki i Ci i Mi i W KV W CV W MV === 2.2: 计算,,;H
H
H
i i Ki i i Ci i i Mi K V W C V W M V W === 2.3: 计算2
()0i i i i i i K C M s θθ++=的目标特征对(,)i i s θ；
2.4: 计算Ritz 向量%;i i i x V s =。