奇异值分解

1 奇异值分解

  在了解特征值分解之后，我们知道，矩阵A不一定是方阵。为了得到方阵，可以将矩阵A的转置乘以该矩阵。从而可以得到公式：

1.3

  现在假设存在M*N矩阵A，我们的目标是在n维空间中找一组正交基，使得经过A变换后还是正交的。假设已经找到这样一组正交基：

1.4

  A矩阵可以将这组正交基映射为如下的形式。

1.5

  要使上面的基也为正交基，即使它们两两正交，那么需要满足下面的条件。

1.6

  如果正交基v选择为$A^{T}A$的特征向量的话，由于$A^{T}A$是对称阵，v之间两两正交，那么

1.7

  由于下面的公式成立

1.8

  所以取单位向量

1.9

  可以得到(下面的公式有误，delta_i 应该等于sqrt(lamda_i))

1.10

  奇异值分解是一个能适用于任意的矩阵的一种分解的方法，它的形式如下：

1.11

  其中，U是一个M*M的方阵，它包含的向量是正交的，称为左奇异向量（即上文的u）。sigma是一个M*N的对角矩阵，每个对角线上的元素就是一个奇异值。V是一个N*N的矩阵，它包含的向量是正交的，称为右奇异向量（即上文的v）。

2 源码分析

  MLlib在RowMatrix类中实现了奇异值分解。下面是一个使用奇异值分解的例子。

  1. import org.apache.spark.mllib.linalg.Matrix
  2. import org.apache.spark.mllib.linalg.distributed.RowMatrix
  3. import org.apache.spark.mllib.linalg.SingularValueDecomposition
  4. val mat: RowMatrix = ...
  5. // Compute the top 20 singular values and corresponding singular vectors.
  6. val svd: SingularValueDecomposition[RowMatrix, Matrix] = mat.computeSVD(20, computeU = true)
  7. val U: RowMatrix = svd.U // The U factor is a RowMatrix.
  8. val s: Vector = svd.s // The singular values are stored in a local dense vector.
  9. val V: Matrix = svd.V // The V factor is a local dense matrix.

2.1 性能

  我们假设n比m小。奇异值和右奇异值向量可以通过方阵$A^{T}A$的特征值和特征向量得到。左奇异向量通过$AVS^{-1}$求得。
ml实际使用的方法方法依赖计算花费。

  • 当n很小（n<100）或者k比n大(k>n/2)，我们会首先计算方阵$A^{T}A$ ，然后在driver本地计算它的top特征值和特征向量。它的空间复杂度是O(n*n)，时间复杂度是O(n*n*k)。

  • 否则，我们用分布式的方式先计算$A^{T}Av$,然后把它传给ARPACK在driver上计算top特征值和特征向量。它需要传递O(k)的数据，每个executor的空间复杂度是O(n),driver的空间复杂度是O(nk)

2.2 代码实现

  1. def computeSVD(
  2.       k: Int,
  3.       computeU: Boolean = false,
  4.       rCond: Double = 1e-9): SingularValueDecomposition[RowMatrix, Matrix] = {
  5.     // 迭代次数
  6.     val maxIter = math.max(300, k * 3)
  7.     // 阈值
  8.     val tol = 1e-10
  9.     computeSVD(k, computeU, rCond, maxIter, tol, "auto")
  10. }

  computeSVD(k, computeU, rCond, maxIter, tol, "auto")的实现分为三步。分别是选择计算模式，$A^{T}A$的特征值分解，计算V,U,Sigma。
下面分别介绍这三步。

  • 1 选择计算模式
  1.  val computeMode = mode match {
  2.       case "auto" =>
  3.         if (k > 5000) {
  4.           logWarning(s"computing svd with k=$k and n=$n, please check necessity")
  5.         }
  6.         if (n < 100 || (k > n / 2 && n <= 15000)) {
  7.           // 满足上述条件，首先计算方阵，然后本地计算特征值，避免数据传递
  8.           if (k < n / 3) {
  9.             SVDMode.LocalARPACK
  10.           } else {
  11.             SVDMode.LocalLAPACK
  12.           }
  13.         } else {
  14.           // 分布式实现
  15.           SVDMode.DistARPACK
  16.         }
  17.       case "local-svd" => SVDMode.LocalLAPACK
  18.       case "local-eigs" => SVDMode.LocalARPACK
  19.       case "dist-eigs" => SVDMode.DistARPACK
  20.  }
  • 2 特征值分解
  1.  val (sigmaSquares: BDV[Double], u: BDM[Double]) = computeMode match {
  2.       case SVDMode.LocalARPACK =>
  3.         val G = computeGramianMatrix().toBreeze.asInstanceOf[BDM[Double]]
  4.         EigenValueDecomposition.symmetricEigs(v => G * v, n, k, tol, maxIter)
  5.       case SVDMode.LocalLAPACK =>
  6.         // breeze (v0.10) svd latent constraint, 7 * n * n + 4 * n < Int.MaxValue
  7.         val G = computeGramianMatrix().toBreeze.asInstanceOf[BDM[Double]]
  8.         val brzSvd.SVD(uFull: BDM[Double], sigmaSquaresFull: BDV[Double], _) = brzSvd(G)
  9.         (sigmaSquaresFull, uFull)
  10.       case SVDMode.DistARPACK =>
  11.         if (rows.getStorageLevel == StorageLevel.NONE) {
  12.           logWarning("The input data is not directly cached, which may hurt performance if its"
  13.             + " parent RDDs are also uncached.")
  14.         }
  15.         EigenValueDecomposition.symmetricEigs(multiplyGramianMatrixBy, n, k, tol, maxIter)
  16.     }

  当计算模式是SVDMode.LocalARPACK和SVDMode.LocalLAPACK时，程序实现的步骤是先获取方阵$A^{T}A$ ，在计算其特征值和特征向量。
