Spark MLlib矩阵分解源码分析

基础知识

特征值分解

如果一个向量 vv 是方阵 AA 的特征向量，可以表示成下面的形式：

Av=λv

Av = \lambda v
其中， λ\lambda 为特征向量 vv 对应的特征值，矩阵 AA 的特征向量是相互正交的。
特征值分解是将矩阵 AA 分解为如下形式：

A=Q∑Q−1

A=Q\sum Q^{-1}
其中，矩阵 QQ 是 AA 的特征向量组成的矩阵， ∑\sum 是对角矩阵。

奇异值分解

如果矩阵 AA 不是方阵，是 m∗nm * n 的矩阵，m≥nm \geq n。奇异值分解是将矩阵 AA 分解成如下形式：

A=U∑VT

A=U\sum V^T
其中，UU 是 m∗mm*m 的方阵，里面的向量为左奇异向量，是相互正交的，VV 是 n∗nn*n 的方阵，里面的向量为右奇异向量，是相互正交的，∑\sum 是 m∗nm*n 的对角矩阵，对角线上的元素为奇异值。
关于SVD的详细解释和实际含义，请参考数据降维–SVD&CUR。

两者之间的关系

(ATA)v=λv

(A^TA)v = \lambda v
通常， ATAA^TA 是 AA 的列向量的格拉姆矩阵，AATAA^T 是 AA 的行向量的格拉姆矩阵。

矩阵分解细节推导，请参考矩阵分解。

求解方法

特征值分解和奇异值分解的计算通常有两种方法：
1、根据对角化矩阵计算所有特征值或者奇异值，计算复杂度为O(n3)O(n^3)；
2、使用迭代方法产生部分特征值或者奇异值，迭代方法主要用 power iteration、针对对称矩阵的 Lanczos iteration、针对非对称矩阵的 Arnoldi iteration。
对于上述的三种迭代方法，原理基本类似，以Lanczos方法为例。Lanczos方法具有 两端的特征值最先收敛的特点, 通过 Lanczos迭代得到一个三对角矩阵，三对角矩阵的特征值通过迭代，最大的特征值先收敛，其他的特征值则会聚集在最大的特征值周围, 继续迭代求解出其他的特征值和特征向量。更多细节，请参考 Lanczos方法：求稀疏矩阵特征值。
利用Lanczos方法计算特征值或者奇异值使用最广的算法包应属 ARPACK，ARPACK使用FORTRAN编写，通过 netlib-java和 breeze接口ARPACK可以在JVM上使用。不过有大神Yixuan Qiu用C++实现了类似于ARPACK功能的算法包 Spectra，用于大规模特征向量的计算，同时基于Spectra实现了R版的算法包 RSpectra ，使用起来非常方便。例如：

library(Matrix)
library(RSpectra)
n = 20
k = 5
set.seed(111)
A1 = matrix(rnorm(n^2), n)  ## class "matrix"
eigs(A1, k)
## Implicitly define the matrix by a function that calculates A %*% x
## 向量x为任意的向量，如果要计算乘积 A*x，name就需要知道A的所有信息，那么函数f就代表矩阵A本身
f = function(x, args)
{as.numeric(args %*% x)
}
eigs(f, k, n = n, args = A1)

ARPACK包

MLlib中SVD的实现

Spark MLlib在RowMatrix类中实现了SVD，使用方法如下：

val mat: RowMatrix = ...
// Compute the top 20 singular values and corresponding singular vectors.
val svd: SingularValueDecomposition[RowMatrix, Matrix] = mat.computeSVD(20, computeU = true)
val U: RowMatrix = svd.U // The U factor is a RowMatrix.
val s: Vector = svd.s // The singular values are stored in a local dense vector.
val V: Matrix = svd.V // The V factor is a local dense matrix.

RowMatrix#computeSVD 将矩阵A(m∗n)A(m *n)的形状分为三种：高瘦型（tall and skinny，m≫nm \gg n）、矮胖型（short and fat，m≪nm\ll n）以及近似方阵（square，m≈nm \approx n），根据矩阵形状的不同，采用不同的算法进行计算SVD。由于矮胖型矩阵通过矩阵转置变为高瘦型矩阵，下面主要讨论高瘦型（tall and skinny，m≫nm \gg n）矩阵和近似方阵（square，m≈nm \approx n）的分布式实现。

Square SVD with ARPACK

对于矩阵 AA 为近似方阵，采用ARPACK算法包计算格拉姆矩阵 ATAA^TA 的奇异值分解。在单机环境下，ARPACK算法包适合稀疏矩阵或者矩阵-向量乘积（matrix-vector product）形式的矩阵，只需要 O(n)O(n) 数量级的浮点数操作和存储。由于ARPACK具有处理任意格式矩阵格式的特性，这样就不用直接操作矩阵，可以通过矩阵-向量乘积预先定义的操作进行处理。调用ARPACK算法包时需要将输入矩阵变成矩阵-向量乘积（matrix-vector product）形式，这样矩阵操作从向量操作中分离出来。当ARPACK需要矩阵操作时，它向调用程序发送矩阵-向量乘积的请求，调用程序执行乘法操作并将结果向量(n∗1n*1)返回给ARPACK。通过使用Spark的分布式计算功能，就可以利用整个集群的计算资源实现分布式矩阵-向量乘法。因而，一方面利用了ARPACK的数值计算功能，另一方面利用了Spark的分布式计算能力。
为了使用ARPACK，需要计算 ATAvA^TAv：具体步骤如下：

