大型稀疏矩阵的LU分解及特征值求解
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
• 直接法基于高斯消元法，即计算A的LU分解。A通常要经过一系列置换排 序，可归并为左置换矩阵P，右置换矩阵Q。基本步骤如下： 1）符号分析： 得到置换排序矩阵 P Q
2）数值分解：
3）前代 回代：
I.S.Duff, A.M.Erisman, and J.K.Reid. Direct Methods for Sparse Matrices. London:Oxford Univ. Press,1986. J.J.Dongarra,I.S.Duff, ... Numerical Linear Algebra for High-Performance Computers. G.H.Golob, C.F.Van loan. Matrix Computations. The Johns Hopkins University Press. 1996
为何需要密集子块(Dense Matrix)
http://eigen.tuxfamily.org/index.php?title=Benchmark
多波前法(multifrontal)简介
• 发展
Duff and Raid [2] J.W.H.Liu等分析，改进 [3] T.A.Davis开发UMFPACK [4]
2 密集子阵
根据符号分析，数值分解算法的不同，直接法有以下几类[15]： 1）general technique(通用方法)： 主要采用消去树等结构进行LU分解。适用于非常稀疏的非结构化矩阵。 2）frontal scheme(波前法)
LU分解过程中，将连续多行合并为一个密集子块(波前)，对 这个子块采用BLAS等高效数学库进行分解。 3）multifrontal method(多波前法)
f,g,h gg·
h ··
{h,i,j}
消去树
bej e ·· j ··
e,i,j i ·· j ··
h,i, j ii· j ·j
d,i,j i ·· j ··
消去树
消去树是符号分析的关键结构，其每个 节点对应于矩阵的一列（行），该节点只 与其父节点相连，父节点定义如下：
将符号结构相同的多行(列) 合并为多个密集子块，以这些密 集子块为单位进行LU分解。这些子块的生成，消去，装配，释放等都需 要特定的数据结构及算法。
稀疏排序
排序算法的作用是减少矩阵LU分解过程中产生的注入元，求解 矩阵的最优排序方法是个NP完全问题（不能够在合理的时间内进行 求解），对具体矩阵而言，目前也没有方法或指标来判定哪种算法 好。因此实测不同的算法，对比产生的注入元，是确定排序算法的 可靠依据。
GSS1.2
15秒
2012
3.0G 4核
GSS 2.3
4秒
2014
i7 8核
GSS 2.4
1秒
• 硬件的发展 CPU，内存等
• 稀疏算法逐渐成熟
稀疏排序，密集子阵 multifrontal ,supernodal…
• 数学库
BLAS，LAPACK
www.netlib.org
稀疏LU分解算法的关键
稀疏LU分解 （不考虑零元的）LU分解 1 零元是动态的概念，需要稀疏排序,减少注入元(fill-in)
LU分解形成frontal
10阶矩阵。 蓝点代表非零元。红点表示分解产生的注入元(fill-in) Frontal划分{a}, {b}{c}{d} {e} {f,g}{h,i,j}
Frontal的装配，消去，更新过程
{a}
ach c ·· h ··
{c} {b}
c,g,h g ··
h ··
{f,g} {e} {d}
主要的排序算法有最小度排序（MMD，minimum degree ordering algorithm）和嵌套排序（nested dissection）两种。矩阵排序方面相应 的成熟软件库有AMD[12] 、COLAMD、METIS[13]等。
Yannokakis M. Computing the minimum fill-in in NP-Complete. SIAM J.Algebraic Discrete Methods, 1981, 2:77~79
求解器的飞速发展
• BBMAT
http://www.cise.ufl.edu/research/sparse/matrices/Simon/bbmat.html
38744阶，分解后元素超过四千万.
1988 巨型机cray-2上
>1000秒
2003
4G
umfpack4
32.6秒[4]
2006
3.0G
稀疏矩阵复杂、多变
• 基本参数
对称性，稀疏性，非零元分布
• 敏感性，病态矩阵
条件数
• 格式多变
Harwell-Boeing Exchange Format 。。。
• 测试集
Harwell-Boeing Sparse Matrix Collection UF sparse matrix collection
近似最小度排序算法 – 商图
近似最小度排序（AMD，Approximate Minimum Degree OrderingAlgorithm）算法于 1996年左右由Patrick R. Amestoy等学者提出
AMESTOY, P. R., DAVIS, ... 1996a. An approximate minimum degree ordering algorithm. SIAM J. Matrix Anal. Applic. 17, 4, 886–905.
• 基本算法
利用稀疏矩阵的特性，得到一系列密集子阵（波前）。将LU分解 转化为对这些波前的装配，消去，更新等操作。
• 多波前法的优点
波前是dense matrix ,可直接调用高性能库（BLAS等） 密集子阵可以节省下标存储 提高并行性
• 目前主要的求解器
UMFPACK,,GSS,MUMPS,PARDISO, WSMP,HSL MA41等
大型稀疏矩阵 的LU分解及特征值求解
陈英时 2016 . 1. 9
2013. 7. 20
稀疏矩阵求解的广泛应用
• 矩阵求解是数值计算的核心[1] • 稀疏矩阵求解是数值计ຫໍສະໝຸດ Baidu的关键之一
偏微分方程，积分方程，特征值，优化… 万阶以上dense matrix不可行
• 稀疏矩阵求解往往是资源瓶颈
时间瓶颈，内存，外存等瓶颈
2）数值分解：
3）前代 回代：
I.S.Duff, A.M.Erisman, and J.K.Reid. Direct Methods for Sparse Matrices. London:Oxford Univ. Press,1986. J.J.Dongarra,I.S.Duff, ... Numerical Linear Algebra for High-Performance Computers. G.H.Golob, C.F.Van loan. Matrix Computations. The Johns Hopkins University Press. 1996
为何需要密集子块(Dense Matrix)
http://eigen.tuxfamily.org/index.php?title=Benchmark
多波前法(multifrontal)简介
• 发展
Duff and Raid [2] J.W.H.Liu等分析，改进 [3] T.A.Davis开发UMFPACK [4]
2 密集子阵
根据符号分析，数值分解算法的不同，直接法有以下几类[15]： 1）general technique(通用方法)： 主要采用消去树等结构进行LU分解。适用于非常稀疏的非结构化矩阵。 2）frontal scheme(波前法)
LU分解过程中，将连续多行合并为一个密集子块(波前)，对 这个子块采用BLAS等高效数学库进行分解。 3）multifrontal method(多波前法)
f,g,h gg·
h ··
{h,i,j}
消去树
bej e ·· j ··
e,i,j i ·· j ··
h,i, j ii· j ·j
d,i,j i ·· j ··
消去树
消去树是符号分析的关键结构，其每个 节点对应于矩阵的一列（行），该节点只 与其父节点相连，父节点定义如下：
将符号结构相同的多行(列) 合并为多个密集子块，以这些密 集子块为单位进行LU分解。这些子块的生成，消去，装配，释放等都需 要特定的数据结构及算法。
稀疏排序
排序算法的作用是减少矩阵LU分解过程中产生的注入元，求解 矩阵的最优排序方法是个NP完全问题（不能够在合理的时间内进行 求解），对具体矩阵而言，目前也没有方法或指标来判定哪种算法 好。因此实测不同的算法，对比产生的注入元，是确定排序算法的 可靠依据。
GSS1.2
15秒
2012
3.0G 4核
GSS 2.3
4秒
2014
i7 8核
GSS 2.4
1秒
• 硬件的发展 CPU，内存等
• 稀疏算法逐渐成熟
稀疏排序，密集子阵 multifrontal ,supernodal…
• 数学库
BLAS，LAPACK
www.netlib.org
稀疏LU分解算法的关键
稀疏LU分解 （不考虑零元的）LU分解 1 零元是动态的概念，需要稀疏排序,减少注入元(fill-in)
LU分解形成frontal
10阶矩阵。 蓝点代表非零元。红点表示分解产生的注入元(fill-in) Frontal划分{a}, {b}{c}{d} {e} {f,g}{h,i,j}
Frontal的装配，消去，更新过程
{a}
ach c ·· h ··
{c} {b}
c,g,h g ··
h ··
{f,g} {e} {d}
主要的排序算法有最小度排序（MMD，minimum degree ordering algorithm）和嵌套排序（nested dissection）两种。矩阵排序方面相应 的成熟软件库有AMD[12] 、COLAMD、METIS[13]等。
Yannokakis M. Computing the minimum fill-in in NP-Complete. SIAM J.Algebraic Discrete Methods, 1981, 2:77~79
求解器的飞速发展
• BBMAT
http://www.cise.ufl.edu/research/sparse/matrices/Simon/bbmat.html
38744阶，分解后元素超过四千万.
1988 巨型机cray-2上
>1000秒
2003
4G
umfpack4
32.6秒[4]
2006
3.0G
稀疏矩阵复杂、多变
• 基本参数
对称性，稀疏性，非零元分布
• 敏感性，病态矩阵
条件数
• 格式多变
Harwell-Boeing Exchange Format 。。。
• 测试集
Harwell-Boeing Sparse Matrix Collection UF sparse matrix collection
近似最小度排序算法 – 商图
近似最小度排序（AMD，Approximate Minimum Degree OrderingAlgorithm）算法于 1996年左右由Patrick R. Amestoy等学者提出
AMESTOY, P. R., DAVIS, ... 1996a. An approximate minimum degree ordering algorithm. SIAM J. Matrix Anal. Applic. 17, 4, 886–905.
• 基本算法
利用稀疏矩阵的特性，得到一系列密集子阵（波前）。将LU分解 转化为对这些波前的装配，消去，更新等操作。
• 多波前法的优点
波前是dense matrix ,可直接调用高性能库（BLAS等） 密集子阵可以节省下标存储 提高并行性
• 目前主要的求解器
UMFPACK,,GSS,MUMPS,PARDISO, WSMP,HSL MA41等
大型稀疏矩阵 的LU分解及特征值求解
陈英时 2016 . 1. 9
2013. 7. 20
稀疏矩阵求解的广泛应用
• 矩阵求解是数值计算的核心[1] • 稀疏矩阵求解是数值计ຫໍສະໝຸດ Baidu的关键之一
偏微分方程，积分方程，特征值，优化… 万阶以上dense matrix不可行
• 稀疏矩阵求解往往是资源瓶颈
时间瓶颈，内存，外存等瓶颈