稀疏特征问题求解 赵永华 中国科学院计算机网络信息中心 超级计算中心 yhzhao@sccas.cn

提要 稀疏特征值问题背景 稀疏矩阵特征值问题软件包列表 稀疏矩阵特征问题求解方法 单/多向量（single and multiple vector）迭代法 Lanczos 和Arnoldi方法 Jacobi-Davidson方法 采用预处理技术的预条件CG（Conjugate-Gradient）法 稀疏矩阵的存储格式 并行软件包设计问题

特征值问题背景 许多大规模科学应用和工程计算中，（包括材料科学、振动力学、电磁学、量子化学等），一个重要的任务是计算特征值问题： Ax=λMx 的部分特征对。通常情况下A和M是大规模、稀疏对称正定矩阵。 对大规模和稀疏特征问题至今没用最好的求解方法。对这些复杂的问题，需要在矩阵的数学特性、算法的求解效率和精度要求等中做出权衡。

电磁场计算 一切电磁波在宏观媒体中都服从Maxwell方程。这些方程通过有限元离散后得到一个广义特征值问题 随着问题规模的增大、有限元基函数解的提高以及有限元网格的加细，离散后的特征值问题的阶数可高达千万阶到上亿阶。 电子结构计算 电子结构中的许多模型对应的方程经过有限元离散化后得到下列形式 其中 是电子密度(特征向量)， 是特征值组成的对角矩阵。

Nano/Material-Science Accelerator Modeling

稀疏特征问题求解方法: Lanczos和Arnoldi方法 单/多向量（single and multiple vector）迭代法 幂(Power)迭代方法、逆迭代方法、Rayleigh商迭代(RQI)方法 Lanczos和Arnoldi方法 Lanczos和Arnoldi过程主要构造Krylov子空间的一个标准正交基： 然后在这个子空间上计算近似的Ritz特征对。该过程是一个三项递推求解过程,在每步仅需要一个矩阵-向量积。 基于不同的正交化方法：完全重正交、选择性重正交和局部重正交等Lanczos和Arnoldi方法； 基于不同的重启动策略：显示重启动和隐式重启动Lanczos和Arnoldi方法； 位移逆转换或谱转换：使用基本的Lanczos和Arnoldi方法总是收敛于最大特征值，对于中间特征值需要使用位移-逆转换。另外，为了计算广义特征问题，需要将 转化为： 为了避免迭代收敛到已求得的特征对，这时需要采用紧缩处理。

稀疏特征问题求解方法:Jacobi-Davidson方法 基于块策略：块Lanczos和Arnoldi算法和band Lanczos。 Jacobi-Davidson方法 Jacobi-Davidson方法的思想来源于JOCC(Jacobi orthogonal Component Correction)和Davidson两种方法。由G. L.G.Sleijpen等于1996年提出。 对于Lanczos和Arnoldi方法而言，当求解问题规模增加时，出现的特征值簇将严重影响它们的性能和健壮性。另外，由于矩阵A的逆不能被直接求解，这样预处理成为弥补迭代法的必然选择。这使得Davidson和Jacobi-Davidson方法成为了Arnoldi方法的扩充。 同Lanczos和Arnoldi方法类似，Jacobi-Davidson也是从一个初始向量开始，连续构造由一组正交基构成的子空间，并由此得到所求解的近似特征对，而这些特征对的求解通常可通过Rayleigh商方法得到。但Jacobi-Davidson 构造的标准正交基不同于Krylov子空间的基 。

稀疏特征问题求解方法:Jacobi-Davidson方法 算法思想 假设当前正交基矩阵V=[u1,u2,…,uk],得到一个投影矩阵 1.计算:Bs=θs 得到特征对（θj, sj）,由此得到A的Ritz向量uj=Vsj 2.寻找一个修正向量 ，使得: 3.从校正方程 中求解v。 4. 对v做正交化并扩充当前正交基，再计算新的Ritz对。 如此循环，这个算法就是基本的Jacobi-Davidson方法。

稀疏特征问题求解方法:Jacobi-Davidson方法 预处理近似计算校正方程 校正方程可以使用线性系统的近似迭代法（像GMRES或BCGSTAB）求解。如果结合适当的预条件子 可以加速迭代收敛过程。记 如果采用左预处理，在迭代中需要对某个向量 计算 上述运算可通过两步完成： 首先计算: 然后使用左预处理,从 计算z

稀疏特征问题求解方法:Jacobi-Davidson方法 重启动 JD方法在计算过程中需要保留标准正交基，当迭代步比较多时对存储要求太高，因此需要重启动。 一个相应的重启动策略为:Dynamic thick 重启动策略。 紧缩 为了避免迭代收敛到已求得的特征对，这时需要采用紧缩处理。假如已有Ritz向量u1,u2,…,uk ，带紧缩的校正方程为: 其中Q=（ u1,u2,…,uk ）

