空间迭代法

也有人从其它方面改进子空间迭代法，例如
¾ 改变投影矩阵[12]，是用求解投影标准特征值问题来替代求解投影广义特征值问题，从而避免求解
投影刚度矩阵。
从文献上看，大多数改进方案主要是针对矩阵阶数比较小的广义特征值问题，并且只和传统的标
准子空间迭代法相比较，没有各个改进方案之间的相互比较。
用本文较大规模的算例对自适应多重逆迭代、超松弛幂迭代、大移频和改变投影矩阵进行比较后
理想，而移频能很好的解决高阶重频特征向量的收敛速度慢的情况。因此，可以把超松弛幂迭代改进 方案用于移频的子空间迭代法。
subspace iteration)[7]的步骤：
1)选取初始向量矩阵 X (0)
2)计算 P(k +1) = MX (k ) , k =0,1,2,…
(k +1)
3)解出 X
(k +1)
KX
= P(k+1)
(k +1)
(k)
4) 计算 P = M X
5) 解出 X (k +1)
( k +1)
K X (k+1) = P
子空间迭代法的一种新改进方案
宫玉才，陈璞
（湍流与复杂系统国家重点实验室，北京大学工学院力学与空天技术系 北京 100871）
摘 要：对于结构固有振动所对应的大型广义特征值问题的求解，子空间迭代法是最可靠解法之一，并 在通用有限元软件中得到了广泛应用。为了进一步提高子空间迭代法的效率，近十年来在又发展了自适 应多重逆迭代、超松弛幂迭代、大移频和改变投影矩阵法等一些新的改进方案。本文总结了子空间迭代 法近十年来的进展，把大移频和超松弛幂迭代两种改进方案结合在一起，提出了一种新的改进方案。该 改进方案很好地把超松弛、矩阵幂迭代和移频三种思想结合在一起，使子空间迭代法的计算效率有了进 一步提高。较大规模的例题表明，新的改进方案效率把大移频的效率平均提高了 20％～30％。 关键字：子空间迭代 广义特征值问题 固有振动 有限元分析
％。
文章分为五部分，第一部分为引言，简单介绍了近十年来子空间迭代法的进展；第二章介绍了
超松弛幂迭代和大移频；第三章把移频与超松弛幂迭代结合在一起，并从子空间维数选取等方面进行
源自文库
了研究；第四部分为算例和讨论；最后一部分为结论。
2 子空间迭代算法
子空间迭代法是反幂法的推广。反幂法的矩阵迭代每次只计算一个特征值和特征向量，从最低 阶特征值开始，逐渐推进到高阶特征值。后来开始用多个向量同时进行迭代(simultaneous iteration)， 并引入了 Gram-Schmidt 过程来正交化向量，能同时求出几个特征值和特征向量。上世纪七十年代初， 在其中加入了子空间上的 Rayleigh-Ritz 过程[1]，明显改善了收敛速度。
X = X Q (k +1)
(k +1) (k +1)
9)计算特征值的相对误差
tol
(k i
+1)
，并判断是否全部收敛。若全部收敛，则退出循环；否则，跳
到 2)继续进行迭代。 矩阵幂迭代加速收敛的效果明显，但稳定性有些差。在有些情况下，在一次循环内进行两次逆
迭代之后， X (k +1) 会变为几乎线形相关的，使投影后的刚度矩阵和质量矩阵变成奇异的，导致数值不
当 K 和 M 规模比较大，例如矩阵阶数超过 1000 的大型特征值问题，子空间迭代法无疑是最受
青睐的方法之一。许多有限元软件, 例如 ABAQUS、ADINA、ANSYS 和 NASTRAN，早已把子空间 迭代法作为它们的广义特征值问题求解器。与迭代 Lanczos 法和迭代 Ritz 向量法相比，子空间迭代法 的速度慢一些，但稳定性要好很多，并且算法易于实现。
发现，超松弛幂迭代和大移频是表现最好的两种加速方案。
与标准字空间迭代相比较，超松弛幂迭代的效率平均提高了 50%，大移频的效率则高于超松弛幂迭代
的效率，并随着所求模态个数增加而提高。
本文把超松弛幂迭代和大移频结合在一起，提出了一种新的改进方案，并从子空间维数选取、
每次三角分解的最大迭代次数等方面进行了研究。新的改进方案把大移频的效率平均提高了 20%～30
其中 λ (k +1) =
λ (k ) q
1 − δ (k +1)
82
1
δ (k +1) 的经验取值为 δ (k +1)
= 0.01δq(k )
=
0.01
⎡ ⎢
(λq(k
)
)
2
(k
xq
+1) T
(k +1)
pq
−
2λq(k )
(k +1)T
xq
p(k +1) q
⎤2 + 1⎥
⎢⎣
⎥⎦
具体的证明过程可参考[9]。
μ = λs+1 + λs+2
(5)
2
这时需要放松 Sturm 序列校核，即需要检查已收敛的特征值个数加 1 与因子矩阵 D 中负对角元的个 数是否相符。大移频移位方案比传统的移位方案平均快 17%[11]。
3 一种新的改进方案
超松弛幂迭代[9]改进方案的加速效果很有成效，缺点是在 λp 附近的高阶重频特征向量收敛速度不
1. 将 X (k ) 进行 M-正交归一化
2. 解试向量矩阵 (LDLT ) X (k+1) = MX (k ) ；
3. 计算 K 和 M 在 X (k+1) 上的投影
K (k +1)
=
T (k +1)
XK
X M (k +1) ， (k +1)
=
T (k +1)
XM
(k +1)
X
83
4.
求解 q 阶广义特征值问题
¾ 矩 阵 幂 迭 代 (matrix power accelerated subspace iteration)[7] ， 是 用 X (k+1) = (K −1M )2 X (k ) 代 替
X (k+1) = K −1MX (k ) ，也就是每进行两次逆迭代，才有一次 Rayleigh-Ritz 过程。
6)计算 K 和 M 在 X (k ) 上的投影
K (k +1)
=
T (k +1)
XK
X M (k +1) ， (k +1)
=
T (k +1)
XM
(k +1)
X
7)求解投影后的 q 阶广义特征值问题
K Q = M Q Λ (k +1) (k +1)
(k +1) (k +1) (k +1)
8)形成新的近似特征向量
子空间迭代法假设 q 个线形无关的初始向量同时进行迭代，求得前 p 个特征值和特征向量。传 统上，q = min(2p, p+8)，用两次迭代之间的特征值的相对误差
( k +1)
tol i
=
λ λ − (k +1)
(k)
i
i
λ (k +1) i
< tol , i =1, 2, …, p
(2)
来控制收敛，tol 通常取 10-6。
的增多抵消了一部分用移频所减少的工作量。
移频子空间迭代法中，子空间维数可取 q = max( s , 4) ，其中 s 为 L 中一行的平均非零元个数，
Imax 由第 II 与 III 步计算量之比确定[19]
I max
=
⎛ max ⎜
⎝
0.5Ns2 2Nsq + 3Nq2 + 2Nq
+ 10q3
⎞ ,4⎟
⎠
以下是移频子空间迭代算法的主要步骤： I．初始化 1. 确定子空间的维数 q
2. 选取初始向量矩阵 X 0
3. 设定每次移轴的最大迭代次数 Imax II．移轴与 Sturm 序列校核 1. 计算移轴 μ，应设法保证它不是特征值
2. 分解移轴刚度矩阵 K − μM＝LDLT
3. Sturm 序列校核 III．迭代 Imax 次，完成后转向 I
的比值越来越小，收敛速度越来越慢。特征值移动 μ 之后，收敛速度为 (λq+1 − μ) /(λi − μ) 。如果 μ
的取值合适，可以使 (λq+1 − μ) /(λi − μ) λq+1 / λi ，从而加速迭代过程的收敛。
每次特征值的移位，都会使矩阵 K − μM 发生变化，需要重新对其进行三角分解。三角分解次数
广义特征值问题求解方法的选择取决于系统自由度大小、矩阵的稀疏性、所求解特征值的个数和 其在特征值谱上的位置。在结构动力学中需要求解
Kϕ − λMϕ = 0
(1)
的 p 个低阶特征值与特征向量( λi , ϕi )，i=1,2,…, p。其中，刚度矩阵 K 为稀疏的对称正定对称
（Symmetrical Positive Definite, SPD）矩阵，质量矩阵 M 一般为半正定的对角矩阵。但矩阵束正定， 即对任意的正数 μ，矩阵 K+μM 正定。在实际工程计算中，矩阵 K 和 M 的阶数可以从几十到几百万。
稳定。
(k +1)
为了避免这种情况，可以用 X 和 X (k +1) 的线形组合
(k +1)
(k +1)
X
=X
α1
+
X
α (k +1) 2
(4)
来代替矩阵幂迭代中的 X (k +1) 。α1 和α2 为系数矩阵，其分量的取值要使收敛速度尽量快。
(k +1)
(k +1)
当X
=X
− λ (k+1) X (k +1) 时，收敛速度比应用矩阵幂迭代的子空间迭代法要快。
(k +1)
(k +1)
(k +1)
的线形组合 X
=X
λ − X (k+1) (k +1) 来 计 算 试 向 量 。 其 中 ， X
= K −1MX (k ) ，
(k +1)
X (k +1) = K −1M X
= (K −1M )2 X (k ) 。
¾ 大移频(large shifting)[10]，是把 μ 放到下两个待收敛特征值之间。
1 引言
在结构动力学、电磁学、声学、量子化学等计算领域中经常会求解广义特征值问题。广义特征值 问题作为科学计算和工程应用中的一个重要研究课题，已经有许多比较成熟的算法。在现今应用中经 常需要求解几十到几百个的特征值，矩阵规模也不断增加，因此需要寻找有效、稳定和快速的算法来 提高单机上的解题规模和速度。
2.2 大移频
移频，又称谱变换(Spectrum Transformation)或移位，相当于在频率谱上对频率作了一维坐标变换， 从而加速待求的特征值和特征向量的收敛速度。移频[2]最早由 Bathe 在上世纪 70 年代提出，近十年来 随着求解特征模态个数的增加才广泛应用到特征值问题解法中。
在标准的子空间迭代法中，第 i 个特征值的收敛速度可以用 λq+1 / λi 来衡量。随着 i 的增加，λq+1 / λi
子 空 间 迭 代 法 从 提 出 [1] 到 至 今 ， 人 们 一 直 都 在 研 究 和 改 进 它 ， 发 展 了 超 松 弛 因 子 (overrelaxation)[2]、移频(shifting)[2]、Chebyshev 多项式(Chebyshev polynomials)加速[5]、选择性二次逆迭代 (selective repeated inverse iteration)[6]等很多加速方案。近十年来，在这些加速方案的基础上又发展了 许多改进方案，主要有
K Q = M Q Λ (k +1) (k +1)
(k +1) (k +1) (k +1)
5. 形成新的近似特征向量 X (k +1) = X Q (k +1) (k +1)
6. 按模态误差判断特征值和特征向量的收敛，移出已收敛的特征向量,并在 X 中加入随机向量 传统的移位方案[2,20]是相当保守的，它只用来加速数轴上其右侧模态的收敛速度，因为它左侧的 模态已经收敛。大移频(large shifting)[10]是将移位移到下两个待收敛特征值之间，即
ε |
Kϕ
(k i
)
−
λi(
k
)
M
ϕ
(k i
)
|2
<
|
Kϕ
(k i
)
|2
ϕ
(3)
来控制收敛。 εϕ 可取 10-3 或 10-4，对于矩阵阶数比较大且求解模态数比较多时，最好取 10-4。
81
2.1 超松弛幂迭代 为了方便介绍超松弛幂迭代[9]改进方案，下面先给出矩阵幂迭代(matrix power accelerated
¾ 自适应多重逆迭代(adaptive multiple inverse iteration)[8]，是在迭代前对每个试向量 xi 进行对应的
多次逆迭代，从而加速高阶特征向量的收敛。
80
(k +1)
¾ 超松弛幂迭代(over-relaxed and matrix power accelerated subspace iteration)[9]，是用 X 和 X (k +1)
子空间迭代法自从上世纪 70 年代提出就得到广泛的应用，并在相当长的时间内成为有限元分析
的首选特征值问题解法。当时，需要求解的特征值个数较少(少于 20 个) ，但在现今应用中，抗震规
范要求 90％的质量参与系数经常导致求解几十甚至几百个振型。
当求解的模态数较多(大于 40)时，用传统上两次迭代之间的特征值的相对误差来控制收敛不能很 好控制高阶特征向量的收敛，需要用模态误差[10]