大数据计算方法 课件 chap4-稀疏矩阵解法

1大数据计算方法Chap4：大型稀疏线性方程组的求解稀疏矩阵基本概念稀疏矩阵的分解算法高效率的稀疏矩阵直接解法迭代解法与共轭梯度法广义最小残量法(GMRES)Outline2稀疏矩阵基本概念基本概念存储格式简单的处理算法3稀疏矩阵(sparsematrix)存在大量零元素的矩阵Wilkinson’sdefinition:“matricesthat

allowspecialtechniquestotakeadvantage

ofthelargenumberofzeroelements.”很多实际仿真问题(电路网络,偏微分方程...)矩阵都是稀疏的非零元数、稀疏度、稠密度带状矩阵(bandmatrix):基本概念>>density=nnz(A)/prod(size(A))>>sparsity=1-density

带宽:非零元到主对角线的最大距离>>[i,j]=find(A)>>bandwidth=max(abs(i-j))4(也称之为”半带宽”)非零元数目为O(n)(不处理零元素)spy看非零元分布情况存储格式一般稀疏矩阵的存储结构三元组(COO格式)非零元的排列顺序可任意;压缩稀疏行(CSR格式)按第1行,第2行,…顺序存非零元;压缩稀疏列

aa1253467891011row11222333445coa1253467891011corow13691112

方便做逐行逐行

执行的操作5冗余:可能连续出现相同行编号prow为各行第一个元素位置aa1364927510811row12324132435pcol14681012

(CSC格式)用二维数组存储的矩阵不是稀疏矩阵!Matlab中的稀疏矩阵采用压缩稀疏列(CSC)方式存储与稠密矩阵(二维数组)转换若阶数很大,稠密格式A无法存

可直接用sparse命令生成

spdiags:“带状”稀疏矩阵sprand,sprandn生成随机的稀疏阵>>whosNameSizeBytesClassAttributesA20x203200doubleS20x20528doublesparse基本的Matlab函数>>S=sparse(i,j,x,m,n)S=spdiags([abc],[-1,0,1],n,n)6>>S=sparse(A)>>A=full(S)(相当于用COO格式输入数据)生成的S满足whos看变量类型(相当于用DIA格式输入数据)带状矩阵的分解与矩阵向量乘基本原则通过不操作0减少计算量(与0的加减乘除都没有增加新信息)例如：A+B的浮点运算次数不超过LU/Chol等复杂的矩阵分解,应设计考虑稀疏矩阵的特殊算法三对角矩阵方程Ax=b

tridisolve.m,用3个一维数组存A71.Fork=1ton-12.Fori=k+1ton3.4.Forj=k+1ton5.6.End7.End8.End主元乘子ak,k+2=0复杂度O(n)!min(nnz(A),nnz(B))LU分解带状矩阵的分解与矩阵向量乘带状矩阵的LU分解仍然保持原始带宽

很多应用问题中的矩阵A对角占优/对称正定，不需要选主元倘若需要选主元，带宽增大，最多使带宽达2

+181.Fork=1ton-12.Fori=k+1tok+

-13.4.Forj=k+1tok+

-15.6.End7.End8.End~线性复杂度(

<<n时)

00乘法次数:带状矩阵的分解与矩阵向量乘矩阵与向量乘法的算法矩阵存储采用CSC或其变种C语言算法实现(y=Ax)9jrow[k]

x

yj

Avoid

matvecC(cscptr

mata,double*x,double*y){

intn=mata->n,j,k,*ki;

double*kv;

for(j=0;j<n;j++)y[j]=0;

for(j=0;j<n;j++){

jval=mata->value[j];

jrow=mata->rowi[j];for(k=0;k<mata->nzcount[j];k++)y[jrow[k]]+=jval[k]*x[j];}

return;}structcsc{intn;int*nzcount;int**rowi;double**value;}*cscptr//第j列非零元行号数组//第j列的非零元素数组复杂度O(nnz(A))!j稀疏矩阵的直接解法直接解法概述Lsolve算法稀疏矩阵的LU分解稀疏矩阵的Cholesky分解消去树与并行算法1011高斯消去过程中的“填入”1.满足何条件时(i,j)处会发生fill-in?2.若非零元分布对称,发生fill-in的位置也对称?XXXXXXXXXX填入给LU/Chol分解带来麻烦？!nz=4991nz=348902基本问题——填入(fill-in)稀疏方程的直接解法概述高斯消元法(LU/Chol分解)针对稀疏矩阵的改进如何利用稀疏性减少内存用量和计算时间？减少内存用量稀疏矩阵存储格式(CSR,CSC,…)考虑矩阵非零元分布的变化(fill-in)，精确分配存储空间减少计算量忽略与0元素有关的运算填入元(fill-in)使矩阵“不稀疏”(也影响内存量)矩阵重排序(reordering)、预排序(preordering)难点--如何同时考虑减少填入与数值选主元?12Fill-reductionnumericalpivotingLsolve算法

13

xjx,b

l:,j=(须反映计算xi的先后依赖关系)算x3依赖x1算x4依赖x2算x8依赖x3,x5......