获取方阵无需赘述，我们只需要注意它无法处理列大于65535的矩阵。我们分别看这两种模式下，如何获取特征值和特征向量。

  在SVDMode.LocalARPACK模式下，使用EigenValueDecomposition.symmetricEigs(v => G * v, n, k, tol, maxIter)计算特征值和特征向量。在SVDMode.LocalLAPACK模式下，直接使用breeze的方法计算。

  在SVDMode.DistARPACK模式下，不需要先计算方阵，但是传入EigenValueDecomposition.symmetricEigs方法的函数不同。

  1.  private[mllib] def multiplyGramianMatrixBy(v: BDV[Double]): BDV[Double] = {
  2.     val n = numCols().toInt
  3.     //v作为广播变量
  4.     val vbr = rows.context.broadcast(v)
  5.     rows.treeAggregate(BDV.zeros[Double](n))(
  6.       seqOp = (U, r) => {
  7.         val rBrz = r.toBreeze
  8.         val a = rBrz.dot(vbr.value)
  9.         rBrz match {
  10.           //计算y += x * a
  11.           case _: BDV[_] => brzAxpy(a, rBrz.asInstanceOf[BDV[Double]], U)
  12.           case _: BSV[_] => brzAxpy(a, rBrz.asInstanceOf[BSV[Double]], U)
  13.           case _ => throw new UnsupportedOperationException
  14.         }
  15.         U
  16.       }, combOp = (U1, U2) => U1 += U2)
  17.   }

  特征值分解的具体分析在特征值分解中有详细分析，请参考该文了解详情。

  • 3 计算U,V以及Sigma
  1.     //获取特征值向量
  2.     val sigmas: BDV[Double] = brzSqrt(sigmaSquares)
  3.     val sigma0 = sigmas(0)
  4.     val threshold = rCond * sigma0
  5.     var i = 0
  6.     // sigmas的长度可能会小于k
  7.     // 所以使用 i < min(k, sigmas.length) 代替 i < k.
  8.     if (sigmas.length < k) {
  9.       logWarning(s"Requested $k singular values but only found ${sigmas.length} converged.")
  10.     }
  11.     while (i < math.min(k, sigmas.length) && sigmas(i) >= threshold) {
  12.       i += 1
  13.     }
  14.     val sk = i
  15.     if (sk < k) {
  16.       logWarning(s"Requested $k singular values but only found $sk nonzeros.")
  17.     }
  18.     //计算s，也即sigma
  19.     val s = Vectors.dense(Arrays.copyOfRange(sigmas.data, 0, sk))
  20.     //计算V
  21.     val V = Matrices.dense(n, sk, Arrays.copyOfRange(u.data, 0, n * sk))
  22.     //计算U
  23.     // N = Vk * Sk^{-1}
  24.     val N = new BDM[Double](n, sk, Arrays.copyOfRange(u.data, 0, n * sk))
  25.     var i = 0
  26.     var j = 0
  27.     while (j < sk) {
  28.         i = 0
  29.         val sigma = sigmas(j)
  30.         while (i < n) {
  31.           //对角矩阵的逆即为倒数
  32.           N(i, j) /= sigma
  33.           i += 1
  34.         }
  35.         j += 1
  36.     }
  37.     //U=A * N
  38.     val U = this.multiply(Matrices.fromBreeze(N))

参考文献

【1】强大的矩阵奇异值分解(SVD)及其应用

【2】奇异值分解(SVD)原理详解及推导

【3】A Singularly Valuable Decomposition: The SVD of a Matrix