基于Spark的巨型矩阵分布式LU计算求逆【第一篇】

概述

本文将介绍如何利用Spark解决巨型矩阵分布式的LU法求逆的问题。本篇则将对LU求逆的前半部分——分布式LU分解做介绍。

【后半部分再我的队友整理完后，会以链接的形式补在这里】

 

正文

首先，我们看一下本地LU分解时有哪些特点，如下图

 

常规LU分解式为 A = L*U，而我们将A、L、U都按照相同的规格划分，不难看出这个矩阵等式可以得到

① A11 = L11*U11

② A21 = L21*U11

③ A12 = L11*U12

④ A22 = L21*U12 + L22*U22

发现式子①也是一个LU分解，而②、③可以结合式子①得出L21和U12。

那么现在就剩下一个恶心的式子④ 了，整理一下 ⑤ A22 - L21*U12 = L22*U22，其实也是个LU分解。

当一个矩阵很大时，经常采用分块的方式进行计算（比如大型矩阵乘法），所以LU也可以像上图这样分块。

而LU分块后的求解又包含了LU计算，那么迭代就可以解决了。

那么迭代哪里呢？

当然是最恶心的部分啦！

 

按照本地可算的矩阵大小，从左上角开始截取，那么式子①对A11的LU则本地可解[1]，

式子② -> L21 = A21*(U11^-1)

       ③ -> U12 = (L11^-1)*A12

对L11的求逆本地可解，对U11的逆我们采用 :上三角矩阵的转置的逆=他逆的转置

U11^-1 = ((U1^t)^-1)^t

所以他们同为下三角矩阵求逆[2] 。

之后的相乘，虽然A21、A12也很大但我们发现，因为你已经将整个矩阵分块，这两个乘法相当于

[1]

[2] * [5]    ， [5] * [1，2，3]   就是把A11这个单位块和A21、A12两个长条矩阵的每个单位块相乘。

[3]

这个Map一下就行了。

至于最后的A22 ,拿去迭代啦，不过要注意，你迭代的是A22 - L21*U12，

因为你下一次迭代是为了将剩下的大矩阵再分布式LU填充上一步，所以第一他满足LU，第二他的LU是现在缺的LU块。

 

激动人心的代码部分

本地矩阵的 LU 及 下三角求逆

def toLU(matrix:DenseMatrix): Tuple2[DenseMatrix,DenseMatrix] ={
    //初始化必要元素
    val N = matrix.numCols
    val A = new linalg.DenseMatrix[Double](N,N,matrix.toArray)
    if (N==1){
      return (new DenseMatrix(N,N,A.data),new DenseMatrix(N,N,A.data))
    }
    val LT = new linalg.DenseMatrix[Double](N,N)
    val UT = new linalg.DenseMatrix[Double](N,N)

    //初始化 LT对角线
    for(i <- 0.to(N-1,1))
      LT(i,i) = 1
    //求 LT UT
    for (i <- 0.to(N-1,1)){
      for(j <- i.to(N-1,1)){
        var temp = 0.0
        for(k <- 0.to(i-1,1))
          temp += LT(i,k) * UT(k,j)
        UT(i,j) = A(i,j) - temp //UT
      }
      for(j <- (i+1).to(N-1,1)){
        var temp2 = 0.0
        for(k <- 0.to(i-1,1))
          temp2 += LT(j,k) * UT(k,i)
        LT(j,i) = (A(j,i)-temp2) / UT(i,i)//LT
      }
    }
    (new DenseMatrix(N,N,LT.data),new DenseMatrix(N,N,UT.data))
  }

//本地下三角矩阵求逆
  def toINV(matrix:DenseMatrix): DenseMatrix = {
    //初始化必要元素
    val N = matrix.numCols
    val L = new linalg.DenseMatrix[Double](N,N,matrix.toArray);
    val Linv = new linalg.DenseMatrix[Double](N,N)

    for(i <- 0.to(N-1,1)){
      for (j <- 0.to(N-1,1)){
        if(i<j) Linv(i,j)=0
        else if(i==j) Linv(i,j) = 1.0 / L(i,i)
        else {
          var temp = 0.0
          for(k <- j.to(i-1,1))
            temp += L(i,k)*Linv(k,j)
          Linv(i, j) = (-1 / L(i, i)) * temp
        }
      }
    }
    new DenseMatrix(N,N,Linv.data)
  }

