奇異值分解與主成分分析，一文帶你理解Spark分佈式降維方法

​原理 Spark是一個極為優秀的大數據框架，在大數據批處理上基本無人能敵，流處理上也有一席之地，機器學習則是當前正火熱AI人工智能的驅動引擎，在大數據場景下如何發揮AI技術成為優秀的大數據挖掘工程師必備技能。本文結合機器學習思想與Spark框架代碼結構來實現分佈式機器學習過程，希望與大家一起學習進步～

本文采用的組件版本為：Ubuntu 19.10、Jdk 1.8.0_241、Scala 2.11.12、Hadoop 3.2.1、Spark 2.4.5，老規矩先開啟一系列Hadoop、Spark服務與Spark-shell窗口：

降維是減少所考慮變量數量的過程。它可用於從原始和嘈雜的特徵中提取潛在特徵，或者在保持結構的同時壓縮數據。spark.mllib為RowMatrix類提供降維支持。

1.SVD介紹

奇異值分解（SVD）將矩陣分解為三個矩陣：U，Σ和V，使得：

這裡U是一個正交矩陣，其列稱為左奇異向量，Σ是對角矩陣，其中非對角線按降序排列，其對角線稱為奇異值，V是一個正交矩陣，其列稱為右奇異向量。

對於大型矩陣，通常不需要完整的因式分解，而僅需要頂部奇異值及其關聯的奇異矢量。這樣可以節省存儲空間，降低噪聲並恢復矩陣的低階結構。如果我們保留前k個奇異值，那麼所得的低秩矩陣的維將為：

• U：m*k
• Σ：k*k
• V：n*k

我們假設n小於m。奇異值和右奇異向量是從Gramian矩陣ATA的特徵值和特徵向量得出的。如果用戶通過computeU參數請求，則通過矩陣乘法將存儲左奇異矢量Ui的矩陣計算為U = A（VS-1）。實際使用的方法是根據計算成本自動確定的：w

• 如果n小（n <100）或k與n（k> n / 2）相比較大，我們首先計算Gramian矩陣，然後在驅動程序上局部計算其最高特徵值和特徵向量。這需要在每個執行器和驅動程序上進行一次O（n2）存儲操作，並在驅動程序上進行O（n2k）時間處理。
• 否則，我們將以分佈式方式計算（ATA）v並將其發送到ARPACK，以計算驅動程序節點上（ATA）的最高特徵值和特徵向量。這需要O（k）次通過，每個執行程序上的O（n）存儲以及驅動程序上的O（nk）存儲。

2.SVD實例

spark.mllib為RowMatrix類中提供的面向行的矩陣提供SVD功能。

<code>import org.apache.spark.mllib.linalg.Matrix
import org.apache.spark.mllib.linalg.SingularValueDecomposition
import org.apache.spark.mllib.linalg.Vector
import org.apache.spark.mllib.linalg.Vectors
import org.apache.spark.mllib.linalg.distributed.RowMatrix
// 定義數組
val data = Array(
  Vectors.sparse(5, Seq((1, 1.0), (3, 7.0))),
  Vectors.dense(2.0, 0.0, 3.0, 4.0, 5.0),
  Vectors.dense(4.0, 0.0, 0.0, 6.0, 7.0))
val rows = sc.parallelize(data)
val mat: RowMatrix = new RowMatrix(rows)
// 計算前5個奇異值和相應的奇異向量。
val svd: SingularValueDecomposition[RowMatrix, Matrix] = mat.computeSVD(5, computeU = true)
val U: RowMatrix = svd.U  // U因子
val s: Vector = svd.s     // 奇異值存儲在一個本地dense向量中
val V: Matrix = svd.V     // V因子/<code>

3.SVD源碼分析

計算SVD的源碼如下：

<code>def computeSVD(
      k: Int,
      computeU: Boolean = false,
      rCond: Double = 1e-9): SingularValueDecomposition[RowMatrix, Matrix] = {
    // 迭代次數
    val maxIter = math.max(300, k * 3)
    // 閾值
    val tol = 1e-10
    computeSVD(k, computeU, rCond, maxIter, tol, "auto")
}/<code>

computeSVD(k, computeU, rCond, maxIter, tol, "auto")的實現分為三步。分別是選擇計算模式，特徵值分解，計算V,U,Sigma。下面分別介紹這三步。首先是選擇計算模式：