Lsolve算法14

(1)(2)

记为X=ReachG(B)顶点j顶点iLsolve算法15有向边说明了xi计算的依赖关系Lsolve算法16遍历所有路径，标记经过的点有向边说明了xi计算的依赖关系Lsolve算法17遍历所有路径，标记经过的点遍历某顶点的出边，相当于遍历矩阵的某列非零元(矩阵采用CSC存储)有向边说明了xi计算的依赖关系Lsolve算法18注意:采用DFS进行后序遍历(栈),得到X的顺序正确(DAG的拓扑顺序)有向边说明了xi计算的依赖关系Lsolve算法191234B={1}X={4}X={2,4}X={3,2,4}X={1,3,2,4}另一个例子DFS后序遍历得出的顺序是正确的(不违反xi之间的依赖关系)Lsolve算法20时间复杂度为O(flop);

flop可能比n和nnz(L)小得多j的邻点集合

(CSC存储结构)DFS后序遍历和栈稀疏LU分解算法左视LU分解算法(jki版本)21kij版本LU分解Right-looking

Left-looking被乘数乘子

22稀疏LU分解算法

只更新排列向量稀疏Cholesky分解算法23逐列计算的Cholesky分解算法

LLTAj由Left-lookingLU分解改造得来，只用、也只算下三角部分

利用Lsolve算法

或者，非零元位置预测可从循环中提出来

符号分解符号分解中的消去树24

顶点i顶点j顶点k

(j是i的双亲节点，节点n是树根)25符号分解中的消去树消去树的例子26符号分解中的消去树消去树的例子

etree_demo.m注意:对选主元的GPLU算法，由于交换行,不能用消去树27

28并行LU分解算法

快速计算符号分解反映了依赖关系的上限

(不受选主元影响)计算顺序反映依赖关系(比Lsolve算法中xi间依赖关系复杂)(串行)高效稳定的稀疏矩阵直接解法减少填入的矩阵重排稀疏对称正定阵的解法高效稳定的稀疏矩阵LU分解29集成电路供电线网仿真DC分析：求解电阻网络

电路中的节点电压一个IBM公司芯片设计的例子n=474,524nnz=2,020,882在Matlab中约占36MB内存“\”求解(SPD阵,内部用chol)

t

2s“R=chol(A)”:时间长,或内存

不够用!

?集成电路供电线网的一部分30?减少填入的重排—实例如何进行重排序？两个行/列对称重排的例子31Fill-ins交换1,2

行/列Nofill-in!(1,2,…,n-1,n)(n,2,…,n-1,1)减少填入的重排Nofill-in!找填入元最少的排序,

NP-hard!只有一些启发式算法

Permutationsforsparsity

“Iobservedthatmostofthecoefficientsinourmatriceswerezero;i.e.,thenonzeroswere‘sparse’inthematrix,andthattypicallythetriangularmatricesassociatedwiththeforwardandbacksolutionprovidedbyGaussianeliminationwouldremainsparseifpivotelementswerechosenwithcare”-HarryMarkowitz,describingthe1950sworkonportfoliotheorythatwon

the1990NobelPrizeforEconomicsCourtesyofJ.GilbertatUCSB32Markowitz算法1957年,定义填入元数目的上限:Markowitz积对未分解部分,计算任一非

零元素的Markowitz积=将该元素交换到主元位置

后,其Markowitz积不变将Markowitz积最小的换到主元位置例子:33乘

子

部

分已分解部分可能fill-in0000(当前行非零元数-1)

(当前列非零元数-1)列重排LU分解算法在akl(l

k)中寻找Markowitz积最小的akj若j

k,交换A的第k,j列用新的第k行执行消元操作减少填入的重排Markowitz算法类似可得行重排、对称重排Markowitz算法近似、局部最优,无法保证全局最优!对称稀疏阵的图(行~节点,非零元~边)图的邻接矩阵?无具体数值的稀疏矩阵34若不重排,5fill-in’s(执行列重排的例子)重排后,0fill-in!(从对角元中选)减少填入的重排Markowitz积=节点度212435对称稀疏正定阵的解法消去一行,图中增加团(clique)—fill-in对称Markowitz重排(MD算法):每次找图中度最小的节点AMD:近似最小度排序算法,效率高,效果类似(T.Davis提出)对称正定矩阵不需要数值选主元,重排序仅依赖于非零元结构。可以先做

对称重排，然后稀疏Chol分解算法

(符号分解+数值分解)35XXXXXXXXXXXXXXXXXXXXX预排序集成电路供电线网仿真DC分析：求解电阻网络

电路中的节点电压一个IBM公司的芯片设计例子n=474,524nnz=2,020,882在Matlab中,约占36MB内存“\”求解(对称正定矩阵),

t

2s“R=chol(A)”:时间长,或内存

不够用!

试[R,flag,p]=chol(A,'vector')集成电路供电线网的一部分36指示分解是否成功?

>>y=R’\b(P);>>

y=R\y;>>x(P)=y;对称稀疏正定阵的解法缺省amd重排