上述步骤描述来自Stanford大学公开课： CME 323: Distributed Algorithms and Optimization。
重复上述步骤，直到有足够多的向量，ARPACK在单机上可以用来计算 ATAA^TA 中k个最大的特征值。

broadcast v, compute x = A %*% v
broadcast x, compute y = A^T %*% x
store y on driver

Tall and Skinny SVD

矩阵 AA 的奇异值和右奇异向量可以通过格拉姆矩阵 ATAA^TA中获得：

ATA=(U∑VT)TU∑VT=V∑UTU∑VT=V∑∑VT

A^TA = (U \sum V^T)^TU \sum V^T = V \sum U^TU\sum V^T = V\sum \sum V^T
左奇异向量通过如下方式获得：

U=AVsolve(∑)

U = AV solve(\sum)
其中， solve(∑)solve(\sum)为矩阵 ∑\sum 的逆矩阵。
对于高瘦型的矩阵 AA，当 m≫nm \gg n，格拉姆矩阵 ATAA^TA 比较小，在单个机器上计算即可。通常，计算矩阵 AA 秩为 rr 的SVD需要 O(mnr)O(mnr)，但是对于 ATAA^TA 只需要 O(n2r)O(n^2r) 次操作。
算法如下：

上述算法描述来自Stanford大学公开课： CME 323: Distributed Algorithms and Optimization。

MLlib SVD源码详解

在MLLib中，RowMatrix#computeSVD分为三步：

1、确定计算模式，即矩阵形状的确定
有三种模式：LocalARPACK, LocalLAPACK, DistARPACK。

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 n is small or k is large compared with n, we better compute the Gramian matrix first// and then compute its eigenvalues locally, instead of making multiple passes.if (k < n / 3) {SVDMode.LocalARPACK} else {SVDMode.LocalLAPACK}} else {// If k is small compared with n, we use ARPACK with distributed multiplication.SVDMode.DistARPACK}case "local-svd" => SVDMode.LocalLAPACKcase "local-eigs" => SVDMode.LocalARPACKcase "dist-eigs" => SVDMode.DistARPACKcase _ => throw new IllegalArgumentException(s"Do not support mode $mode.")
}

2、计算奇异值和右奇异向量
根据不同的矩阵形状，采用不同的算法求解。

val (sigmaSquares: BDV[Double], u: BDM[Double]) = computeMode match {case SVDMode.LocalARPACK =>require(k < n, s"k must be smaller than n in local-eigs mode but got k=$k and n=$n.")val G = computeGramianMatrix().toBreeze.asInstanceOf[BDM[Double]]// tol: termination tolerance EigenValueDecomposition.symmetricEigs(v => G * v, n, k, tol, maxIter)case SVDMode.LocalLAPACK =>// breeze (v0.10) svd latent constraint, 7 * n * n + 4 * n < Int.MaxValuerequire(n < 17515, s"$n exceeds the breeze svd capability")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.")}require(k < n, s"k must be smaller than n in dist-eigs mode but got k=$k and n=$n.")EigenValueDecomposition.symmetricEigs(multiplyGramianMatrixBy, n, k, tol, maxIter)
}

3、计算左奇异向量

if (computeU) {// N = Vk * Sk^{-1}val N = new BDM[Double](n, sk, Arrays.copyOfRange(u.data, 0, n * sk))var i = 0var j = 0while (j < sk) {i = 0val sigma = sigmas(j)while (i < n) {N(i, j) /= sigmai += 1}j += 1}val U = this.multiply(Matrices.fromBreeze(N))SingularValueDecomposition(U, s, V)
} else {SingularValueDecomposition(null, s, V)
}

格拉姆矩阵 ATAA^TA 计算

计算矩阵 AA 任意两列之间的相似性，如推荐中找相似的电影，本质上是计算格拉姆矩阵 ATAA^TA，它的元素 (i,j)(i ,j) 为矩阵 AA 中的列 cic_i 和 cjc_j 的点积 (cTi,cj)(c_i^T, c_j)。因而需要分布式计算 ATAA^TA。

普通方式

由于ATA=∑mi=0rirTiA^TA=\sum_{i=0}^m r_i r_i^T，rir_i 为矩阵 AA 的第i列，因此可以用如下MapReduce方式计算：

设矩阵 AA 的每一行至少有 LL 个非零元素，那么MapReduce shuffle大小为O(mL2)O(mL^2)，reduce-key最大为O(m)O(m)，而 mm 通常都比较大（10810^8），因此需要一种好的采样算法解决问题复杂性。

采样方式

DIMSUM: Dimension Independent Matrix Square Using MapReduce

DIMSUM Mapper与Naive Mapper很相似，只不过每个元素是以某概率的情形下发送出去而非全部元素发送，这样通过采样减少了计算代价。 DIMSUM Reduce汇总Mapper发送过来的数据，用Mapper的发送概率归一（scale）结果。发送概率有由可调参数 γ\gamma 控制：
- γ\gamma 较小时，保留 ATAA^TA 相似性；
- γ\gamma 较大时，保留 ATAA^TA 的奇异值；
shuffle大小变为O(nLγ)O(nL\gamma)，reduce-key最大为O(γ)O(\gamma)

参考资料：
1. CME 323: Distributed Algorithms and Optimization
2. Spectra
3. RSpectra
4. DIMSUM: Dimension Independent Matrix Square Using MapReduce