<code> val computeMode = mode match {
      case "auto" =>
        if (k > 5000) {
          logWarning(s"computing svd with k=$k and n=$n, please check necessity")
        }
        if (n < 100 || (k > n / 2 && n <= 15000)) {
          // 滿足上述條件，首先計算方陣，然後本地計算特徵值，避免數據傳遞
          if (k < n / 3) {
            SVDMode.LocalARPACK
          } else {
            SVDMode.LocalLAPACK
          }
        } else {
          // 分佈式實現
          SVDMode.DistARPACK
        }
      case "local-svd" => SVDMode.LocalLAPACK
      case "local-eigs" => SVDMode.LocalARPACK
      case "dist-eigs" => SVDMode.DistARPACK
 }/<code>

特徵值分解：

<code> val (sigmaSquares: BDV[Double], u: BDM[Double]) = computeMode match {
      case SVDMode.LocalARPACK =>
        val G = computeGramianMatrix().toBreeze.asInstanceOf[BDM[Double]]
        EigenValueDecomposition.symmetricEigs(v => G * v, n, k, tol, maxIter)
      case SVDMode.LocalLAPACK =>
        // breeze (v0.10) svd latent constraint, 7 * n * n + 4 * n < Int.MaxValue
        val G = computeGramianMatrix().toBreeze.asInstanceOf[BDM[Double]]
        val brzSvd.SVD(uFull: BDM[Double], sigmaSquaresFull: BDV[Double], _) = brzSvd(G)
        (sigmaSquaresFull, uFull)
      case SVDMode.DistARPACK =>
        if (rows.getStorageLevel == StorageLevel.NONE) {
          logWarning("The input data is not directly cached, which may hurt performance if its"
            + " parent RDDs are also uncached.")
        }
        EigenValueDecomposition.symmetricEigs(multiplyGramianMatrixBy, n, k, tol, maxIter)
    }/<code>

計算U,V以及Sigma：

<code>//獲取特徵值向量
    val sigmas: BDV[Double] = brzSqrt(sigmaSquares)
    val sigma0 = sigmas(0)
    val threshold = rCond * sigma0
    var i = 0
    // sigmas的長度可能會小於k
    // 所以使用 i  
< min(k, sigmas.length) 代替 i < k.
    if (sigmas.length < k) {
      logWarning(s"Requested $k singular values but only found ${sigmas.length} converged.")
    }
    while (i < math.min(k, sigmas.length) && sigmas(i) >= threshold) {
      i += 1
    }
    val sk = i
    if (sk < k) {
      logWarning(s"Requested $k singular values but only found $sk nonzeros.")
    }
    //計算s，也即sigma
    val s = Vectors.dense(Arrays.copyOfRange(sigmas.data, 0, sk))
    //計算V
    val V = Matrices.dense(n, sk, Arrays.copyOfRange(u.data, 0, n * sk))
    //計算U
    // N = Vk * Sk^{-1}
    val N = new BDM[Double](n, sk, Arrays.copyOfRange(u.data, 0, n * sk))
    var i = 0
    var j = 0
    while (j < sk) {
        i = 0
        val sigma = sigmas(j)
        while (i < n) {
          //對角矩陣的逆即為倒數
          N(i, j) /= sigma
          i += 1
        }
        j += 1
    }
    //U=A * N
    val U = this.multiply(Matrices.fromBreeze(N))/<code>

4.PCA介紹

主成分分析是最常用的一種降維方法。我們首先考慮一個問題：對於正交矩陣空間中的樣本點，如何用一個超平面對所有樣本進行恰當的表達。容易想到，如果這樣的超平面存在，那麼他大概應該具有下面的性質。 基於最近重構性和最大可分性，能分別得到主成分分析的兩種等價推導。

• 最近重構性：樣本點到超平面的距離都足夠近
• 最大可分性：樣本點在這個超平面上的投影儘可能分開

主成分分析（PCA）是一種統計方法，用於查找旋轉，以使第一個座標具有最大的方差，而每個後續座標又具有最大的方差。旋轉矩陣的列稱為主成分。PCA被廣泛用於降維。spark.mllib支持將PCA用於以行格式和任何Vector存儲的高而瘦的矩陣。

