基础知识
特征值分解
如果一个向量
其中,
特征值分解是将矩阵
其中,矩阵
奇异值分解
如果矩阵
其中,
关于SVD的详细解释和实际含义,请参考 数据降维–SVD&CUR。
两者之间的关系
通常,
矩阵分解细节推导,请参考矩阵分解。
求解方法
特征值分解和奇异值分解的计算通常有两种方法:
1、根据对角化矩阵计算所有特征值或者奇异值,计算复杂度为
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
将矩阵
Square SVD with ARPACK
对于矩阵
为了使用ARPACK,需要计算
上述步骤描述来自Stanford大学公开课: CME 323: Distributed Algorithms and Optimization。
重复上述步骤,直到有足够多的向量,ARPACK在单机上可以用来计算
broadcast v, compute x = A %*% v
broadcast x, compute y = A^T %*% x
store y on driver
Tall and Skinny SVD
矩阵
左奇异向量通过如下方式获得:
其中,
对于高瘦型的矩阵
算法如下:
上述算法描述来自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.LocalLAPACK
case "local-eigs" => SVDMode.LocalARPACK
case "dist-eigs" => SVDMode.DistARPACK
case _ => 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.MaxValue
require(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 = 0
var j = 0
while (j < sk) {
i = 0
val sigma = sigmas(j)
while (i < n) {
N(i, j) /= sigma
i += 1
}
j += 1
}
val U = this.multiply(Matrices.fromBreeze(N))
SingularValueDecomposition(U, s, V)
} else {
SingularValueDecomposition(null, s, V)
}
格拉姆矩阵
ATA
计算
计算矩阵
普通方式
由于
设矩阵
采样方式
DIMSUM: Dimension Independent Matrix Square Using MapReduce
DIMSUM Mapper与Naive Mapper很相似,只不过每个元素是以某概率的情形下发送出去而非全部元素发送,这样通过采样减少了计算代价。 DIMSUM Reduce汇总Mapper发送过来的数据,用Mapper的发送概率归一(scale)结果。发送概率有由可调参数
-
-
shuffle大小变为
参考资料:
1. CME 323: Distributed Algorithms and Optimization
2. Spectra
3. RSpectra
4. DIMSUM: Dimension Independent Matrix Square Using MapReduce