37高效稳定的稀疏矩阵LU分解TheHSLMathematicalSoftwareLibrary:

最大匹配问题(行交换)对电路仿真等领域,这可能是最好的方法各种减少填入的排序算法嵌套割集(Nesteddissection):矩阵图节点集合的割集

为S(尽量小),分割出两部分为P,Q,按这

三部分对节点编号(S最后);如此反复受的启发;AMD速度快，对中/大规模问题效果较好amd:快;symamd:慢带宽缩小排序(ReverseCuthill-McKee)使非零元向对角线集中;对细长结构的分析有效colamd:列重排AMD；colperm:基于各列非零元数目38Thebestgeneral-purposeorderingsareND/AMDhybridsPQS对非常大的例子效果较好symrcm命令;

见rcm_demo.m(METISdissect_demo.m高效稳定的稀疏矩阵LU分解(colamd_demo.m)(amd_demo.m)

39预排序数值稳定,填入减少

高效稳定的稀疏矩阵LU分解>>maxNumCompThreads(1)ans=4设置Matlab单线程矩阵A

(按如下顺序)求解算法稀疏的对角阵右端项元素除以矩阵对角元较稠密的带状方阵部分主元的带状矩阵LU分解(LAPACK包)上三角或下三角矩阵回代法或前代法(若是稀疏的,有特殊技巧)对三角矩阵作行排列形成的矩阵重排序后,用回代或前代法对称矩阵,且对角线元素大于零尝试Cholesky分解算法(稠密矩阵:LAPACK包,稀疏矩阵:CHOLMOD包);

若稠密、且不正定使用选主元的对称矩阵LDL分解算法(LAPACK包)稠密上Hessenberg矩阵消为上三角矩阵再回代求解一般的稀疏方阵针对稀疏矩阵的直接解法(UMFPACK包)一般的稠密方阵部分主元LU分解(LAPACK包)不是方阵通过矩阵的QR分解得到最小二乘解Matlab“\”的内部算法40TAMU,TimDavisSuiteSparse,Matrixcollection试试命令：dA=decomposition(A,type)x=dA\b41\--“It’samazinghowmuchnumericallinearalgebraiscontainedinthatsinglecharacter”Matlab中处理稀疏矩阵的命令生成或转换成稀疏矩阵sparse(X)orsparse(i,j,s,m,n);sparse(m,n)构造特殊的稀疏矩阵speye(n,m);

spones(pattern);spdiags(B,d,m,n)随机稀疏矩阵生成器sprand(S)orsprand(m,n,density);sprandn统计非零元数目，获取非零元信息nnz(X);nonzeros(A),find(X);nzmax

(存储空间);

稀疏矩阵直接解法有关symbfact:

SPD阵的符号分解(返回非零元数目);etree:消去树

amd/colamd/symamd,dissect:重排序

\,

lu,chol,qr,svd等各种函数都支持稀疏矩阵42迭代解法与共轭梯度法基本思想Krylov子空间迭代法最速下降法共轭梯度法实用技术与比较43迭代解法的基本思想直接解法的缺点对于大的、稀疏矩阵A,计算代价大无法灵活地“tradeoff”准确度与计算时间迭代解法固定格式迭代Krylov子空间迭代代数多重网格(AMG)等固定格式迭代法误差变化：预条件矩阵的选择

44

(基本迭代法)

Jacobi:P=diag(A);G-S:P=tril(A);…迭代解法的基本思想

45

改善它的收敛性Krylov子空间迭代法

46

...最速下降法

47

几何意义?

最速下降法

或相对残差48

最速下降法

49推不出最优步长最优步长也没必要能解一般优化问题共轭梯度法

改进

之处:

50

共轭梯度法

51

共轭梯度法算法描述收敛性:

一次矩阵向量乘法

两次向量内积收敛所需步数

n

52

共轭梯度法

53Y.Saad,“IterativeMethodsforSparseLinearSystems,”SIAMPress,2003共轭梯度法54

PCG算法CG算法方案1方案2

实用技术

opt.type='ict';55

opt.type=‘nofill';opt.droptol=?广义最小残量法(GMRES)投影方法Arnoldi过程GMRES算法改进与辅助技巧56Krylov子空间与投影方法

57

Arnoldi过程

58=

A

O

即Arnoldi过程

59如何得到

的基向量V?(类似于Gram-Schmidt正交化过程)可用Arnoldi-MGS代替

Arnoldi过程矩阵形式Arnoldi过程每一循环步生成矩阵的一列60......结果矩阵维度为n

m扩一行的上Hessenberg阵m+1nm

Hm为方阵，是上Hessenberg矩阵

GMRES算法

61Dountil足够小

End如何简化计算？m该如何取值？采用Arnoldi过程构造Km的基向量矩阵Vm

,

其中Solve,其中

(最小残量法)求解对应的系数矩阵右端项n

m(m+1)

mGMRES算法

62

扩大的上Hessenberg矩阵GMRES算法

63

**00残差为0，得到准确解！WenjianYu