大型稀疏矩阵的LU分解及特征值求解
稀疏矩阵方程算法
稀疏矩阵方程算法稀疏矩阵是指矩阵中绝大多数元素为0的矩阵。
在实际问题中，很多矩阵都是稀疏的，例如图像处理、自然语言处理等领域。
由于稀疏矩阵的特殊性，传统的矩阵运算方法效率较低，因此需要设计高效的算法来解决稀疏矩阵方程。
稀疏矩阵方程是指形如Ax=b的线性方程，其中A是一个稀疏矩阵，b是一个向量。
解决稀疏矩阵方程的一种常用方法是使用迭代算法，例如共轭梯度法（Conjugate Gradient，CG）和广义最小残差法（Generalized Minimal Residual，GMRES）等。
共轭梯度法是一种迭代法，它可以用来解决对称正定稀疏矩阵方程。
该方法的基本思想是通过最小化残差的二次范数来逼近方程的解。
具体而言，共轭梯度法通过迭代计算一个与残差正交的搜索方向，并在该方向上进行搜索，直到找到方程的解。
广义最小残差法是一种迭代法，它可以用来解决非对称稀疏矩阵方程。
该方法的基本思想是通过最小化残差的二范数来逼近方程的解。
与共轭梯度法不同的是，广义最小残差法使用Krylov子空间来进行搜索，并在该子空间上进行最小化残差的计算。
除了迭代算法外，还可以使用直接求解方法来解决稀疏矩阵方程。
其中一种常用的方法是LU分解。
LU分解是将稀疏矩阵A分解为一个下三角矩阵L和一个上三角矩阵U的乘积，即A=LU。
通过LU分解，可以将原始方程Ax=b转化为Ly=b和Ux=y两个方程，进而求解出x的值。
稀疏矩阵方程的求解算法还有很多，例如Jacobi迭代法、高斯-赛德尔迭代法等。
这些算法在不同的问题和应用中具有不同的优势和适用性。
在实际应用中，稀疏矩阵方程的求解是一个复杂且关键的问题。
通过选择合适的算法和优化技术，可以提高求解的效率和精度。
同时，还可以利用稀疏矩阵的特殊性质，例如压缩存储和并行计算等，进一步提高算法的性能。
稀疏矩阵方程是一类特殊的线性方程，传统的矩阵运算方法在处理稀疏矩阵时效率较低。
针对稀疏矩阵方程，可以采用迭代算法和直接求解方法来求解。
矩阵的lu分解计算量
矩阵的LU分解是一种将一个矩阵分解为一个下三角矩阵（L）和一个上三角矩阵（U）的乘积的方法。
LU分解是线性代数中一种重要的矩阵分解方法，它在求解线性方程组、计算行列式、计算矩阵的逆等方面都有广泛的应用。
对于一个$n \times n$的矩阵A，其LU分解的计算量主要取决于以下几个因素：
1. 存储量：LU分解需要存储原始矩阵A、下三角矩阵L和上三角矩阵U，因此需要额外的$3n^2$个浮点数存储空间。
2. 计算量：LU分解需要进行一系列的矩阵乘法和行交换操作，因此计算量相对较大。
具体来说，LU分解的计算量主要包括以下几部分：
* 计算L矩阵：需要执行$n(n+1)/2$次乘法操作，其中$n$是矩阵的阶数。
* 计算U矩阵：需要执行$n^3/6$次乘法操作，其中$n$是矩阵的阶数。
* 行交换操作：需要执行$n-1$次行交换操作，其中$n$是矩阵的阶数。
因此，LU分解的计算量大约为$O(n^3)$，其中$n$是矩阵的阶数。
在实际应用中，为了提高计算效率，通常会采用一些优化算法和并行计算技术来加速LU分解的计算过程。
lu分解方法求矩阵特征值
lu分解方法求矩阵特征值
矩阵特征值可以通过LU分解方法来求解。
首先，我们需要将矩阵A进行LU分解，得到一个下三角矩阵L和一个上三角矩阵U，使得LU=A。
接下来，我们可以利用LU分解的结果来求解特征值。
假设我们已经得到了LU分解后的矩阵L和U，那么我们可以进行如下步骤来求解矩阵A的特征值：
1. 首先，求解U矩阵的特征值。
U矩阵是一个上三角矩阵，其特征值就是它的对角线元素。
2. 接下来，求解L矩阵的特征值。
L矩阵是一个下三角矩阵，其特征值也是它的对角线元素。
3. 最后，将U矩阵和L矩阵的特征值组合起来，就得到了矩阵A的特征值。
需要注意的是，LU分解方法求解特征值的过程相对比较复杂，尤其是对于大型矩阵而言。
在实际应用中，可以借助计算机软件来进行LU分解和特征值求解，以提高计算的准确性和效率。
总之，通过LU分解方法可以求解矩阵的特征值，但需要注意计算的复杂性和精度的要求。
希望这个回答能够帮助到你理解如何利用LU分解方法来求解矩阵的特征值。
lu分解算法
lu分解算法
LU分解算法是一种将一个非奇异矩阵分解成一个下三角矩阵L和一个上三角矩阵U的方法，它可以用于解线性方程组以及求矩阵的逆等计算中。
具体的LU分解算法如下：
输入：一个n×n的非奇异矩阵A
输出：下三角矩阵L和上三角矩阵U
1. 初始化一个n×n的下三角矩阵L和一个n×n的上三角矩阵U，使它们的所有对角元素为1。
2. 对于矩阵A的第一行，将其作为矩阵U的第一行。
3. 对于矩阵A的第一列，将其除以矩阵U的第一个元素得到矩阵L的第一列。
4. 对于矩阵A的剩余行，以及对应的列，进行如下操作：
- 计算当前元素的值，即A(i, j)减去矩阵L的第i行与矩阵U的第j列的内积。
- 如果i小于等于j，将计算得到的值赋给矩阵U的第i行第j列元素。
- 如果i大于j，将计算得到的值除以矩阵U的第j列第j个元素，然后赋给矩阵L的第i行第j列元素。
5. 返回矩阵L和矩阵U作为结果。
通过LU分解算法，可以将解线性方程组的计算转化为简单的矩阵乘法和求解步骤。
此外，通过求解LU分解后的矩阵，还可以求矩阵的逆和行列式等相关计算。
大型稀疏矩阵的LU分解及特征值求解-bestparallelsparsesolver
利用稀疏矩阵的特性，得到一系列密集子阵（波前）。将LU分解 转化为对这些波前的装配，消去，更新等操作。
• 多波前法的优点
波前是dense matrix ,可直接调用高性能库（BLAS等） 密集子阵可以节省下标存储 提高并行性
• 目前主要的求解器
UMFPACK,,GSS,MUMPS,PARDISO, WSMP,HSL MA41等
GSS1.2
15秒
2012
3.0G 4核
GSS 2.34秒ຫໍສະໝຸດ 2014i7 8核
GSS 2.4
1秒
• 硬件的发展 CPU，内存等
• 稀疏算法逐渐成熟
稀疏排序，密集子阵 multifrontal ,supernodal…
• 数学库
BLAS，LAPACK
稀疏LU分解算法的关键
稀疏LU分解 （不考虑零元的）LU分解 1 零元是动态的概念，需要稀疏排序,减少注入元(fill-in)
• 直接法基于高斯消元法，即计算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
8、矩阵特征值问题计算
对应的特征向量x1, x2 ,, xm线性无关.
定理7(对称矩阵的正交约化 ) 设A R nn为对称矩阵 , 则
(1) A的特征值均为实数； (2) A有n个线性无关的特征向量； (3) 存在正交矩阵P使得
1 2 ， P 1 AP n 且i (i 1,2,, n)为A的特征值, 而P (u1,u2 , ,un )的列 向量u j为对应于 j的特征向量.
k
k
k A v0 max(vk ) max max(Ak 1v ) 0 k k 2 n 1 maxa1 x1 a2 x2 an xn 1 1 k 1 k 1 2 n maxa1 x1 a2 x2 an xn 1 1 1 (k )
k k 1
lim
vk
a1 x1.
即vk 是1的近似的特征向量. 而主特征值 (vk 1 ) j 1 n (vk 1 ) j 1 , 或1 . (v k ) j n j 1 (v k ) j
定理12 设A R nn有n个线性无关的特征向量, 其特征值
1 2 n ,
并设A的主特征值是实根，且满足
1 2 n ,
现在讨论求1及x1的基本方法.
(2.1)
v0 a1 x1 a2 x2 an xn , (设a1 0)
v1 Av0 a11 x1 a22 x2 ann xn ,
k k 2 n k vk Avk 1 1 a1 x1 a2 x2 an xn . 1 1 k 当k很大时，k 1 a1 x1, vk 1 1vk , Avk 1vk， v
大型稀疏矩阵直接求解算法的研究及实现共3篇
大型稀疏矩阵直接求解算法的研究及实现共3篇大型稀疏矩阵直接求解算法的研究及实现1大型稀疏矩阵直接求解算法的研究及实现随着计算机技术的不断发展和数学建模需求的增加，大型稀疏矩阵直接求解算法的研究和实现日益受到人们的关注。
在实际应用中，大型稀疏矩阵经常出现在各种科学计算、工程计算以及机器学习等领域。
因此，如何高效地求解大型稀疏矩阵成为了一个十分重要的问题。
一般来说，大型稠密矩阵的求解可以使用各种经典算法，如高斯消元、LU分解等。
然而，大型稀疏矩阵的求解却需要特殊的算法和数据结构。
传统的直接求解方法存在着效率低下和存储空间过大等问题，因此研究者们提出了许多改进方法和优化方案。
稀疏矩阵存储结构是求解算法中的重要问题之一。
目前，广泛应用的稀疏矩阵存储格式包括压缩列（Compressed Column，CC）、压缩行（Compressed Row，CR）以及双重压缩（Double Compressed）等。
这些存储格式各有优缺点，具体用哪一种存储格式取决于矩阵的具体特点和求解算法的需求。
比如，在随机梯度下降等机器学习算法中，常常使用压缩行存储方式来优化矩阵乘法操作的速度。
多核并行、GPU加速等技术也被广泛应用于大型稀疏矩阵的求解算法中，以提高计算效率。
并行求解算法可以将巨大的计算任务划分成多个子任务，并分配给多个核心同时执行，充分利用计算机的计算资源。
而GPU加速则充分利用了GPU的特殊架构，通过将计算任务映射到各个流处理器上并行执行，进一步提高求解效率。
除了以上所述的算法优化和技术应用，近年来还出现了一些新的求解算法。
比如，基于埃米尔特矩阵分解的求解算法，具有比传统LU分解更快的求解速度；基于内点法的求解算法，在高稀疏性的情况下，具有比其他算法更优的求解速度和精度。
综上所述，大型稀疏矩阵直接求解算法的研究和实现是一个充满挑战的领域。
在实际应用中，选择适合的算法和存储结构，并结合多核并行、GPU加速等技术，可以有效提高求解速度和精度。
[全]矩阵LU分解的几种算法
![[全]矩阵LU分解的几种算法](https://img.taocdn.com/s3/m/890d0cf2f111f18582d05a7c.png)
矩阵LU分解的几种算法Doolittle分解将矩阵A分解为单位下三角矩阵L和上三角矩阵UCrout 分解将矩阵A分解为下三角矩阵L和单位上三角矩阵UCholesky分解Doolittle分解和Crout 分解适于一般非奇异的矩阵，但对于一些更特殊的矩阵，我们有更好的分解方法。
基础概念矩阵A对称：A^T=A矩阵A正定：A的各阶顺序主子式大于0，对于实对称矩阵A正定的等价条件是A的特征值全为正假设矩阵A是对称正定矩阵，则可以分解为：其中L为下三角矩阵注：这里不给出证明，具体的分解过程，大部分数学软件都有相应的函数，我们更关心如何应用这样可以将求解线性方程组的过程看做两个步骤由于L为下三角矩阵，所以x,y都很好求解，简化了运算过程。
现在假设A为对称矩阵，去掉正定的条件，但是规定矩阵A的各阶顺序主子式不为0那么矩阵A可以做如下分解其中D为对角阵，L为下三角矩阵这样我们可以将求解线性方程组的过程同样看做两个步骤由于D为对角阵，它的逆就是它的倒数，其余的矩阵都是三角矩阵，所以计算也十分简便。
值得注意的是，显然，如果矩阵A是对称正定的，那么也是可以分解为LDL^T的，但如果矩阵A 不是正定的，那么不能分解为LL^T。
补充知识：一个三角矩阵的逆，也是三角矩阵且对角线上元素是倒数关系，但其余位置不是的。
例如：追赶法追赶法是针对带状矩阵（尤其是三对角矩阵）这一大稀疏矩阵的特殊结构，得出的一种保带性分解的公式推导，实质结果也是LU分解可以将一个三对角的稀疏矩阵分解为如下形式：其中三对角矩阵A为：最后提及一句：mathematica中提供LUDecomposition,CholeskyDecomposition 两个函数实现矩阵的LU分解。
LU分解法
具体的计算结果比较，取一正定对称矩阵,比较不 同算法 的结果，以及考虑稀疏和不考虑稀疏计算 时间的差别。
lu分解法是直接分解法中的一种算法将方程组axb中的稀疏矩阵a分解为一个上三角矩阵和一个下三角矩阵其中alu令yux那么在方程租的运算中可以先解lyb再解uxy在编程过程中分两步进行先对矩阵a进行lu分解然后再解方程组由lua及对l和u的要求可以得到分解的计算公式
LU分解法
线性方程组的解法通常分为两大类：直接接法和 迭代解法。
直接法是在没有舍入误差的情况下，通过有限步 四则运算可以求得方程组。但是，在实际计算时， 由于初始数据变为机器数而产生的误差以及计算 过程中所产生的舍入误差等都对解的精确度产生 影响，因此直接法实际上也只能算出方程组真解 的近似值。
目前较实用的直接法是高斯消去法和一些变形，例 如选主元的高斯消去法和矩阵的三角分解法，它们 都是目前计算机上常用的有效方法。直接法的优点 是计算量小，可事先估计，缺点是所需存贮单元较 多，编写程序较复杂，计算程序所需时间较长。
n2 11 12 13 21 31 22 23 32 33
u u U
(n-1)n nn
1n 2n
3n
…
n1
…
n3
…
a a
…
n2
a … a
nn
LU分解的函数程序为
• • • • • • • • • • • • • • function [L,U]=LU(A) r=length(A(1,:)) x=zeros(r,1) m=zeros(r,r) for i=1:r-1 for j=i+1:r m(j,i)=A(j,i)/A(i,i) for k=i:r A(j,k)=A(j,k)-m(j,i)*A(i,k) end end end L=m+eye(r) U=A
稀疏特征问题求解方法
紧缩 为了避免迭代收敛到已求得的特征对，这时需要采用紧缩处理。假
如已有Ritz向量u1,u2,…,uk ，带紧缩的校正方程为:
(I uiui*) I QkQk* (A iI ) I QkQk* (I uiui*)ti ri
0
nE 0
n
(
1
E)
0
Kx Mx
稀疏特征问题求解方法: Lanczos和Arnoldi方法
单/多向量（single and multiple vector）迭代法 幂(Power)迭代方法、逆迭代方法、Rayleigh商迭代(RQI)方法
Lanczos和Arnoldi方法 Lanczos和Arnoldi过程主要构造Krylov子空间的一个标准正交基：
预处理近似计算校正方程
校正方程可以使用线性系统的近似迭代法（像GMRES或
BCGSTAB）求解。如果结合适当的预条件子K A I 可以加速迭
代收敛过程。记
K (I ukuk* )K (I ukuk* ),
A (I ukuk* )( A k I )(I ukuk* )
m ( A, v) Span{v, Av, A2v,..., Am1v} 然后在这个子空间上计算近似的Ritz特征对。该过程是一个三项递 推求解过程,在每步仅需要一个矩阵-向量积。 基于不同的正交化方法：完全重正交、选择性重正交和局部重正交 等Lanczos和Arnoldi方法； 基于不同的重启动策略：显示重启动和隐式重启动Lanczos和 Arnoldi方法； 位移逆转换或谱转换：使用基本的Lanczos和Arnoldi方法总是收敛 于最大特征值，对于中间特征值需要使用位移-逆转换。另外，为
矩阵的LU分解
实验二矩阵的LU分解一、题目：求矩阵[2 1 1 2;1 2 3 2;2 4 1 1;3 1 2 3]的Doolittle分解二、方法：Doolittle分解，矩阵的紧凑格式的LU分解法三、程序：(1) 紧凑格式的LU分解法function a=nalupad(a)n=length(a);a(2:n-1)=a(2:n-1)/a(1,1);for k=2:na(k,k:n)=a(k,k:n)-a(k,1:k-1)*a(1:k-1,k:n);a(k+1:n,k)=(a(k+1:n,k)-a(k+1:n,1:k-1)*a(1:k-1,k))/a(k,k);end(2) LU分解function f=LU_decom(A)[m,n]=size(A)L=eye(n);U=zeros(n);flag='ok';for i=1:nU(1,i)=A(1,i);endfor r=2:nL(r,1)=A(r,1)/U(1,1);endfor i=2:nfor j=i:nz=0;for r=1:i-1z=z+L(i,r)*U(r,j);endU(i,j)=A(i,j)-z;endif abs(U(i,i))<epsflag='failure'return;endfor k=i+1:nm=0;for q=1:i-1m=m+L(k,q)*U(q,i);endL(k,i)=(A(k,i)-m)/U(i,i);endendLUend四、结果:(1) 紧凑格式的LU分解法>> nalupad([2 1 1 2;1 2 3 2;2 4 1 1;3 1 2 3])ans =2.0000 1.0000 1.0000 2.00000.5000 1.5000 2.5000 1.00001.00002.0000 -5.0000 -3.00003.0000 -1.3333 -0.4667 -3.0667(2) LU分解>> LU_decom([2 1 1 2;1 2 3 2;2 4 1 1;3 1 2 3])m =4n =4L =1.0000 0 0 00.5000 1.0000 0 01.00002.0000 1.0000 01.5000 -0.3333 -0.2667 1.0000U =2.0000 1.0000 1.0000 2.00000 1.5000 2.5000 1.00000 0 -5.0000 -3.00000 0 0 -0.4667五、拓展：矩阵分解成LU形式是有条件的，首先矩阵必须是非奇异的矩阵，其次矩阵的全部顺序主子式非零的时候才能完全保证矩阵可分解成LU且分解唯一。
复杂代数方程组的高效解法
QR分解法
1. 基本原理：介绍基本思想和步骤，让听众了解是如何进行的。 2. QR分解法的优点：客观分析QR分解法相较于其他代数方程组解法的优势和适用范围，如在处理大规模方程组时的高效性等等。 3. QR分解法的应用案例：通过实际案例，说明QR分解法在科学研究和工程领域的应用，如在计算机图像处理和金融统计建模中的应用等等。
矩阵分解算法概述
我们可以进一步探讨：
1. 基于LU分解的求解方法：利用高斯消元法将增广矩阵化为上三角矩阵， 再通过回代得到方程组的解；也可以通过LU分解得到两个三角矩阵L和U， 进而求解方程组。此方法对于系数矩阵为稠密矩阵且系数矩阵的行列式不为 的方程组求解效果最佳。 2. 基于QR分解的求解方法：通过将系数矩阵分解为正交矩阵Q和上三角矩 阵R来求解方程组。此方法适用于系数矩阵为满秩矩阵且系数矩阵的列数远 大于行数的方程组求解，例如线性回归问题。QR分解方法的求解速度一般 比LU分解方法快，但在解稠密方程组时相对较慢。
特征值分解的概念
特征值分解是一种用于矩阵处理的方法，它有以下几个方面的特点： 1. 定义和求解：特征值分解是将一个矩阵分解为若干特征向量和特征值的线性 组合，其中特征向量和特征值是通过求解线性方程组得出的。 2. 应用领域：特征值分解在众多的科学领域中得到了应用，包括计算机图形学、 量子力学、信号处理、机器学习等。在这些领域中，特征向量和特征值在对矩阵 进行变换和特征提取时起着至关重要的作用。 3. 数学意义：特征值分解的数学意义很重要，它将矩阵的变化和特征向量和特 征值联系起来，使得我们更好地理解和解释矩阵的行为和性质。同时，特征向量 和特征值的比较也可以帮助我们分类和推断数据中的模式。
Efficient Solution of Complex algebraic equation
大型稀疏矩阵的求解
• 引言 • 大型稀疏矩阵的存储方法 • 大型稀疏矩阵的求解算法 • 大型稀疏矩阵求解的优化技术 • 实际应用案例 • 总结与展望
01
引言
背景介绍
01
大型稀疏矩阵在科学计算、工程 技术和经济领域中广泛存在，如 有限元分析、电路分析、交通流 建模等。
02
由于其非零元素数量远小于总元 素数量，因此存储和计算效率成 为关键问题。
3
矩阵压缩存储的另一个重要原理是利用矩阵的稀 疏性，通过一些优化技术来进一步提高存储效率。
矩阵压缩存储的优点和局限性
01
矩阵压缩存储的主要优点是能够显著 减少存储空间的需求，特别是对于那 些稀疏度较高的矩阵。
02
此外，由于只存储了非零元素，因此 在某些计算过程中也可以减少计算量 ，提高计算效率。
03
算矩阵的逆或矩阵的行列式值来求解线性方程组。
02
直接法适用于小规模稠密矩阵，但对于大型稀疏矩阵
，由于计算量巨大，效率较低。
03
常见的直接法求解算法包括高斯消元法、LU分解法
等。
迭代法求解
迭代法是通过不断迭代逼近解的一种方法，适用 于大规模稀疏矩阵的求解。
迭代法利用矩阵元素的稀疏性，减少了计算量， 提高了求解效率。
预处理技术
预处理技术是指在求解稀疏矩阵之前，对原始数据进行一系列处理，以改善矩阵的性质，提高求解效率。常见的预处理方法 包括对角化预处理、不完全LU分解预处理等。
预处理技术可以减小矩阵的条件数，从而减小迭代方法的收敛速度，提高计算精度。同时，预处理技术还可以减少迭代过程 中的数值不稳定性，提高计算结果的可靠性。
理，从而显著缩短计算时间。
03
软件库和工具的开发
浅谈矩阵的LU分解和QR分解及其应用
浅谈矩阵的LU 分解和QR 分解及其应用基于理论研究和计算的需要，往往有必要把矩阵分解为具有某种特性的矩阵之积，这就是我们所说的矩阵分解.本文将介绍两种常用的矩阵分解方法，以及其在解线性方程组及求矩阵特征值中的应用.1.矩阵的LU 分解及其在解线性方程组中的应用 1.1 高斯消元法通过学习，我们了解到利用Gauss 消去法及其一些变形是解决低阶稠密矩阵方程组的有效方法.并且近些年来利用此类方法求具有较大型稀疏矩阵也取得了较大进展.下面我们就通过介绍Gauss 消去法，从而引出矩阵的LU 分解及讨论其对解线性方程组的优越性. 首先通过一个例子引入：例1,解方程组(1.1)(1. 2)(1.3)解. 1Step (1.1)(2)(1.3)⨯-+ 消去(1.3)中未知数,得到23411x x --=- (1.4)2Shep . (1.2)(1.4)+ 消去(1.4)中的未知数2x有12323364526x x x x x x ++=-=-=-⎧⎪⎨⎪⎩ 显然方程组的解为*x =123⎛⎫ ⎪ ⎪ ⎪⎝⎭上述过程相当于 111604152211⎛⎫ ⎪- ⎪ ⎪-⎝⎭~111604150411⎛⎫ ⎪- ⎪ ⎪---⎝⎭~111604150026⎛⎫⎪- ⎪ ⎪--⎝⎭2-（）+ ()i i r 表示矩阵的行由此看出，消去法的基本思想是：用逐次消去未知数的方法把原方程化为与其等价的三角方程组.下面介绍解一般n 阶线性方程组的Gauss 消去法.设111n n1nn a a a a A ⎛⎫ ⎪=⎪ ⎪⎝⎭ 1n x X x ⎛⎫ ⎪= ⎪ ⎪⎝⎭ 1n b b b ⎛⎫ ⎪= ⎪ ⎪⎝⎭则n 阶线性方程组AX b = (1.5)并且A 为非奇异矩阵.通过归纳法可以将AX b =化为与其等价的三角形方程，事实上： 及方程(1.5)为()()11A X b =,其中()1A A = ()1b b =(1) 设(1)110a≠，首先对行计算乘数()()11i1111i a m m =.用1i m -乘(1.5)的第一个方程加到第()2,3,,i i n =⋯个方程上.消去方程(1.5)的第2个方程直到第n 个方程的未知数1x .得到与(1.5)等价的方程组()()()11n 12n 111nn 0a a x x a ⎛⎫⋯⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⋯⎝⎭⎝⎭=()()112n b b ⎛⎫ ⎪⎪ ⎪⎝⎭简记作()()22Ab = (1.6)其中()()()()()()211211111 ijij i ij i i i a m b b m a a b =-=- (2) 一般第()11k k n ≤≤-次消去，设第1k -步计算完成.即等价于()()k k AX b = (1.7)且消去未知数121,,,k x x x -⋯.其中()()()()()()()()()()1111112122222k k k k kk knk nknna n n a a a a a A a a a a ⎛⎫⎪ ⎪⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭设()0k kk a ≠计算()() (i=/1,,)k k ik ikkkaa k m n =+⋯,用()()(1,,)a n ikik k kki k n a m ==+⋯消去第1k +个方程直到第n 个方程的未知数k x .得到与(1.7)等价的方程组()()1k 1k A X b ++= 故由数学归纳法知，最后可以把原方程化成一个与原方程等价的三角方程组.但是以上分析明显存在一个问题，即使A 非奇异也无法保证()0i ii a ≠,需要把非奇异的条件加强.引理1 约化主元素()01,,)i ii a k ≠=⋯（i 的充要条件是矩阵A 的顺序主子式0i D ≠.即1111110,0ikk kkk a a D a a D a =≠=≠⋯证明 利用数学归纳法证明引理的充分性.显然，当1k = 时引理的充分性是成立的，现在假设引理对1k -是成立的，求证引理对k 亦成立.有归纳法，设()()01,21iii i a k ≠=⋯-于是可用Gauss 消去法将中，即()()()()()()()()()()()11111121n22222n 1k k k k k kk knnknn a a a a a A a a a a A ⎛⎫ ⎪ ⎪⎪ ⎪→= ⎪ ⎪ ⎪ ⎪⎝⎭即()()()()()()()()11121231112211223112233222a a D a D a a a a a ===()()()()()()()()()11111122212222k 11122k k k kkk kka a a a a a D a a a ==⋯ (1.8) 由设0(1,,)i i D k ≠=⋯及式(1.8)有()0k kk a ≠显然，由假设()()01,2iiii k a ≠=⋯，利用(1.8)亦可以推出0(1,,)i i D k ≠=⋯ 从而由此前的分析易得；定理1 如果n 阶矩阵A 的所有顺序主子式均不为零，则可通过Gauss 消去法（不进行交换两行的初等变换），将方程组(1.5)约化成上三角方程组，即()()()()()()()()()1111111121122222222b b b n n n n nn n n a a a x x a a x a ⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪⎝⎭⎝⎭⎝⎭ (1.9) 1.2 矩阵LU 分解从而由以上讨论即能引出矩阵的LU 分解，通过高等代数我们得知对A 施行行初等变换相当于用初等矩阵左乘A ，即()()()()121211L A Lb A b == 其中 211n11101L m m ⎛⎫⎪- ⎪= ⎪⎪-⎝⎭一般第k 步消元，，相当于()()()()11k kk k k kL A A L b b ++==重复这一过程，最后得到()()()()11211121n n n n Ab L L L A L L L b --⎧⋯=⎪⎨⋯=⎪⎩ (1.10) 其中k 1,111m 1n k k km L +⎛⎫ ⎪ ⎪ ⎪=⎪- ⎪ ⎪ ⎪ ⎪-⎝⎭将上三角形矩阵()n A U 记作,由式(1.9)得到111121=U n A L L L LU ----⋯=,其中211111211211m 1n n n m L L m L L ----⎛⎫⎪⎪=⋯= ⎪⎪⎝⎭由以上分析得；定理2 （LU 分解） 设A 为n 阶矩阵，如果A 的顺序主子式i 0(1,2,,1)D i n ≠=-.则A 可分解为一个单位下三角矩阵L 和一个上三角矩阵U的乘积，且这种分解是唯一的.证明 由先前的分析得出存在性是显然的，即A LU =.下证唯一性,设A LU CD == 其中L , C 为单位下三角矩阵，U ,D 为上三角矩阵.由于1D -11D C L U --=上式右端为上三角矩阵，左端为单位下三角矩阵，从而上式两端都必须等于单位矩阵，故U D =,L C =.证毕.例2 对于例子1 系数矩阵矩阵111041221A ⎛⎫ ⎪=- ⎪ ⎪-⎝⎭由Gauss 消去法，得结合例1，故100111010041211002A LU ⎛⎫⎛⎫⎪⎪==- ⎪⎪ ⎪⎪--⎝⎭⎝⎭对于一般的非奇异矩阵，我们可以利用初等排列矩阵kki I (由交换单位矩阵I的第k 行与第k i 行得到)，即()()()()()()()()111212111111,,kk k k k ki k i k k i i k A L b L I A I b L I A I b A L b++⎧==⎪⎨==⎪⎩ (1.11) 利用(1.11)得()1111,11n nn n i i L I L U I A A ---==.简记做.其中下面就n 情况来考察一下矩阵.()()4321444343544332211443443243)(i i i i i i i i i i I L I L I L I A I L I I L I I A A L L I ===⨯4324324321432()i i i i i i I I I L I I I 43214321 )(i i i i I I I I A从而记从而容易的为单位下三角矩阵,总结以上讨论可得如下定理.定理3 如果A 非奇异矩阵，则存在排列矩阵P 使PA LU = 其中L 为单位下三角矩阵，U 为上三角矩阵.1.3 矩阵LU 分解的应用以上对非奇异矩阵A 的LU 分解进行了全面的讨论，一下我们就简单介绍一下应用.对于矩阵A 一旦实现了LU 分解，则解线性方程的问题,便可以等价于：(1)Ly b = 求y (2)=Ux y , 求x (1.12)即，设A 为非奇异矩阵，且有分解式A LU =，其中L 为单位下三角矩阵，U 为上三角矩阵。
大型稀疏线性方程组符号LU分解法
Ｃ ｍ ｕ ｒＥ ＃ｎ ｅ ｎ ｎ ｐ ｌａｉ ｓ 算机 工程 与应用 ｏ ｐ ｔ ｎ ｅｈ ｇ ａｄ Ａ ｐｉ ｔ ｎ 计 ｅ ｃ ｏ
大型稀疏线性 方程 组符号 Ｌ 分解 法 Ｕ
张永 杰 ， 孙 秦
Ｚ Ａ Ｇ Ｙ ｎ －ｉ，Ｕ ｉ Ｈ Ｎ ｏｇｊ Ｓ Ｎ Ｑｎ ｅ
西北工业大学 航空学院 ， 西安 ７ ０ ７ ０２ １
Ｓ ｈ ｏ ｆ Ａｅ ｏ ａ ｔ ｓ ＮＰ Ｘｉ ａ １ ０ ２ Ｃ ｉａ ｃ ｏ ｌ ｏ ｒ ｎ ｕ ｉ ， Ｕ， ’ ｎ ７ ０ ７ ， ｈ ｎ ｃ
Ｅ ｍａ ：ｙｌ １ １ ｗ ｕｅ ｕｃ － ｉ ｚｊ９ ９ ＠ｎ ｐ ． ． ｌ ｄ ｎ Ｚ ＮＧ Ｙ ｎ －ｉ，Ｕ Ｑｉ￣ｙ ｏ Ｕ ｅｏ ｏｉｏ ｔｏ ｆｌｒｅｓａｅｓａ ｓ ｎ ａ ｑ ａ ｏ ｓ ｏ ｕｅ ｎ ｉｅｒ ｇ ＨＡ ｏ ｇ ｊ Ｓ Ｎ ｎ ｍｂ ｌ ｅ Ｌ ｄ ｃｍｐ ｓ ｉｎ ｍｅｈ ｄ ｏ ｇ ｃｌ ｐ ｒｅｌ ｅｒｅ ｕ ｔ ｎ． ｍｐ ｔｒＥ ｇｎ ｅｉ ｔ ａ ｉ ｉ Ｃ ｎ ａ ｄ Ａｐ ｌａｉｎ ，０ ７ ４ （８）２ － ０ ｎ ｐｉ ｔ ｓ ２ ０ ，３ ２ ：９ ３ 。 ｃ ｏ
有限元稀疏线性方程组的求解效率很大程度上取决于稀疏矩阵的存贮策略为了达到节省存贮空间和便于存取运算的目的先后出现了对角元存贮法等带宽存贮法变带宽存贮法坐标存贮法ellpackitpack存贮法csr存贮法shermans存贮法是一种结构简单的链表式存贮充分利用了矩阵的稀疏性对称性特点不仅存贮规模小而且便于元素查询插入和删除操作对于提高方程组的求解效率十分有利
关键词 ： 大型稀疏线性方程组 ； 全稀疏存贮策略 ； 符号 Ｌ Ｕ分解
文章编号 ：０ ２ ８ ３ （０ ７２ —０ ９ ０ 文献标识码 ： 中图分类号 ：Ｐ ０ ． １０ — ３ １２ ０ ）８ ０ ２ — ２ Ａ Ｔ３１ ６
增量法求解大规模稀疏矩阵方程组
增量法求解大规模稀疏矩阵方程组近年来，随着人工智能、大数据分析和科学计算等领域的快速发展，大规模稀疏矩阵方程组的求解变得尤为重要。
在实际应用中，由于矩阵规模大、非零元素分布稀疏等特点，传统的直接或间接法求解效率低下，甚至无法胜任。
增量法成为了求解大规模稀疏矩阵方程组的一种有效手段。
在本篇文章中，我们将从简入深地介绍增量法的基本原理、常见算法以及应用实例，帮助读者全面、深刻地理解增量法在求解大规模稀疏矩阵方程组中的重要作用。
1. 基本原理大规模稀疏矩阵方程组中的“大规模”指的是矩阵的规模非常大，通常是上万甚至上亿级别；“稀疏”则表示矩阵的非零元素分布非常稀疏，大部分元素为零。
在实际问题中，这样的矩阵可能来自于网络分析、有限元法、信号处理等领域。
传统的直接法求解（如高斯消元法）由于需要处理大量的零元素而效率较低，而间接法求解（如迭代法）则常常陷入收敛速度慢、计算精度不高等问题。
增量法的基本原理是将原始的大规模稀疏矩阵方程组逐步化为一个个规模较小的子问题，并利用迭代或直接法求解这些子问题，最终得到原问题的解。
通过这种方式，增量法在一定程度上克服了传统方法的缺点，提高了求解效率和精度。
2. 常见算法在增量法求解大规模稀疏矩阵方程组中，有几种常见的算法被广泛应用，包括但不限于以下几种：2.1 LU分解LU分解是一种经典的增量法求解大规模稀疏矩阵方程组的方法。
它将原始矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积形式，从而将原问题转化为两个规模较小的子问题。
可以通过前代和回代的方式分别求解这两个子问题，从而得到原问题的解。
2.2 Cholesky分解Cholesky分解是一种针对对称正定矩阵的增量法求解方法。
与LU分解类似，Cholesky分解也是将原始矩阵分解为一个下三角矩阵和其转置矩阵的乘积形式。
通过这种分解，可以将原问题转化为一个规模较小的子问题，并利用直接法求解得到原问题的解。
2.3 前代和回代前代和回代是一种常见的增量法求解大规模稀疏矩阵方程组的迭代方法。
大型稀疏矩阵的求解
稀疏矩阵的求解方法及优缺点
各方法的优缺点：
高斯消去法与直接三角分解法：计算量大致相同，约为
1 3 n 3
乘除法运
算。三角分解的优越性在于：如果实现了A的三角分解A=LU，则解方程组Ax=b等
价于解方程组LUx=b。令Ux=y，则解Ax=b等价于解Ly=b，Ux=b。这是两个三角形 方程组，很容易求解。直接法精确性较高，对矩阵本身性质要求较低，所以适用 性比较广，但要防止舍入误差过大地影响结果，导致解的失真。 Jacobi迭代法，Gauss-Seidel迭代法，SOR迭代法：迭代法对矩阵的本身性质要 求比较高，但对于大型矩阵，由于计算机有效位数的限制，直接法中的误差，消 元中有效位数的损失等将会影响方程求解精度，所以这时迭代法比较适用。
，则得解，输出
x
(k )
若在k步时无法使
x ( k ) x ( k 1) ，则返回第二步。
4）若总步数k到N时还得不到满足误差的解，则算法失败，k步内无法得 精确解。
2、软件的应用
本次使用的软件是我们组员共同用VB语言自己编写的程序，我负 责的部分是Jacobi法，具体操作界面如下：
aii 0 ，i=1，2， ，n。）
( k 1) (b1 a12 x2 ( k ) a13 x3 ( k ) a1n xn ( k ) ) x1 a11 ( k 1) 1 (b2 a21 x1( k ) a23 x3 ( k ) a2 n xn ( k ) ) x2 a22 1 ( k 1) xn (bn an1 x1( k ) an 2 x2 ( k ) an ,n1 xn1( k ) ) ann
大型稀疏矩阵的概念及其特点
几类稀疏矩阵特征值的算法实现
几类稀疏矩阵特征值的算法实现吴承逊;谌稳【摘要】总结了几类适用于用迭代法求解稀疏矩阵特征值的算法.文章用到了幂法和反幂法,并在幂法算法的基础上对其进行规范化,说明了二分法求特征值的原理;实现了用二分法求对称三对角矩阵所有特征值的算法;讨论了能把大部分稀疏矩阵变成对称三对角矩阵的Lanczos方法,且对不同类稀疏矩阵使用Lanczos方法进行探讨,把二分法和Lanczos方法结合到一个算法中;并通过数值实验验证了这些算法的有效性.【期刊名称】《湖南城市学院学报（自然科学版）》【年(卷),期】2018(027)005【总页数】5页(P51-55)【关键词】稀疏矩阵;特征值;幂法;Lanczos法;算法实现【作者】吴承逊;谌稳【作者单位】湖南城市学院理学院,湖南益阳 413000;湖南城市学院理学院,湖南益阳 413000【正文语种】中文【中图分类】O242迭代法可以用来近似解决稀疏矩阵的特征值问题，它有存储空间小、程序简单的特点，是求解稀疏矩阵特征值的有效方法，其主要方法包括幂法及反幂法、二分法和Lanczos方法．幂法及反幂法是求稀疏矩阵特定特征值的常用方法，它适用性强，但是不能求出全部的特征值．随着J. W. Givens和A. S. Householder所研究的Givens变换和Householder变换的出现，使得任何实对称矩阵都可以通过这2种正交变换约化变成对称的三对角矩阵，于是二分法[1]得以出现，它的理论基础是Sturm序列．二分法是解决对称的三对角矩阵全部特征向量和特征值的有效方法．但是由于稀疏矩阵的阶数过于庞大，约化矩阵变得十分困难和复杂，且二分法不能很好的求解非对称矩阵．Lanczos[2]研究出了一种Lanczos方法，把稀疏矩阵化成对称三对角矩阵，这样稀疏矩阵的特征值问题就变成了转化的对称三对角矩阵的特征值问题，然后用QR方法或者二分法很容易就求解出来了．然而，实际计算中Lanczos向量很容易失去正交性，且长期是不稳定的．1970年之后，Peters Wilkinson、Golub等人重新研究了该方法，对Lanczos方法进行改进，于是该方法成为能够求解大型的稀疏矩阵特征值问题的有效方法[3]．近年来，关于Lanczos方法的研究主要在其应用方面，曾玮、赵永华[4]基于谱分割并行求解稀疏矩阵特征值，并将矩阵的特征值求解区间划分为多个独立的子区间，分别对各个子区间内的特征值进行独立求解﹒焦宝宝[5]用重正交化Lanczos法求解大型非正交归问题，计算结果表明了Lanczos法在Matlab编程和大型壳模型计算中的精确性和可行性﹒幂法是用来求解矩阵的主特征值及其对应的特征向量的迭代方法，其基本思路是：设Ａ是一个有着n个线性无关特征向量的实矩阵，于是它的特征值是，对应的特征向量是．假设Ａ的主特征值是实根，且有，任取1个非零初始向量，用A构造的一个向量序列为可称作是迭代向量．幂法算法的步骤如下：(1)初始化设置，随机生成稀疏矩阵A，是第1个元素为1，其余为0的且和A矩阵同阶的列向量；(2)计算，；(3)计算的绝对值，如果满足误差精度，就输出主特征值m；否则，下标加1，回到步骤(2)．运用幂法求解的主特征值m和它对应的特征向量v，若m的绝对值大于1或者小于1，就会使得迭代向量v的每个不等于０的分量随着k趋近于无穷大时而趋向于无穷大(或者趋向于0)，这样在电脑上用算法去实现时就可能会“溢出”，为了解决这个缺点，需要将幂法中的迭代向量进行规范化．对于迭代向量，将它规范化，可得到其中，max{}是指v的绝对值最大的分量．于是，规范化幂法就是按照公式(3)进行构造的向量序列，即由式(3)得到的最大特征值的近似值为，它对应的近似特征向量的迭代向量为．规范化幂法算法的步骤如下：(1)初始化设置，随机生成稀疏矩阵A，是第1个元素为1，其余为0的且和A矩阵同阶的列向量；(2)计算，；(3)计算的绝对值，如果满足误差精度，就输出主特征值；否则，下标加1，回到步骤(2)．反幂法是幂法的一种变形，它用来求解矩阵模最小的特征值及其对应的特征向量，也能用于计算给定的近似特征值的特征向量．幂法是求A的最大特征值的问题，而反幂法则是先求解的最大特征值，进而求得A的最小特征值的问题．反幂法的迭代公式可描述为：对于任意的初始向量，可以构造向量序列，即其中，迭代向量能够求解线性方程组得出；最小特征值为．当然，反幂法是可以用于计算出给定特征值的特征向量的．例如给定A的1个特征值a，其对应的特征向量为x，则可以对(p约等于a)应用反幂法，它的1个特征值是，这个值是所有特征值中最小的，于是就能够求出这个特征向量．同样，由于特征值值小，迭代次数往往也不大，几次就能迭代出特征向量．反幂法求最小特征值与幂法算法相似，只是在迭代公式上有些许不同，用反幂法求给定特征近似值的特征向量步骤如下：(1)初始化设置，输入稀疏矩阵A和给定的特征近似值p，是第1个元素为1，其余为0的且和A矩阵同阶的列向量；(2)计算出，以及；(3)求出，tol是误差精度，如果tol满足规定的条件，则上述循环结束，输出所求的；否则，下标加1，回到步骤(2)．此算法设计并非最优，因为算法中涉及到求矩阵的逆，可以对其进行如下改进：先把A进行LU分解，再通过求解，来进行算法的迭代．幂法及反幂法只能够求出可对角化矩阵A的最大和最小特征值，无法求出A的所有特征值，于是需要用到二分法．二分法的理论基础是Sturm序列，下面介绍下什么是Sturm序列．对于n阶对称的三对角矩阵用来表示的k阶顺序主子式，即令，且满足式(5)的关系：则称是B的Sturm序列，同时从公式(5)中很明显的看出来，是B的特征多项式．把中每相邻2元素符号相同的次数称为同号数，记为．那么，应用二分法的步骤就可以分为以下3步：(1)用Householder变换和Givens变换把所求矩阵约化为对称三对角矩阵；(2)再求出三对角矩阵的顺序主子式符号，用来确定中根的个数，式中分别表示A 的r阶顺序主子阵和单位矩阵的r阶主子，然后再用区间的二分法隔离出所需要的根；(3)计算所求矩阵的相应特征向量，可以通过反幂法求出．二分法的关键是可以根据Sturm序列的同号数来判断特定区间的根的个数，也就是说就是在区间中根的个数，即B中不小于的特征值的个数．特征值问题就是求解式(6)的问题，即其中，A是一个非对称实矩阵；、分别是A的特征值所对应的右特征向量和左特征向量．Lanczos方法就是生成了2组双正交向量和，使得，其中是指n阶的单位矩阵，满足，且由上面2个公式可以得出将式(7)以列的形式描述出来，得到对其进行重新排序，得到递推公式当然，系数不是唯一的，只要两者之一满足即可，再引入双正交的参数，取、以及剩下的系数，其中可以由方程得到．上述Lanczos方法有着显著的优势，它把非对称实矩阵约化成了对称三对角矩阵，上面描述的迭代过程最多进行n次就终止了，除了元素间可能出现的正负号之差外．本文对Lanczos方法的算法[6]，作出一点改进，把它和二分法结合在一起，设计出一个新算法，它们两者是承上启下的，这个算法的优势在于最后输出结果是对称三对角阵、最大和最小特征值，具体算法步骤如下：(1)初始化设置，，；(2)从循环计算，且循环次数为n次(n为矩阵的阶数)．这样，矩阵就被上述循环算法化成了对称的三对角矩阵；(3)对转化的三对角矩阵使用二分法，然后按照公式(5)计算和同号数；再判断是否，如果是，则，否则；继续判断的值是否不小于误差，若是，则输出，否则就重新进行这个循环．可以使用Matlab随机生成稀疏矩阵，如生成1个100阶的稀疏矩阵，为了保证这个稀疏矩阵是一个非奇异矩阵，随机生成100阶矩阵后，取这个矩阵加上1个100阶的单位矩阵的和矩阵．随机选取4组不同阶数的稀疏矩阵来进行幂法的算法运算，稀疏矩阵的非0元素分布密度为0.3，运算结果如表1所示．表1的二范数q是的二范数，其中，A是所求的稀疏矩阵；是A的特征值；是对应的特征向量．显而易见，规范化幂法的算法是更优于幂法算法的，其改进程度也十分显著．从表1中可看出，幂法的缺点十分明显，它在前3组的稀疏矩阵迭代中迭代次数比规范化幂法的迭代次数多，差向量二范数也十分的大，说明它的特征向量中元素的值很大，这样计算会出现舍入误差，造成数值的不准确；在第4组稀疏矩阵的迭代中，它的特征值是1，明显是错误的，因为特征值为1，说明在特别稀疏的情况下，如果主特征值是很小的一个数，趋近于零，那么幂法就会失败．然而，规范化幂法就没有这样的缺点，它是一个更优于幂法的改进算法，能够适用于绝大多数的具有完备特征向量组的稀疏矩阵．随机选出了3个的10阶稀疏矩阵，用Lanczos算法来运算，运算结果如表2所示．从表2的运算情况来看，这个Lanczos算法较为完美的把稀疏矩阵化成了对称三对角矩阵，再通过把最大及最小值用幂法和反幂法验证，发现相差不大，证明了该方法有效．稀疏矩阵的特征值问题长期是矩阵研究中的一项重要任务，求解稀疏矩阵特征值的方法有很多是不稳定的、有误差的，虽然有着各种近似迭代方法的出现和发展，但是在每一种方法都有对应的特定情况，所以要证明哪种方法是最行之有效的，这是十分困难的．本文研究了幂法、规范化幂法、反幂法、二分法和Lanczos方法，发现能够解决稀疏矩阵特征值问题的主要方法还是Lanczos方法；而且，当用Lanczos方法运算后，用二分法求出的最大及最小特征值，有时候会和幂法与反幂法求出的特征值有一些误差．这是因为Lanczos方法迭代的双正交向量在多次迭代后会出现失去正交性的问题，求出的特征值不太理想，这个问题还有待进一步的研究和探讨．【相关文献】[1]黄云清, 舒适, 陈艳萍, 等. 数值计算方法[M]. 北京: 科学出版社, 2009.[2]LANCZOS C. An iteration method for the solution of the eigenvalue problem of linear differential and intergral operators[J]. Journal of Research of National Bureau of Standards, 1950, 45(4): 255-282.[3]周树荃. 计算大型稀疏矩阵部分特征值问题的有效方法[J]. 南京航空学院学报, 1981(1): 1-18.[4]曾玮, 赵永华. 基于谱分割的稀疏矩阵特征值问题并行求解[J]. 数值计算与计算机应用, 2015,36(2): 132-146.[5]焦宝宝. 用重正交化Lanczos法求解大型非正交归一基稀疏矩阵的特征值问题[J]. 物理学报, 2016, 65(19): (192101)1-8.[6]邓健新. 用Lanczos方法解高阶稀疏矩阵广义特征值问题的某些经验[J]. 数值计算与计算机应用, 1980, 1(3): 173-180.。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Frontal的装配，消去，更新过程
{a}
{c} {f,g} {b} {e} {d}
ach c ·· h ··
c,g,h g ·· h ··
bej e ·· j ··
d,i,j i ·· j ··
f,g,h gg· h ··
e,i,j i ·· j ··
{h,i,j}
消去树
h,i, j i i · j ·j
条件数
求解器的飞速发展
• BBMAT
/research/sparse/matrices/Simon/bbmat.html
38744阶，分解后元素超过四千万. 1988 巨型机cray-2上 2003 4G umfpack4 2006 3.0G GSS1.2 2012 3.0G 4核 GSS 2.3 2014 i7 8核 GSS 2.4
稀疏矩阵复杂、多变
• 基本参数 • 敏感性，病态矩阵 • 格式多变 • 测试集
Harwell-Boeing Exchange Format 。。。 Harwell-Boeing Sparse Matrix Collection UF sparse matrix collection
对称性，稀疏性，非零元分布
/index.php?title=Benchmark
多波前法(multifrontal)简介
• 发展
Duff and Raid [2] J.W.H.Liu等分析，改进 [3] T.A.Davis开发UMFPACK 一系列密集子阵（波前）。将LU分解 转化为对这些波前的装配，消去，更新等操作。
AMESTOY, P. R., DAVIS, ... 1996a. An approximate minimum degree ordering algorithm. SIAM J. Matrix Anal. Applic. 17, 4, 886 –905.
为何需要密集子块(Dense Matrix)
• 多波前法的优点
波前是dense matrix ,可直接调用高性能库（BLAS等） 密集子阵可以节省下标存储 提高并行性
• 目前主要的求解器
UMFPACK,,GSS,MUMPS,PARDISO, WSMP,HSL MA41等
LU分解形成frontal
10阶矩阵。 蓝点代表非零元。红点表示分解产生的注入元(fill-in) Frontal划分{a}, {b}{c}{d} {e} {f,g}{h,i,j}
• 直接法基于高斯消元法，即计算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
Yannokakis M. Computing the minimum fill-in in NP-Complete. SIAM J.Algebraic Discrete Methods, 1981, 2:77~79
近似最小度排序算法 – 商图
近似最小度排序（AMD，Approximate Minimum Degree OrderingAlgorithm）算法于 1996年左右由Patrick R. Amestoy等学者提出
大型稀疏矩阵 的LU分解及特征值求解
陈英时 2016 . 1. 9
2013. 7. 20
稀疏矩阵求解的广泛应用
• • 矩阵求解是数值计算的核心[1] 稀疏矩阵求解是数值计算的关键之一
偏微分方程，积分方程，特征值，优化… 万阶以上dense matrix不可行
• 稀疏矩阵求解往往是资源瓶颈
时间瓶颈，内存，外存等瓶颈
2 密集子阵
根据符号分析，数值分解算法的不同，直接法有以下几类[15]： 1）general technique(通用方法)： 主要采用消去树等结构进行LU分解。适用于非常稀疏的非结构化矩阵。 2）frontal scheme(波前法) LU分解过程中，将连续多行合并为一个密集子块(波前)，对 这个子块采用BLAS等高效数学库进行分解。 3）multifrontal method(多波前法) 将符号结构相同的多行(列) 合并为多个密集子块，以这些密 集子块为单位进行LU分解。这些子块的生成，消去，装配，释放等都需 要特定的数据结构及算法。
>1000秒 32.6秒[4] 15秒 4秒 1秒
• 硬件的发展 CPU，内存等 • 稀疏算法逐渐成熟
稀疏排序，密集子阵 multifrontal ,supernodal…
• 数学库
BLAS，LAPACK
稀疏LU分解算法的关键
稀疏LU分解 （不考虑零元的）LU分解
1 零元是动态的概念，需要稀疏排序,减少注入元(fill-in)
稀疏排序
排序算法的作用是减少矩阵LU分解过程中产生的注入元，求解 矩阵的最优排序方法是个NP完全问题（不能够在合理的时间内进行 求解），对具体矩阵而言，目前也没有方法或指标来判定哪种算法 好。因此实测不同的算法，对比产生的注入元，是确定排序算法的 可靠依据。 主要的排序算法有最小度排序（MMD，minimum degree ordering algorithm）和嵌套排序（nested dissection）两种。矩阵排序方面相应 的成熟软件库有AMD[12] 、COLAMD、METIS[13]等。