稀疏特征问题求解方法:Jacobi-Davidson方法 块Jacobi-Davidison实现 JD的块算法采用具有k个初始向量的块作为算法的开始，每次扩展的基也是一个k个向量的块。这些向量是k个校正方程的近似解，而每个用于不同的Ritz对求解。 JDCG方法 Y.Notay2002年针对实对称矩阵的最小特征对，结合JD方法和预条件共轭梯度（CG）法，提出JDCG方法。 JDQR方法和JDQZ方法 并行化 JD方法使用数据并行方法比较容易并行实现。每个处理器可以保留每个基向量和特征向量的子集，所有的向量修正可在本地进行。同串行算法的不同在于内积运算需要一个全局归约操作。而且矩阵-向量积和预处理操作必须被提供。 软件包：PRIMME 求解问题：电磁场计算；电子线路稳定性

稀疏特征问题求解方法:预条件CG方法 预条件CG方法 CG方法是Chebyshev方法的改进，基于这一基本方法产生了几个不同的CG求解方法，主要包括:LOBPCG、DACG方法。 LOBPCG(Locally Optimal Block ，Preconditioned Conjugate Gradient) 该方法是一种快速求解对称特征问题的方法（A.VKayazev,2001）。由于具有简单性、健壮性和快速收敛性，该方法作为Lanczos和Davidson方法的潜在 竞争方法，引起人们的关注。 该方法是一种梯度法，在选取合适的预条件下，对于求解大规模稀疏问题比较有效，而且它可以块方式求解多个特征值。 LOBPCG也是一种块投影算法，它的子空间由上一步向量、当前向量和局部梯度方向组成，然后在子空间上应用Rayleigh-Ritz方法。 LOBPCG方法 典型的软件包是BLOPEX ,MATLAB

稀疏特征问题求解方法:预条件CG方法 是一种采用最优化方法的CG方法。通过Rayleigh商的CG最小化，该方法可以计算部分特征对。 DACG(Deflation-Accelerated Conjugate Gradient ) 是一种采用最优化方法的CG方法。通过Rayleigh商的CG最小化，该方法可以计算部分特征对。 当使用有效的预处理时，DACG的效率要好于ARPACK软件包。并且当计算少量特征对时， DACG在数值和效率方面比得上Jacobi-Davidson方法。 比较有效的预条件包括块Jacobi、近似逆预条件子:FSAI和AINV,由此开发出Jacobi-DACG，FSAI-DACG和AINA-DACG方法。 这些方法当前用于有限元、混合有限元和有限差分的特征问题中。 相关的软件包PDACG

稀疏矩阵存储格式 绝大多数迭代求解的效率取决于矩阵-向量积，而矩阵-向量积取决于矩阵所使用的存储策略。存储策略的选取需要考虑具体的应用问题。 几种常用的存储格式： 压缩行(列)存储策略(RS:Compress Row storage) 块压缩行(列)存储策略(BCRS:Block Compress Row storage) 压缩对角存储格式(CDS:Compressed Diagonal Storage ) 锯齿对角存储格式(JDS:Jagged Diagonal Sstorage)

稀疏矩阵存储格式 创建三个向量： val:向量以行方式存储矩阵的非零元素值； col_ind:存储val中元素的列标； row_ptr:存储行中第一个元素在val中位置。 Val=[10,-2,3,9,3,7,8,7,3,8,7,5,8,9,9,13,4,2,-1] Col_ind=[1,5,1,2,6,2,3,4,1,3,4,5,2,4,5,6,2,5,6] Row_ptr=[1,3,6,9,13,17,20]

稀疏矩阵存储格式 稀疏块压缩行存储BCRS(Block compress row storage) 一个更有效的存储方法是稀疏块压缩行存储格式。 这样的块矩阵典型的出现在偏微分方程的离散。 与CRS存储格式类似，BCRS存储格式可通过四个数组表示： val=[11,12,0,22,31,32,41,42,33,0,43,44,55,56,0,66,0,0,67,0,77,78,87,88] AN=[1,5,9,13,17,21] AJ=[1,1,2,3,4,4] AI=[1,2,4,6,7]

稀疏矩阵存储格式 压缩对角存储(the compressed diagonal storage:CDS ) 压缩对角存储是另一种稀疏矩阵存储格式，该存储方式利用了特殊矩阵的结构形式。这种结构与有限元或有限差分分解有关。 使用数组Val(6,-1:1)存储右边矩阵。

稀疏矩阵存储格式 锯齿状存储

并行软件包设计问题 用户数据接口 并行软件包的一个重要问题是应用程序同软件包的接口，特别是由问题产生的矩阵如何被描述。基于这一考虑，软件包提供具有各种灵活度的不同机制。 逆通信是解决软件包独立于矩阵描述的一种很好解决办法。 将矩阵输出到文件是一种接口方式，但灵活性最差。 并行方法 特征问题并行软件通常采用基于数据并行的SPMD模式，支持MPI并行。并要求用户提供某些操作的MPI-并行实现，例如矩阵-向量积。 另一种并行方法采用了OpenMP共享内存模式。 极少的软件包采用了JAVA的线程级并行。

并行软件包设计问题 程序语言 C，C++和FORTRAN77/90 支持的标量类型和精度 实数和复数标量类型，单精度和双精度。 健壮性 对并行软件包而言，算法的健壮性是一个重要问题。为此需要使用一些通用的库函数和稳定的算法，例如LAPACK,BLAS, PETSC，MATLAB等。

谢谢！