-
基于Spark框架的大型分布式矩阵求逆运算实现(二)——大型下三角矩阵求逆运算
2018-02-25 11:19:38先说一下原理: 对于一个可逆矩阵A,必然会得到它的唯一LU分解,即分解为一个下三角矩阵L和一个上三角矩阵U使得 A=L*U; 我们需要求得的问题是A的逆矩阵A`,已知 A=LU,A*A`=E,所以 A` = U`*L`; ...基于实际需要,需要对五百万阶的方阵进行求逆运算,但查看Spark(v. 2.2.0)的官方api并没有此方面的信息,就自己尝试着实现了一个;
先说一下原理:
对于一个可逆矩阵A,必然会得到它的唯一LU分解,即分解为一个下三角矩阵L和一个上三角矩阵U使得 A=L*U;
我们需要求得的问题是A的逆矩阵A`,已知 A=LU,A*A`=E,所以 A` = U`*L`;
由线性代数知识可知,一个下三角矩阵的逆必然也是下三角矩阵;
又因为,逆的转置等于转置的逆,即,(U`)T = (UT),而UT和L一样是下三角矩阵,可以实现代码复用;
所以问题便转化成了,
(一)大型可逆矩阵的LU分解
(二)大型下三角矩阵的求逆
第一部分由我的同学实现,之后会放出链接;这里主要讲一下大型下三角矩阵求逆的方法和实现;
大型矩阵运算,因为数据量过大,无法在单台计算机上进行,故需要进行并行化处理,这里采用分块矩阵乘法的思想。
首先设定一个步长 s,使得阶数为s的方阵可以在单台节点上进行求逆运算。
根据分块矩阵乘法,
L1*L1`=E
L2*L1`+L3*L2`=0
L3*L3`=E
化简可知,
L1`=(L1)`
L2`=-(L3)`*L2*(L1)`
L3`=(L3)`
如此一来。只要知道(L3)`,便可以知道整个矩阵的逆,而(L3)`同样是下三角矩阵,如此一来便可以进行迭代,当迭代到L3的步长不大于s时,便可以在单台节点上进行计算,如此一来,便可以反推回整个矩阵的逆;
下面进入实际实现部分:
1.基于Spark的api,将HDFS上的矩阵加载到内存中,类型为 BlockMatrix
2.调用BlockMatrix.blocks方法得到底层RDD,过滤出行标不大于列标的分块(下三角矩阵上半部分全是0,减少运算量)
3.首先得到原矩阵右下角的分块,求逆得到(L3)`,行标-1,得到L1和L2,得到(L1)`和L2,如此一来便可以拼凑出原矩阵右下角分块为2的矩阵的逆,迭代运算便可得到最终结果;
期间遇到的难点:
1.矩阵加载,Spark提供的原生api无法加载CSV文件直接转成BlockMatrix,所以此处进行了封装:
new IndexedRowMatrix(spark.sparkContext.textFile(path,Main.excutors).map(UDF.line2IndexedRow)) .toCoordinateMatrix().toBlockMatrix(steps, steps)
/* *输入一行以逗号(英文 , )分割的浮点数,最开始的数字作为索引 *返回一个IndexRow */ def line2IndexedRow(line:String): IndexedRow ={ val arrayBuffer = line.split(",").map(_.toDouble).toBuffer val index = arrayBuffer.head.toLong arrayBuffer.trimStart(1) val vector = Vectors.dense(arrayBuffer.toArray) new IndexedRow(index,vector) }
2.在计算 L2`=-(L3)`*L2*(L1)`时,由于直接调用矩阵分块乘法api会导致分块最终位置与算法设想不同,需要自行解决;
3.在本地运行时结果与集群运行结果不一致:由于算法全程使用尾递归进行迭代,有部分全局变量需要广播到各个节点;
4.性能优化,在矩阵运算过程中,由于是懒执行,部分变量会重复计算造成计算资源浪费,需要在SparkUI上查看,逐项调校;
5.Spark的persist机制:在调用RDD的persist方法后,RDD并不会马上被缓存,而是要等到第一个action调用时才会执行,但实际上
本算法中action的调用距离RDD首次生成相隔甚远,所以,需要在persist方法后接一个action来进行显示缓存;由于缓存项目过多
可能造成大量IO操作,需要及时进行unpersist操作;
优化后的RDD DAG截图如下:
可以看到,大部分的RDD操作由于缓存,节省了大量计算资源;
测试结果表明,在计算20阶,步长为5的矩阵运算时,优化前的计算时间为36.39秒;
优化后,将时间缩减到10.809秒,优化成果显著;
-
C语言之单位下三角矩阵求逆
2015-10-26 17:51:46//矩阵第一列求逆 for(i = 1; i ; i++ ) { l[i][0]=-l[i][0]; } //bs表示矩阵的维数 for(n=1;n;n++) { for(i = n+1; i ; i++) { for(j = 0; j ...#include<stdio.h>
int main(){
//矩阵保存在二位数组也可以随机生成
double l[4][4]={1,0,0,0,2,1,0,0,3,2,1,0,5,4,2,1};
long bs=4,i,j,n;
//矩阵第一列求逆
for(i = 1; i < bs; i++ )
{
l[i][0]=-l[i][0];
}
//bs表示矩阵的维数
for(n=1;n<bs-1;n++)
{
for(i = n+1; i < bs; i++)
{
for(j = 0; j < n; j++)
{
l[i][j]-=l[n][j]*l[i][n];
}
l[i][n]=-l[i][n];
}
}//求逆后保存在原矩阵
//打印矩阵
for(i=0;i<bs;i++)
{
for(j=0;j<bs;j++)
{
printf("%.2f ",l[i][j]);
}
printf("\n");
}
return 0;
} -
下三角矩阵的逆矩阵_上三角或下三角矩阵的逆矩阵能否简便方法求出??只有主副对角线不为0的矩阵能否直接写...
2020-12-18 12:48:59答2、下三角矩阵的逆矩阵将下三角矩阵划分成块矩阵,如上图所示,则其逆矩阵结果如下图。3、只有主对角线不为零的矩阵主对角元素取倒数,原位置不变。4、只有副对角线不为零的矩阵副对角元素取倒数,位置颠倒。示例...1、上三角矩阵的逆矩阵
将上三角矩阵划分成块矩阵,如上图所示,则其逆矩阵结果如下回图。答
2、下三角矩阵的逆矩阵
将下三角矩阵划分成块矩阵,如上图所示,则其逆矩阵结果如下图。
3、只有主对角线不为零的矩阵
主对角元素取倒数,原位置不变。
4、只有副对角线不为零的矩阵
副对角元素取倒数,位置颠倒。
示例如下:
扩展资料
矩阵求逆的求法
(1)初等变换法,通过初等变换将A矩阵变换成单位矩阵,则对应的单位矩阵变换成B矩阵,B矩阵即为A矩阵的逆矩阵,(A I)->(I B);
(2)伴随阵法,公式为:
;
(3)利用定义求逆矩阵
设A、B都是n阶方阵,如果存在n阶方阵B使得AB=BA=E,则称A为可逆矩阵,而称B为A的逆矩阵。
(4)恒等变形法
恒等变形法求逆矩阵的理论依据为逆矩阵的定义,此方法也常用与矩阵的理论推导上,就是通过恒等变形把要求的值化简出来,题目中的逆矩阵可以不求,利用
,把题目中的逆矩阵化简掉。
参考资料来源:百度百科--矩阵求逆
-
求下三角矩阵的逆矩阵的详细算法
2012-05-21 16:44:38矩阵计算中第一次实验题,计算下三角矩阵的逆矩阵的详细算法,可以正常运行,有所有的测试数据与运行结果 -
python矩阵求逆算法_09-30:Python矩阵求逆
2021-01-12 09:01:35A的逆等于U的逆乘于L的逆,L的逆就利用下三角矩阵求逆算法进行求解,U的逆可以这样求:先将U转置成下三角矩阵,再像对L求逆一样对U的转置求逆,再将得到的结果转置过来,得到的就是U的逆。因此,关键是下三角矩阵的...旁听了今天的上机课,收获良多。
方阵A求逆,先做LU分解。A的逆等于U的逆乘于L的逆,L的逆就利用下三角矩阵求逆算法进行求解,U的逆可以这样求:先将U转置成下三角矩阵,再像对L求逆一样对U的转置求逆,再将得到的结果转置过来,得到的就是U的逆。
因此,关键是下三角矩阵的求逆。
1.下三角矩阵求逆算法
我利用的公式计算公式如下:
对角元素.png
对角元素以下的元素.png
我的代码如下:
def triInverse(matA):
'''
@author:zengwei
输入:
matA:一个等待求逆的下三角矩阵,大小为n*n,并且希望里面的元素值是浮点数
输出:
matInv:matA的逆矩阵
'''
numRows = matA.shape[1]
matL = matA.copy()
matInv = np.zeros((numRows,numRows))
for row in np.arange(0,numRows):
matInv[row,row] = 1/matL[row,row]
for k in np.arange(row-1,-1,-1):
matInv[row,k] = -(np.dot(matInv[row,k+1:row+1],matL[k+1:row+1,k]))/matL[k,k]
return matInv
同样,生成一个矩阵来测试一下上面的函数,并得到输出。
import numpy as np
test_mat = np.array([[1,0,0,0],[2,3,0,0],[4,5,6,0],[7,8,9,10]],dtype=float)
'''
test_mat = array([[ 1., 0., 0., 0.],
[ 2., 3., 0., 0.],
[ 4., 5., 6., 0.],
[ 7., 8., 9., 10.]])
'''
Inv = triInverse(test_mat)
'''
Inv = array([[ 1. , 0. , 0. , 0. ],
[-0.66666667, 0.33333333, 0. , 0. ],
[-0.11111111, -0.27777778, 0.16666667, 0. ],
[-0.06666667, -0.01666667, -0.15 , 0.1 ]])
'''
显然,我不知道解出来对不对,但是我可以用这个逆矩阵Inv和测试矩阵test_mat相乘看看是不是单位矩阵来判断。
np.dot(test_mat,Inv)
'''
array([[1., 0., 0., 0.],
[0., 1., 0., 0.],
[0., 0., 1., 0.],
[0., 0., 0., 1.]])
'''
是单位矩阵,说明下三角矩阵求逆成功。接下来,利用上面的函数来进行矩阵的求逆。
2.矩阵求逆
首先,先贴出我LU分解的函数:
def getLandU(A):
'''
@author:zengwei
输入:
A:系数矩阵;array格式,且数值类型必须为浮点数,否则会出现截断误差。
输出:
LU分解中的L和U矩阵
'''
matA = A.copy()
if matA[0,0]==0:
print('不符合要求!需要换行操作。')
else:
numRows = matA.shape[0]
matL = np.identity(numRows) # 准备给L矩阵用
for row in np.arange(0,numRows-1): # 前n-1行,row代表第几行
for k in range(row+1,numRows): # 在第row行的情况下,遍历剩下的行
matL[k,row] = matA[k,row]/matA[row,row] # 获得L矩阵
if matA[k,row] != 0:
matA[k,:]=matA[k,:] - (matA[k,row]/matA[row,row])*matA[row,:]
return matL,matA # matL为L矩阵,matA变为U矩阵
然后我们利用getLandU函数和triInverse函数来写矩阵求解函数。如下:
def matInverse(A):
'''
@author:zengwei
输入:
A:想要求逆的系数矩阵,n*n,希望里面的数值是浮点数
输出:
A的逆矩阵
'''
L,U = getLandU(A) #LU分解
L_inv = triInverse(L) #L求逆
U_inv = triInverse(U.T).T #U求逆
return np.dot(U_inv,L_inv)
下面,我们生成一个随机矩阵来测试,并验证求得的逆矩阵和原矩阵相乘是否为单位矩阵。
TestMat = np.float64(np.random.randint(0,10,(4,4)))
'''
TestMat = array([[3., 7., 4., 0.],
[8., 3., 9., 3.],
[5., 4., 2., 4.],
[7., 0., 7., 6.]])
'''
TestMatInv = matInverse(TestMat)
isI = np.int64(np.dot(TestMatInv,TestMat)) # 验证
'''
isI = array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 0]], dtype=int64)
'''
可见,逆矩阵乘与原矩阵是单位矩阵,这里用了一下np.int64()函数,是为了让结果好看。否则原本是0的地方就会是非常非常非常小的数,这是由于计算机用浮点数运算得到的效果。
放假。
-
python矩阵求逆函数_09-30:Python矩阵求逆
2020-12-09 07:41:02A的逆等于U的逆乘于L的逆,L的逆就利用下三角矩阵求逆算法进行求解,U的逆可以这样求:先将U转置成下三角矩阵,再像对L求逆一样对U的转置求逆,再将得到的结果转置过来,得到的就是U的逆。因此,关键是下三角矩阵的... -
三角矩阵的逆矩阵怎么求_对称矩阵、对角矩阵与三角矩阵
2020-12-12 18:35:00例如: 可以看到,对称矩阵的转置等于其自身,即: 对角矩阵对角矩阵(Diagonal Matrix)是指除主对角线之外其他元素都为0的矩阵,例如: 三角矩阵三角矩阵(Triangular Matrix)分为上三角矩阵和下三角矩阵。... -
【线性代数】上三角矩阵/下三角矩阵
2018-07-20 09:50:20定义 性质 求逆 ...主对角线以下都是零的方阵称为上三角矩阵。...主对角线以上都是零的方阵称为下三角矩阵。...上(下)三角矩阵间的加减法和乘法运算的结果仍是上(下)三角矩阵。...一般化的上下三角矩阵求逆 ... -
Math01: 两种方法求下三角矩阵的逆
2015-09-18 00:13:00方法一:单位矩阵消元 \begin{equation}\left(A|I\right)\rightarrow\left(I|B\right)\end{equation} 1 clear; 2 n = 500; 3 A = zeros(n,n); 4 for j = 1:n 5 for i = j:n 6 A(i,j) = (i + j)^2; ..... -
矩阵求逆 LU三角分解
2018-09-02 19:18:42L是下三角矩阵,U是上三角矩阵,P是一个置换矩阵(P将在下一篇博客中写出) LUP分解:PA=LU② 由①②可得1.正向替换(设y=Ux):Ly=Pb 2.反向替换:Ux=y 所以 忽略P,下面说明LU的求法: 1.参数矩阵A做如下... -
c++矩阵求逆的lu分解实现.cpp
2020-01-08 12:36:49利用矩阵lu分解的优秀特性进行矩阵求逆的代码,减少求逆计算量,大约200行 求逆矩阵思路: 1.求矩阵的crout(LU)分解,其中L为下三角矩阵,U为上三角矩阵 2.求L,U矩阵的伴随阵,参考文献:三角形矩阵求伴随矩阵... -
加速LU分解和非奇异矩阵求逆
2020-04-06 10:40:04加速LU分解和非奇异矩阵求逆 1. 背景介绍 上(下)三角矩阵有许多性质使得计算机可以方便地对它们进行各种线性代数基础运算,例如,对三角矩阵的求逆非常直观容易。因此在计算中大型矩阵的逆时常将原矩阵AAA分解为一个... -
c++矩阵求逆的lu分解实现
2019-12-29 15:21:31初学c++,尝试利用基础知识编写矩阵求逆。但发现求解伴随阵的过程非常复杂且难以实现,而我正好看到一篇求三角阵伴随矩阵的文章,故尝试编程实现。在这种方法下,计算量明显减小,实现方法,思路适合初学者。 参考... -
【图像处理】矩阵运算代码实现2-矩阵求逆
2016-02-18 21:43:19这篇总结是《矩阵运算代码...将原矩阵A分解成两个三角矩阵,上三角矩阵L和下三角矩阵U。通过求解U和L对应的逆矩阵,即可求得相应A的逆矩阵。算法分成以下几个步骤: 1)矩阵的LU分解 对于LU上每个位置的值,可以用以 -
Matlab求矩阵的主对角元素、上(下)三角、矩阵的逆、行列式的值、矩阵的秩、矩阵的范数、矩阵的条件数和...
2020-03-24 16:03:55Matlab求矩阵的主对角元素、上(下)三角、矩阵的逆、行列式的值、矩阵的秩、矩阵的范数、矩阵的条件数和矩阵的迹 -
矩阵求逆运算 Python实现
2018-09-11 23:40:24编制下三角部分消元和上三角部分消元的代码 a. 从对角线元素往下比较取得这一列的最大值所在的行,与对角线元素所在行进行交换。 b. 从对角线所在行往下,利用矩阵的行变换将这一列下所有元素消为0,存储消元过程... -
LU 分解法求矩阵的逆
2010-01-05 17:09:49LU分解是矩阵分解的一种,可以将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积(有时是它们和一个置换矩阵的乘积)。LU分解主要应用在数值分析中,用来解线性方程、求矩阵的逆或计算行列式。 -
三阶矩阵的lu分解详细步骤_矩阵LU分解求逆(学习笔记)
2020-12-28 20:04:46矩阵的一种有效而广泛应用的分解方法是矩阵的LU三角分解,将一个n阶矩阵A分解为一...所以,矩阵求逆的算法流程可表述如下:1)进行LU分解;2)对分解后的L阵(下三角矩阵)和U阵(上三角矩阵)进行求逆;3)L阵的逆矩阵和U阵... -
分块矩阵乘法以及求逆应用
2020-03-08 18:07:39分块矩阵求逆,对角 分块矩阵乘 下三角或上三角用高斯消元求逆,相当比较快,一般的还要先化为下三角或者上三角。 用这个分块的方法计算量最复杂的地方就是要算3次2*2矩阵,好处是分块计算,总有对的... -
矩阵LU分解求逆(学习笔记)
2014-09-25 16:51:10矩阵的一种有效而广泛应用的分解方法是矩阵的LU三角分解,将一个n阶矩阵A分解为一个下三角矩阵L和一个上三角矩阵U的乘积。所以首先对矩阵进行三角分解,这里采用Doolittle分解...所以,矩阵求逆的算法流程可表述如下: -
求逆矩阵 —— LU分解法
2019-05-27 23:58:301. 将 A 矩阵分解为 L 下三角矩阵和 U 上三角矩阵。 step1.L对角线填充为1 step2. step3. step4. U是按行迭代计算,L是按列迭代计算,UL交错计算,且U先L一步 for k = 1 to m-1: { } 2. 分别对 L 和 U ... -
PLU-分解以及求逆矩阵
2021-02-03 11:17:40PLU-分解 PLU-分解是对LU分解的一种改进,其增加了选主元的操作...P为置换矩阵,L为下三角矩阵,U为上三角矩阵。选主元操作在计算过程中以一个一维数组保存代替n×nn\times nn×n的矩阵,以下为此算法的Fortran代码, -
利用dgetrf和dgetri进行LU分解求逆
2020-08-03 08:40:07下三角和上三角矩阵关于矩阵的求逆和乘积都具有封闭性,即下三角的矩阵的逆矩阵是下三角矩阵,两个下三角矩阵的乘积还是下三角矩阵;而上三角矩阵的逆矩阵是上三角矩阵,连个上三角矩阵的乘积还是上三角矩阵。 此... -
cholesky求逆
2020-10-21 15:16:36矩阵求逆经常会用到cholesky分解,接下来写一下主要的计算公式 设对称正定矩阵X=【】5*5矩阵 cholesky分解可得一个三角矩阵 如果是上三角矩阵,则计算公式为X=U'U,如果是下三角矩阵,则X=LL' 通过下面的例子可以... -
常见的几种矩阵分解方式
2016-09-25 15:54:231.三角分解(LU分解)矩阵的LU分解是将一个矩阵分解为一个下三角矩阵与上三角矩阵的乘积。本质上,LU分解是高斯消元的一种表达方式。首先,对矩阵A通过初等行变换将其变为一个上三角矩阵。对于学习过线性代数的同学来... -
矩阵分解算法
2018-06-14 21:01:45矩阵的三角分解体现了高斯消元的一个过程,将实矩阵拓展成A:E的形式,可以很方便的求解矩阵的逆,而如果我们将A变成一个下三角矩阵事,在A的所有顺序主子式不为0的情况下总是有E变为一上三角矩阵。不妨认为下三角... -
矩阵分解
2016-12-16 11:25:22矩阵的LU分解是将一个矩阵分解为一个下三角矩阵与上三角矩阵的乘积。本质上,LU分解是高斯消元的一种表达方式。首先,对矩阵A通过初等行变换将其变为一个上三角矩阵。对于学习过线性代数的同学来说,这个过程应该很... -
矩阵分解-LU分解
2020-10-09 15:06:09它的用途主要在简化一个大矩阵的行列式值的计算过程,求逆矩阵,和求解联立方程组。不过要注意这种分解法所得到的上下三角形矩阵并非唯一,还可找到数个不同 的一对上下三角形矩阵,此两三角形矩阵相乘也会得到原... -
矩阵分解-Cholesky分解
2020-10-09 14:32:02它的用途主要在简化一个大矩阵的行列式值的计算过程,求逆矩阵,和求解联立方程组。不过要注意这种分解法所得到的上下三角形矩阵并非唯一,还可找到数个不同 的一对上下三角形矩阵,此两三角形矩阵相乘也会得到原...
-
JAVA 常用类
-
机载激光雷达沙尘探测能量优化配置的统计研究
-
Qt调用libVLC实现播放器
-
精益开发治理的最佳实践,第3部分:角色和政策
-
斜光轴数字强度相关计量的像模糊容限
-
深究字符编码的奥秘,与乱码说再见
-
50个优秀的名片设计作品欣赏
-
PHP 日期 加减 月数,天数,周数,小时,分,秒等等
-
OpenCV-学习历程21- 霍夫变换2-圆检测(HoughCircles)
-
Amoeba 实现 MySQL 高可用、负载均衡和读写分离
-
网页元素轻设计–尊重用户产品体验
-
Python启蒙到架构师的核心技术精讲课程
-
精益用户体验(UX):摆脱只注重结果的工作
-
百瓦级全光纤线偏振激光振荡器
-
Spring实现事务管理
-
MySQL你该了解的那些事【服务端篇】
-
“Operation not permitted”报错
-
紫外区全角度光子晶体反射镜
-
系统设计:准备系统设计面试问题-源码
-
【网络通信与信息安全】之深入解析 TCP 的“拥塞控制”原理