5.PCA實例

以下代碼演示瞭如何在RowMatrix上計算主成分並將其用於將向量投影到低維空間中。

<code>import org.apache.spark.mllib.linalg.Matrix
import org.apache.spark.mllib.linalg.Vectors
import org.apache.spark.mllib.linalg.distributed.RowMatrix
val data = Array(
  Vectors.sparse(5, Seq((1, 1.0), (3, 7.0))),
  Vectors.dense(2.0, 0.0, 3.0, 4.0, 5.0),
  Vectors.dense(4.0, 0.0, 0.0, 6.0, 7.0))
val rows = sc.parallelize(data)
val mat: RowMatrix = new RowMatrix(rows)
// 計算4個主成分
// 主成分存儲在本地dense矩陣中
val pc: Matrix = mat.computePrincipalComponents(4)
// 將行投影到前4個主要成分所跨越的線性空間
val projected: RowMatrix = mat.multiply(pc)/<code>

6.PCA源碼分析

主成分分析的實現代碼在RowMatrix中實現。源碼如下：

<code>def computePrincipalComponents(k: Int): Matrix = {
    val n = numCols().toInt
    //計算協方差矩陣
    val Cov = computeCovariance().toBreeze.asInstanceOf[BDM[Double]]
    //特徵值分解
    val brzSvd.SVD(u: BDM[Double], _, _) = brzSvd(Cov)
    if (k == n) {
      Matrices.dense(n, k, u.data)
    } else {
      Matrices.dense(n, k, Arrays.copyOfRange(u.data, 0, n * k))
    }
  }/<code>

這段代碼首先會計算樣本的協方差矩陣，然後在通過breeze的svd方法進行奇異值分解。這裡由於協方差矩陣是方陣，所以奇異值分解等價於特徵值分解。下面是計算協方差的代碼：

<code>def computeCovariance(): Matrix = {
    val n = numCols().toInt
    checkNumColumns(n)
    val (m, mean) = rows.treeAggregate[(Long, BDV[Double])]((0L, BDV.zeros[Double](n)))(
      seqOp = (s: (Long, BDV[Double]), v: Vector) => (s._1 + 1L, s._2 += v.toBreeze),
      combOp = (s1: (Long, BDV[Double]), s2: (Long, BDV[Double])) =>
        (s1._1 + s2._1, s1._2 += s2._2)
    )
    updateNumRows(m)
    mean :/= m.toDouble
    // We use the formula Cov(X, Y) = E[X * Y] - E[X] E[Y], which is not accurate if E[X * Y] is
    // large but Cov(X, Y) is small, but it is good for sparse computation.
    // TODO: find a fast and stable way for sparse data.
    val G = computeGramianMatrix().toBreeze.asInstanceOf[BDM[Double]]
    var i = 0
    var j = 0
    val m1 = m - 1.0
    var alpha = 0.0
    while (i < n) {
      alpha = m / m1 * mean(i)
      j = i
      while (j < n) {
        val Gij = G(i, j) / m1 - alpha * mean(j)
        G(i, j) = Gij
        G(j, i) = Gij
        j += 1
      }
      i += 1
    }
    Matrices.fromBreeze(G)
  }/<code>

Spark 降維算法的內容至此結束，有關Spark的基礎文章可參考前文：

想要入門大數據？這篇文章不得不看！Spark源碼分析系列

阿里是怎麼做大數據的？淘寶怎麼能承載雙11？大數據之眸告訴你

Spark分佈式機器學習源碼分析：如何用分佈式集群構建線性模型？

高頻面經總結：最全大數據+AI方向面試100題（附答案詳解）

Spark分佈式機器學習系列：一文帶你理解並實戰樸素貝葉斯！

Spark分佈式機器學習系列：一文帶你理解並實戰決策樹模型！

Spark分佈式機器學習系列：一文帶你理解並實戰集成樹模型！

一文帶你理解並實戰協同過濾！Spark分佈式機器學習系列

Spark分佈式機器學習源碼分析：Kmeans族聚類

一文帶你理解並實戰Spark隱式狄利克雷分佈（LDA）

參考鏈接：

http://spark.apache.org/docs/latest/mllib-clustering.html

https://github.com/endymecy/spark-ml-source-analysis