迭代部分

@annotation.tailrec
    def loop(Blocks_tmp:RDD[((Int,Int),DenseMatrix)],block_x:Int,block_y:Int,
             L:RDD[((Int,Int),DenseMatrix)],
             U:RDD[((Int,Int),DenseMatrix)]):
                          Tuple2[RDD[((Int,Int),DenseMatrix)],RDD[((Int,Int),DenseMatrix)]]={
      //获得矩阵 B11, B12, B21, B22
      val _B11 = Blocks_tmp.filter(blockId=>{
        (blockId._1._1==block_x)&&(blockId._1._2==block_y)})
        .map(_._2).collect()(0)
      val _B21  = Blocks_tmp.filter(blockId=>{(blockId._1._1>block_x)&&(blockId._1._2==block_y)})
      val _B12  = Blocks_tmp.filter(blockId=>{(blockId._1._1==block_x)&&(blockId._1._2>block_y)})
      val _B22  = Blocks_tmp.filter(blockId=>{(blockId._1._1>block_x)&&(blockId._1._2>block_y)})


      //求 B11的 LU -> L11, U11
      val B11_LU = LocalMatrixApi.toLU(_B11)
      val L1 = L.union(LocalMatrixApi.getRDD(B11_LU._1,block_x,block_y,spark))
      val U1 = U.union(LocalMatrixApi.getRDD(B11_LU._2,block_x,block_y,spark))

      //求 U12 = L11^-1 * B12
      val L11_inv = LocalMatrixApi.toINV(B11_LU._1)
      val U12 = _B12.map(m=>(m._1,L11_inv.multiply(m._2)))
      val U2 = U1.union(U12)

      //求 L21 = B21 * U11^-1
      val U11_inv = LocalMatrixApi.toINV(B11_LU._2.transpose).transpose
      val L21 = _B21.map(m=>(m._1,m._2.multiply(U11_inv)))
      val L2 = L1.union(L21)

      // 求下次迭代的 L22U22 --> B22-L21U12
      val L21_bm = new BlockMatrix(L21.map(block=>(block._1,Matrices.dense(block._2.numRows,block._2.numCols,block._2.toArray))),steps,steps)
      val U12_bm = new BlockMatrix(U12.map(block=>(block._1,Matrices.dense(block._2.numRows,block._2.numCols,block._2.toArray))),steps,steps)
      val B22_bm = new BlockMatrix(_B22.map(block=>(block._1,Matrices.dense(block._2.numRows,block._2.numCols,block._2.toArray))),steps,steps)
      val next_rdd = B22_bm.subtract(L21_bm.multiply(U12_bm))
        .blocks
        .filter(block=>{(block._1._1>block_x)&&(block._1._2>block_y)})
        .map(block=>(block._1,new DenseMatrix(block._2.numRows,block._2.numCols,block._2.toArray)))

      //跳出判断
      if(step_now == total-1){
        val B22 = next_rdd.map(_._2).collect()(0)
        val B22_LU = LocalMatrixApi.toLU(B22)
        val L3 = L2.union(LocalMatrixApi.getRDD(B22_LU._1,block_x+1,block_y+1,spark))
        val U3 = U2.union(LocalMatrixApi.getRDD(B22_LU._2,block_x+1,block_y+1,spark))
        (L3,U3)
      }
      else {step_now += 1;loop(next_rdd,block_x+1,block_y+1,L2,U2)}
    }
    //启动迭代
    val LU_final = loop(blocks,0,0,L_1,U_1)
    val L = new BlockMatrix(LU_final._1.map(block=>(block._1,Matrices.dense(block._2.numRows,block._2.numCols,block._2.toArray))),steps,steps)
    val U = new BlockMatrix(LU_final._2.map(block=>(block._1,Matrices.dense(block._2.numRows,block._2.numCols,block._2.toArray))),steps,steps)

总结

到这里，剩下的问题就是优化了，针对LU当前的可优化点如下：

1. 迭代中依然使用了大型分块矩阵BlockMatrix的运算——A22 - L21*U12。

    我们发现，如果将迭代过程倒过来，从右下角向左上角扩散则能去掉所有的BlockMatrix的运算，

    只有RDD的运算。【这个马上就能改出来，测试完就发】

２.将频繁使用的结果Persist到内存&硬盘，节省运算

3.减少Shuffle！！！【2、3我还在研究，研究完再上传吧】