Spark 矩阵matrix 原理用法示例源码详解CoordinateMatrix RowMatrix IndexedRowMatrix BlockMatrix
文章目录
CoordinateMatrix
Spark的CoordinateMatrix是用于存储和处理稀疏矩阵的数据结构。它通过以坐标形式存储非零元素的位置和值来表示稀疏矩阵，以支持高效的稀疏矩阵运算。
原理
CoordinateMatrix的数据存储方式是一个RDD，其中每个元素表示一个非零元素的位置和值。每个元素包括行索引、列索引和对应的值。在分布式环境中，可以并行地处理每个非零元素，从而高效地进行矩阵运算。 CoordinateMatrix通过存储非零元素的坐标，可以在大小和稀疏度不确定的情况下高效地表示矩阵。这在处理大规模稀疏矩阵时特别有用。
用法示例
下面是使用CoordinateMatrix的简单示例代码：
import org.apache.spark.sql.SparkSession
object CoordinateMatrixTest {
def main(args: Array[String]): Unit = {
import org.apache.spark.mllib.linalg.distributed.{CoordinateMatrix, MatrixEntry}
val spark = SparkSession.builder()
.appName("CoordinateMatrixTest")
.master("local")
.getOrCreate()
val sc = spark.sparkContext
// 创建一个RDD表示矩阵的非零元素
val entries = sc.parallelize(Seq(
MatrixEntry(0, 0, 1.0),
MatrixEntry(0, 1, 2.0),
MatrixEntry(1, 1, 3.0),
MatrixEntry(1, 2, 4.0),
MatrixEntry(2, 0, 5.0),
MatrixEntry(2, 2, 6.0)
))
// 将RDD转换为CoordinateMatrix
val coordinateMatrix = new CoordinateMatrix(entries)
// 计算转置矩阵
val transposedMatrix = coordinateMatrix.transpose()
// 显示转置矩阵的非零元素
transposedMatrix.entries.foreach(println)
}
}
//MatrixEntry(0,0,1.0)
//MatrixEntry(1,0,2.0)
//MatrixEntry(1,1,3.0)
//MatrixEntry(2,1,4.0)
//MatrixEntry(0,2,5.0)
//MatrixEntry(2,2,6.0)
在示例中，首先创建了一个RDD，表示一个稀疏矩阵的非零元素。然后，使用CoordinateMatrix的构造函数将RDD转换为CoordinateMatrix对象。接下来，使用transpose方法计算转置矩阵。最后，使用foreach方法遍历并打印转置矩阵的非零元素。 CoordinateMatrix提供了其他方法，如计算矩阵乘法、计算行数和列数等，可以根据具体需求进行使用。同样，需要注意在处理大规模矩阵时可能会受到内存限制，需要进行合理分区和计算资源管理。
中文源码
/**
* 表示分布式矩阵中的一个条目.
* @param i 行索引
* @param j 列索引
* @param value 条目的值
*/
@Since("1.0.0")
case class MatrixEntry(i: Long, j: Long, value: Double)
/**
* 以坐标形式表示矩阵.
*
* @param entries 矩阵条目
* @param nRows 行数. 非正数值表示未知, 此时行数将由最大行索引加1决定.
* @param nCols 列数. 非正数值表示未知, 此时列数将由最大列索引加1决定.
*/
@Since("1.0.0")
class CoordinateMatrix @Since("1.0.0") (
@Since("1.0.0") val entries: RDD[MatrixEntry],
private var nRows: Long,
private var nCols: Long) extends DistributedMatrix {
/** 另一种构造函数，将矩阵的维度自动确定. */
@Since("1.0.0")
def this(entries: RDD[MatrixEntry]) = this(entries, 0L, 0L)
/** 获取或计算列数. */
@Since("1.0.0")
override def numCols(): Long = {
if (nCols <= 0L) {
computeSize()
}
nCols
}
/** 获取或计算行数. */
@Since("1.0.0")
override def numRows(): Long = {
if (nRows <= 0L) {
computeSize()
}
nRows
}
/** 转置这个CoordinateMatrix. */
@Since("1.3.0")
def transpose(): CoordinateMatrix = {
new CoordinateMatrix(entries.map(x => MatrixEntry(x.j, x.i, x.value)), numCols(), numRows())
}
/**
* 转换为IndexedRowMatrix. 列数必须在整数范围内.
*/
@Since("1.0.0")
def toIndexedRowMatrix(): IndexedRowMatrix = {
val nl = numCols()
if (nl > Int.MaxValue) {
sys.error(s"无法转换为以行为导向的格式，因为列数 $nl 太大.")
}
val n = nl.toInt
val indexedRows = entries.map(entry => (entry.i, (entry.j.toInt, entry.value)))
.groupByKey()
.map { case (i, vectorEntries) =>
IndexedRow(i, Vectors.sparse(n, vectorEntries.toSeq))
}
new IndexedRowMatrix(indexedRows, numRows(), n)
}
/**
* 转换为RowMatrix, 在按行索引分组后丢弃行索引. 列数必须在整数范围内.
*/
@Since("1.0.0")
def toRowMatrix(): RowMatrix = {
toIndexedRowMatrix().toRowMatrix()
}
/**
* 转换为BlockMatrix. 创建大小为1024 x 1024的`SparseMatrix`块.
*/
@Since("1.3.0")
def toBlockMatrix(): BlockMatrix = {
toBlockMatrix(1024, 1024)
}
/**
* 转换为BlockMatrix. 创建`SparseMatrix`块.
* @param rowsPerBlock 每个块的行数. 底部边界的块可能具有较小的值. 必须为大于0的整数.
* @param colsPerBlock 每个块的列数. 右边边界的块可能具有较小的值. 必须为大于0的整数.
* @return 一个[[BlockMatrix]]
*/
@Since("1.3.0")
def toBlockMatrix(rowsPerBlock: Int, colsPerBlock: Int): BlockMatrix = {
require(rowsPerBlock > 0,
s"rowsPerBlock需要大于0. rowsPerBlock: $rowsPerBlock")
require(colsPerBlock > 0,
s"colsPerBlock需要大于0. colsPerBlock: $colsPerBlock")
val m = numRows()
val n = numCols()
// 由于块矩阵需要整数行和列索引
require(math.ceil(m.toDouble / rowsPerBlock) <= Int.MaxValue,
"行数除以rowsPerBlock不能超过最大整数.")
require(math.ceil(n.toDouble / colsPerBlock) <= Int.MaxValue,
"列数除以colsPerBlock不能超过最大整数.")
val numRowBlocks = math.ceil(m.toDouble / rowsPerBlock).toInt
val numColBlocks = math.ceil(n.toDouble / colsPerBlock).toInt
val partitioner = GridPartitioner(numRowBlocks, numColBlocks, entries.partitions.length)
val blocks: RDD[((Int, Int), Matrix)] = entries.map { entry =>
val blockRowIndex = (entry.i / rowsPerBlock).toInt
val blockColIndex = (entry.j / colsPerBlock).toInt
val rowId = entry.i % rowsPerBlock
val colId = entry.j % colsPerBlock
((blockRowIndex, blockColIndex), (rowId.toInt, colId.toInt, entry.value))
}.groupByKey(partitioner).map { case ((blockRowIndex, blockColIndex), entry) =>
val effRows = math.min(m - blockRowIndex.toLong * rowsPerBlock, rowsPerBlock).toInt
val effCols = math.min(n - blockColIndex.toLong * colsPerBlock, colsPerBlock).toInt
((blockRowIndex, blockColIndex), SparseMatrix.fromCOO(effRows, effCols, entry))
}
new BlockMatrix(blocks, rowsPerBlock, colsPerBlock, m, n)
}
/** 通过计算最大行/列索引确定大小. */
private def computeSize() {
// 如果 `entries` 为空则会抛出异常.
val (m1, n1) = entries.map(entry => (entry.i, entry.j)).reduce { case ((i1, j1), (i2, j2)) =>
(math.max(i1, i2), math.max(j1, j2))
}
// 最右边可能为空列, 最底边可能为空行.
nRows = math.max(nRows, m1 + 1L)
nCols = math.max(nCols, n1 + 1L)
}
/** 收集数据并组装成本地矩阵. */
private[mllib] override def toBreeze(): BDM[Double] = {
val m = numRows().toInt
val n = numCols().toInt
val mat = BDM.zeros[Double](m, n)
entries.collect().foreach { case MatrixEntry(i, j, value) =>
mat(i.toInt, j.toInt) = value
}
mat
}
}
RowMatrix
Spark的RowMatrix是用于存储和处理分布式稀疏矩阵的数据结构。它将稀疏矩阵表示为一组行向量，并以分布式方式存储和操作这些行向量。
原理
RowMatrix将稀疏矩阵分布在不同的分区上，并通过RDD进行处理。每个分区存储一个或多个行向量。在分布式环境中，可以并行地处理行向量，从而高效地进行矩阵运算。 RowMatrix的数据存储方式是一个键值对RDD，其中键是表示行的索引，值是对应的行向量。每个行向量都使用SparseVector或DenseVector表示。
用法示例
下面是使用RowMatrix的简单示例代码：
object RowMatrixTest {
def main(args: Array[String]): Unit = {
import org.apache.spark.mllib.linalg.distributed.{CoordinateMatrix, MatrixEntry}
val spark = SparkSession.builder()
.appName("RowMatrixTest")
.master("local")
.getOrCreate()
val sc = spark.sparkContext
import org.apache.spark.mllib.linalg.Matrix
import org.apache.spark.mllib.linalg.Vectors
import org.apache.spark.mllib.linalg.distributed.RowMatrix
// 创建一个RDD表示行向量
val rows = sc.parallelize(Seq(
Vectors.dense(1.0, 2.0, 3.0),
Vectors.dense(4.0, 5.0, 6.0),
Vectors.dense(7.0, 8.0, 9.0)
))
// 将RDD转换为RowMatrix
val rowMatrix = new RowMatrix(rows)
// 计算协方差矩阵
val covMatrix: Matrix = rowMatrix.computeCovariance()
// 显示协方差矩阵
println(covMatrix)
}
}
//9.0 9.0 9.0
//9.0 9.0 9.0
//9.0 9.0 9.0
在示例中，首先创建了一个RDD，表示一个3x3的稠密矩阵。然后，使用RowMatrix的构造函数将RDD转换为RowMatrix对象。接下来，使用computeCovariance方法计算协方差矩阵。最后，使用println方法显示协方差矩阵。 RowMatrix提供了其他方法，如计算特征向量和特征值等，可以根据具体需求进行使用。值得注意的是，RowMatrix在处理大规模矩阵时可能会受到内存限制，需要进行合理分区和计算资源管理。
中文源码
/**
* 表示一个行向量分布矩阵，没有有意义的行索引。
*
* @param rows 以RDD[Vector]存储的行
* @param nRows 行数。非正值表示未知，此时行数将根据RDD `rows` 中记录的数量来确定。
* @param nCols 列数。非正值表示未知，此时列数将根据第一行的大小来确定。
*/
@Since("1.0.0")
class RowMatrix @Since("1.0.0") (
@Since("1.0.0") val rows: RDD[Vector],
private var nRows: Long,
private var nCols: Int) extends DistributedMatrix with Logging {
/** 使用自动确定矩阵维度的替代构造函数。 */
@Since("1.0.0")
def this(rows: RDD[Vector]) = this(rows, 0L, 0)
/** 获取或计算列数。 */
@Since("1.0.0")
override def numCols(): Long = {
if (nCols <= 0) {
try {
// 如果`rows`为空，则调用`first`会抛出异常。
nCols = rows.first().size
} catch {
case err: UnsupportedOperationException =>
sys.error("无法确定列数，因为未在构造函数中指定且行RDD为空。")
}
}
nCols
}
/** 获取或计算行数。 */
@Since("1.0.0")
override def numRows(): Long = {
if (nRows <= 0L) {
nRows = rows.count()
if (nRows == 0L) {
sys.error("无法确定行数，因为未在构造函数中指定且行RDD为空。")
}
}
nRows
}
/**
* 将Gramian矩阵`A^T A`乘以右侧的密集向量，而不计算`A^T A`。
*
* @param v 一个密集向量，其长度必须与该矩阵的列数匹配
* @return 表示乘积的密集向量
*/
private[mllib] def multiplyGramianMatrixBy(v: BDV[Double]): BDV[Double] = {
val n = numCols().toInt
val vbr = rows.context.broadcast(v)
rows.treeAggregate(BDV.zeros[Double](n))(
seqOp = (U, r) => {
val rBrz = r.asBreeze
val a = rBrz.dot(vbr.value)
rBrz match {
// 使用特殊的axpy方法以获得更好的性能
case _: BDV[_] => brzAxpy(a, rBrz.asInstanceOf[BDV[Double]], U)
case _: BSV[_] => brzAxpy(a, rBrz.asInstanceOf[BSV[Double]], U)
case _ => throw new UnsupportedOperationException(
s"不支持来自类型${rBrz.getClass.getName}的向量操作。")
}
U
}, combOp = (U1, U2) => U1 += U2)
}
/**
* 计算Gramian矩阵`A^T A`。
*
* @note 无法计算超过65535列的矩阵。
*/
@Since("1.0.0")
def computeGramianMatrix(): Matrix = {
val n = numCols().toInt
checkNumColumns(n)
// 计算n*(n+1)/2，避免乘法溢出。
// 这个在n<=65535时成功，以上已经检查过
val nt = if (n % 2 == 0) ((n / 2) * (n + 1)) else (n * ((n + 1) / 2))
// 计算gram矩阵的上三角部分。
val GU = rows.treeAggregate(null.asInstanceOf[BDV[Double]])(
seqOp = (maybeU, v) => {
val U =
if (maybeU == null) {
new BDV[Double](nt)
} else {
maybeU
}
BLAS.spr(1.0, v, U.data)
U
}, combOp = (U1, U2) =>
if (U1 == null) {
U2
} else if (U2 == null) {
U1
} else {
U1 += U2
}
)
RowMatrix.triuToFull(n, GU.data)
}
private def checkNumColumns(cols: Int): Unit = {
if (cols > 65535) {
throw new IllegalArgumentException(s"超过65535列的参数: $cols")
}
if (cols > 10000) {
val memMB = (cols.toLong * cols) / 125000
logWarning(s"$cols 列将需要至少 $memMB 兆字节的内存！")
}
}
}
/**
* 计算该矩阵的奇异值分解。假设该矩阵为A（m x n）。这个方法会计算矩阵U，S，V，使得A ~= U * S * V'，其中S包含前k个奇异值，U和V包含相应的奇异向量。
*
* 最多返回k个非零奇异值及其对应的向量。如果存在k个这样的值，则返回的维度为：
* - U是一个大小为m x k的RowMatrix，满足U' * U = eye(k),
* - s是一个大小为k的向量，按降序排列，存储奇异值，
* - V是一个大小为n x k的矩阵，满足V' * V = eye(k)。
*
* 我们假设n小于m，尽管这并不是严格要求。
* 奇异值和右奇异向量是通过Gram矩阵A' * A的特征值和特征向量导出的。如果用户请求，右奇异向量矩阵U通过矩阵乘法计算为U = A * (V * S^-1^)。实际使用的方法是根据成本自动确定的：
* - 如果n较小（n < 100）或者k与n相比较大（k > n / 2），我们先计算Gram矩阵，然后在驱动程序上本地计算其前k个特征值和特征向量。这需要每个执行器和驱动程序上的一次通过，以及驱动程序上的O(n^2^ k)时间和O(n^2^)存储。
* - 否则，我们以分布式方式计算(A' * A) * v，并将其发送给ARPACK的DSAUPD来在驱动程序节点上计算(A' * A)的前k个特征值和特征向量。这需要O(k)次迭代，每个执行器上的O(n)存储，以及驱动程序上的O(n k)存储。
*
* 一些内部参数已设置为默认值。倒数条件数rCond设置为1e-9。所有小于rCond * sigma(0)的奇异值都被视为零，其中sigma(0)是最大的奇异值。ARPACK的最大Arnoldi更新迭代次数设置为300或k * 3，取其中较大的值。ARPACK的特征值分解的数值容差设置为1e-10。
*
* @param k 保留的前k个奇异值的数目（0 < k <= n）。
* 如果存在数值上为零的奇异值，或者在达到最大Arnoldi更新迭代次数之前没有足够的Ritz值收敛（在矩阵A条件较差的情况下），可能会返回少于k个奇异值。
* @param computeU 是否计算U
* @param rCond 倒数条件数。所有小于rCond * sigma(0)的奇异值都被视为零，其中sigma(0)是最大的奇异值。
* @return SingularValueDecomposition(U, s, V)。如果computeU = false，则U = null。
*
* @note 决定内部使用哪种方法和默认参数的条件可能会发生变化。
*/
@Since("1.0.0")
def computeSVD(
k: Int,
computeU: Boolean = false,
rCond: Double = 1e-9): SingularValueDecomposition[RowMatrix, Matrix] = {
// ARPACK的最大Arnoldi更新迭代次数
val maxIter = math.max(300, k * 3)
// ARPACK的数值容差
val tol = 1e-10
computeSVD(k, computeU, rCond, maxIter, tol, "auto")
}
/**
* SVD算法的实际实现，用于测试。
*
* @param k 保留的前k个奇异值（0 < k <= n）
* @param computeU 是否计算U
* @param rCond 倒数条件数
* @param maxIter 最大迭代次数（如果使用ARPACK）
* @param tol 终止容差（如果使用ARPACK）
* @param mode 计算模式（auto：自动确定使用的模式，
* local-svd：计算格拉姆矩阵并在本地计算其完整的SVD，
* local-eigs：计算格拉姆矩阵并在本地计算其前k个特征值，
* dist-eigs：分布式计算格拉姆矩阵的前k个特征值）
* @return SingularValueDecomposition(U, s, V)。如果computeU=false，则U=null。
*/
private[mllib] def computeSVD(
k: Int,
computeU: Boolean,
rCond: Double,
maxIter: Int,
tol: Double,
mode: String): SingularValueDecomposition[RowMatrix, Matrix] = {
val n = numCols().toInt
require(k > 0 && k <= n, s"要求k个奇异值，但是传入的k=$k numCols=$n。")
object SVDMode extends Enumeration {
val LocalARPACK, LocalLAPACK, DistARPACK = Value
}
val computeMode = mode match {
case "auto" =>
if (k > 5000) {
logWarning(s"计算k=$k和n=$n的SVD，请检查是否必要")
}
// TODO：以下条件尚未完全测试。
if (n < 100 || (k > n / 2 && n <= 15000)) {
// 如果n很小或者k相对于n很大，最好先计算格拉姆矩阵，
// 然后在本地计算其特征值，而不是多次遍历。
if (k < n / 3) {
SVDMode.LocalARPACK
} else {
SVDMode.LocalLAPACK
}
} else {
// 如果k相对于n很小，我们使用分布式乘法的ARPACK算法。
SVDMode.DistARPACK
}
case "local-svd" => SVDMode.LocalLAPACK
case "local-eigs" => SVDMode.LocalARPACK
case "dist-eigs" => SVDMode.DistARPACK
case _ => throw new IllegalArgumentException(s"不支持的模式 $mode。")
}
// 计算A' * A的特征分解。
val (sigmaSquares: BDV[Double], u: BDM[Double]) = computeMode match {
case SVDMode.LocalARPACK =>
require(k < n, s"在local-eigs模式下，k必须小于n，但是传入的k=$k n=$n。")
val G = computeGramianMatrix().asBreeze.asInstanceOf[BDM[Double]]
EigenValueDecomposition.symmetricEigs(v => G * v, n, k, tol, maxIter)
case SVDMode.LocalLAPACK =>
// breeze (v0.10) svd 潜在约束，7 * n * n + 4 * n < Int.MaxValue
require(n < 17515, s"$n 超过breeze的svd能力")
val G = computeGramianMatrix().asBreeze.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("输入数据没有直接缓存，如果其父RDD也没有缓存，那么性能可能受到影响。")
}
require(k < n, s"在dist-eigs模式下，k必须小于n，但是传入的k=$k n=$n。")
EigenValueDecomposition.symmetricEigs(multiplyGramianMatrixBy, n, k, tol, maxIter)
}
val sigmas: BDV[Double] = brzSqrt(sigmaSquares)
// 确定有效的秩。
val sigma0 = sigmas(0)
val threshold = rCond * sigma0
var i = 0
// 如果某些里兹值不满足tol指定的收敛标准，则sigmas的长度可能小于k。
// 因此使用i < min(k, sigmas.length)而不是i < k。
if (sigmas.length < k) {
logWarning(s"要求$k个奇异值，但是只找到${sigmas.length}个收敛的奇异值。")
}
while (i < math.min(k, sigmas.length) && sigmas(i) >= threshold) {
i += 1
}
val sk = i
if (sk < k) {
logWarning(s"要求$k个奇异值，但是只找到$sk个非零奇异值。")
}
// 末尾也要发出警告，以增加可见性。
if (computeMode == SVDMode.DistARPACK && rows.getStorageLevel == StorageLevel.NONE) {
logWarning("输入数据没有直接缓存，如果其父RDD也没有缓存，那么性能可能受到影响。")
}
val s = Vectors.dense(Arrays.copyOfRange(sigmas.data, 0, sk))
val V = Matrices.dense(n, sk, Arrays.copyOfRange(u.data, 0, n * sk))
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)
}
}
/**
* 计算协方差矩阵，将每一行都视为一个观察值。
*
* @return 一个大小为 n x n 的本地密集矩阵
*
* @note 无法对列数大于65535的矩阵进行计算
*/
@Since("1.0.0")
def computeCovariance(): Matrix = {
val n = numCols().toInt
checkNumColumns(n)
val summary = computeColumnSummaryStatistics()
val m = summary.count
require(m > 1, s"RowMatrix.computeCovariance 要求矩阵必须至少有 $m 行." +
" 无法对行数 <= 1 的 RowMatrix 计算协方差.")
val mean = summary.mean
// 我们使用 Cov(X, Y) = E[X * Y] - E[X] E[Y] 的公式来计算协方差，这个公式对于 E[X * Y] 很大但 Cov(X, Y) 很小的情况不准确，
// 但对于稀疏计算是好的。
// TODO: 找到一种对稀疏数据适用的快速且稳定的方法。
val G = computeGramianMatrix().asBreeze
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)
}
/**
* 计算前 k 个主成分以及每个主成分解释的方差比例的向量。
* 行对应观察，列对应变量。
* 主成分以 n x k 的本地矩阵存储。
* 每一列对应一个主成分，列按照分量的方差大小降序排列。
* 行数据在计算之前不需要“居中”，每一列的均值为 0 是不必要的。
*
* @param k 前 k 个主成分数
* @return 一个大小为 n x k 的矩阵，其中的列为主成分，以及一个值向量，指示每个主成分解释的方差比例
*
* @note 无法对列数大于65535的矩阵进行计算
*/
@Since("1.6.0")
def computePrincipalComponentsAndExplainedVariance(k: Int): (Matrix, Vector) = {
val n = numCols().toInt
require(k > 0 && k <= n, s"k = $k 超出了范围 (0, n = $n]")
val Cov = computeCovariance().asBreeze.asInstanceOf[BDM[Double]]
val brzSvd.SVD(u: BDM[Double], s: BDV[Double], _) = brzSvd(Cov)
val eigenSum = s.data.sum
val explainedVariance = s.data.map(_ / eigenSum)
if (k == n) {
(Matrices.dense(n, k, u.data), Vectors.dense(explainedVariance))
} else {
(Matrices.dense(n, k, Arrays.copyOfRange(u.data, 0, n * k)),
Vectors.dense(Arrays.copyOfRange(explainedVariance, 0, k)))
}
}
/**
* 仅计算前 k 个主成分。
*
* @param k 前 k 个主成分数
* @return 一个大小为 n x k 的矩阵，其中的列为主成分
* @see computePrincipalComponentsAndExplainedVariance
*/
@Since("1.0.0")
def computePrincipalComponents(k: Int): Matrix = {
computePrincipalComponentsAndExplainedVariance(k)._1
}
/**
* 计算列统计摘要。
*/
@Since("1.0.0")
def computeColumnSummaryStatistics(): MultivariateStatisticalSummary = {
val summary = rows.treeAggregate(new MultivariateOnlineSummarizer)(
(aggregator, data) => aggregator.add(data),
(aggregator1, aggregator2) => aggregator1.merge(aggregator2))
updateNumRows(summary.count)
summary
}
/**
* 将该矩阵与右边的一个本地矩阵进行乘法运算。
*
* @param B 一个本地矩阵，其行数必须与该矩阵的列数相匹配
* @return 一个 [[org.apache.spark.mllib.linalg.distributed.RowMatrix]]，表示乘积，保留了分区
*/
@Since("1.0.0")
def multiply(B: Matrix): RowMatrix = {
val n = numCols().toInt
val k = B.numCols
require(n == B.numRows, s"维度不匹配：$n vs ${B.numRows}")
require(B.isInstanceOf[DenseMatrix],
s"目前只支持密集矩阵，但发现了 ${B.getClass.getName}.")
val Bb = rows.context.broadcast(B.asBreeze.asInstanceOf[BDM[Double]].toDenseVector.toArray)
val AB = rows.mapPartitions { iter =>
val Bi = Bb.value
iter.map { row =>
val v = BDV.zeros[Double](k)
var i = 0
while (i < k) {
v(i) = row.asBreeze.dot(new BDV(Bi, i * n, 1, n))
i += 1
}
Vectors.fromBreeze(v)
}
}
new RowMatrix(AB, nRows, B.numCols)
}
/**
* 使用计算标准化点积的暴力方法计算该矩阵的列之间的所有余弦相似度。
*
* @return 一个大小为 n x n 的上三角稀疏矩阵，其中的元素为列之间的余弦相似度
*/
@Since("1.2.0")
def columnSimilarities(): CoordinateMatrix = {
columnSimilarities(0.0)
}
/**
* 使用采样方法计算矩阵的列之间的相似度。
*
* threshold 参数是一个在估计质量和计算成本之间的折中参数。
*
* 设置阈值为0可以保证确定性的正确结果，但与蛮力方法的计算成本完全相同。
* 将阈值设置为正值比蛮力方法产生更少的计算成本，但计算出的相似度将是估计值。
*
* 采样对于那些相似度大于给定相似度阈值的列对保证相对误差的正确性。
*
* 为了描述这个保证，我们设置了一些符号:
* 让 A 是这个矩阵中绝对值非零元素的最小值。
* 让 B 是这个矩阵中绝对值非零元素的最大值。
* 让 L 是每行的非零元素的最大数量。
*
* 例如，对于 {0,1} 矩阵: A=B=1。
* 另一个例子，对于 Netflix 矩阵: A=1，B=5
*
* 对于那些高于阈值的列对，计算的相似度正确到相对误差不超过 20% 的概率为
* 最少为 1 - (0.981)^10/B^
*
* 洗牌大小受以下两个表达式中较小的一个限制:
*
* O(n log(n) L / (threshold * A))
* O(m L^2^)
*
* 后者是蛮力方法的计算成本，所以对于非零阈值，计算成本总是比蛮力方法便宜。
*
* @param threshold 设置为0以获得确定性的正确性保证。
* 高于此阈值的相似度是估计的，
* 具体取决于上述折中关系的计算成本和估计质量。
* @return 一个 n x n 的稀疏上三角矩阵，表示该矩阵的列之间的余弦相似度。
*/
@Since("1.2.0")
def columnSimilarities(threshold: Double): CoordinateMatrix = {
require(threshold >= 0, s"Threshold cannot be negative: $threshold")
if (threshold > 1) {
logWarning(s"Threshold is greater than 1: $threshold " +
"Computation will be more efficient with promoted sparsity, " +
"however there is no correctness guarantee.")
}
val gamma = if (threshold < 1e-6) {
Double.PositiveInfinity
} else {
10 * math.log(numCols()) / threshold
}
columnSimilaritiesDIMSUM(computeColumnSummaryStatistics().normL2.toArray, gamma)
}
/**
* 计算[[RowMatrix]]的QR分解。该实现旨在优化[[RowMatrix]]的QR分解（因子分解），适用于“高而瘦”的矩阵。
* 参考文献:
* Paul G. Constantine, David F. Gleich. "Tall and skinny QR factorizations in MapReduce
* architectures" (详见 <a href="http://dx.doi.org/10.1145/1996092.1996103">这里</a>)
*
* @param computeQ 是否计算 Q
* @return QRDecomposition(Q, R)，如果 computeQ = false，则 Q = null。
*/
@Since("1.5.0")
def tallSkinnyQR(computeQ: Boolean = false): QRDecomposition[RowMatrix, Matrix] = {
val col = numCols().toInt
// 将行水平分为较小的矩阵，并对每个矩阵计算 QR 分解
val blockQRs = rows.retag(classOf[Vector]).glom().filter(_.length != 0).map { partRows =>
val bdm = BDM.zeros[Double](partRows.length, col)
var i = 0
partRows.foreach { row =>
bdm(i, ::) := row.asBreeze.t
i += 1
}
breeze.linalg.qr.reduced(bdm).r
}
// 将前面的结果中的 R 部分在垂直方向组合成一个高矩阵
val combinedR = blockQRs.treeReduce { (r1, r2) =>
val stackedR = BDM.vertcat(r1, r2)
breeze.linalg.qr.reduced(stackedR).r
}
val finalR = Matrices.fromBreeze(combinedR.toDenseMatrix)
val finalQ = if (computeQ) {
try {
val invR = inv(combinedR)
this.multiply(Matrices.fromBreeze(invR))
} catch {
case err: MatrixSingularException =>
logWarning("R is not invertible and return Q as null")
null
}
} else {
null
}
QRDecomposition(finalQ, finalR)
}
/**
* 使用DIMSUM抽样算法查找所有相似的列，算法在以下两篇论文中有描述：
*
* http://arxiv.org/abs/1206.2082
* http://arxiv.org/abs/1304.1467
*
* @param colMags 列幅值的向量
* @param gamma 过采样参数。为了得到可证明的结果，设置为 10 * log(n) / s，
* 其中 s 是要估算的最小相似度分值，n 是列的个数
* @return 一个 n x n 的稀疏上三角矩阵，表示该矩阵的列之间的余弦相似度
*/
private[mllib] def columnSimilaritiesDIMSUM(
colMags: Array[Double],
gamma: Double): CoordinateMatrix = {
require(gamma > 1.0, s"过采样参数必须大于 1: $gamma")
require(colMags.size == this.numCols(), "幅值的个数与列的维度不匹配")
val sg = math.sqrt(gamma) // sqrt(gamma) 多次使用
// 对于幅值为 0 的列，不要除以 0
val colMagsCorrected = colMags.map(x => if (x == 0) 1.0 else x)
val sc = rows.context
val pBV = sc.broadcast(colMagsCorrected.map(c => sg / c))
val qBV = sc.broadcast(colMagsCorrected.map(c => math.min(sg, c)))
val sims = rows.mapPartitionsWithIndex { (indx, iter) =>
val p = pBV.value
val q = qBV.value
val rand = new XORShiftRandom(indx)
val scaled = new Array[Double](p.size)
iter.flatMap { row =>
row match {
case SparseVector(size, indices, values) =>
val nnz = indices.size
var k = 0
while (k < nnz) {
scaled(k) = values(k) / q(indices(k))
k += 1
}
Iterator.tabulate (nnz) { k =>
val buf = new ListBuffer[((Int, Int), Double)]()
val i = indices(k)
val iVal = scaled(k)
if (iVal != 0 && rand.nextDouble() < p(i)) {
var l = k + 1
while (l < nnz) {
val j = indices(l)
val jVal = scaled(l)
if (jVal != 0 && rand.nextDouble() < p(j)) {
buf += (((i, j), iVal * jVal))
}
l += 1
}
}
buf
}.flatten
case DenseVector(values) =>
val n = values.size
var i = 0
while (i < n) {
scaled(i) = values(i) / q(i)
i += 1
}
Iterator.tabulate (n) { i =>
val buf = new ListBuffer[((Int, Int), Double)]()
val iVal = scaled(i)
if (iVal != 0 && rand.nextDouble() < p(i)) {
var j = i + 1
while (j < n) {
val jVal = scaled(j)
if (jVal != 0 && rand.nextDouble() < p(j)) {
buf += (((i, j), iVal * jVal))
}
j += 1
}
}
buf
}.flatten
}
}
}.reduceByKey(_ + _).map { case ((i, j), sim) =>
MatrixEntry(i.toLong, j.toLong, sim)
}
new CoordinateMatrix(sims, numCols(), numCols())
}
private[mllib] override def toBreeze(): BDM[Double] = {
val m = numRows().toInt
val n = numCols().toInt
val mat = BDM.zeros[Double](m, n)
var i = 0
rows.collect().foreach { vector =>
vector.foreachActive { case (j, v) =>
mat(i, j) = v
}
i += 1
}
mat
}
/** 更新或验证行数 */
private def updateNumRows(m: Long) {
if (nRows <= 0) {
nRows = m
} else {
require(nRows == m,
s"行数 $m 与之前指定或计算的不同: ${nRows}.")
}
}
@Since("1.0.0")
object RowMatrix {
/**
* 从上三角部分填充一个完整的方阵
*/
private def triuToFull(n: Int, U: Array[Double]): Matrix = {
val G = new BDM[Double](n, n)
var row = 0
var col = 0
var idx = 0
var value = 0.0
while (col < n) {
row = 0
while (row < col) {
value = U(idx)
G(row, col) = value
G(col, row) = value
idx += 1
row += 1
}
G(col, col) = U(idx)
idx += 1
col +=1
}
Matrices.dense(n, n, G.data)
}
}
IndexedRowMatrix
Spark的IndexedRowMatrix是用于存储和处理带有索引的稀疏矩阵的数据结构。它将稀疏矩阵表示为一组带有索引的行向量，并以分布式方式存储和操作这些行向量。
原理
IndexedRowMatrix将稀疏矩阵分布在不同的分区上，并通过RDD进行处理。每个分区存储一个或多个带有索引的行向量。在分布式环境中，可以并行地处理行向量，从而高效地进行矩阵运算。 IndexedRowMatrix的数据存储方式是一个键值对RDD，其中键是表示行的索引，值是对应的带有索引的行向量。每个行向量都使用SparseVector或DenseVector表示。
用法示例
下面是使用IndexedRowMatrix的简单示例代码：
object IndexedRowMatrixTest {
def main(args: Array[String]): Unit = {
import org.apache.spark.mllib.linalg.distributed.{CoordinateMatrix, MatrixEntry}
val spark = SparkSession.builder()
.appName("ColumnNameExample")
.master("local")
.getOrCreate()
val sc = spark.sparkContext
// 创建一个RDD表示带有索引的行向量
val rows = sc.parallelize(Seq(
IndexedRow(0, Vectors.dense(1.0, 2.0, 3.0)),
IndexedRow(1, Vectors.dense(4.0, 5.0, 6.0)),
IndexedRow(2, Vectors.dense(7.0, 8.0, 9.0))
))
// 将RDD转换为IndexedRowMatrix
val indexedRowMatrix = new IndexedRowMatrix(rows)
// 显示结果矩阵
indexedRowMatrix.rows.foreach(println)
}
}
//IndexedRow(0,[1.0,2.0,3.0])
//IndexedRow(1,[4.0,5.0,6.0])
//IndexedRow(2,[7.0,8.0,9.0])
在示例中，首先创建了一个RDD，表示一个3x3的稀疏矩阵，每个行向量都带有索引。然后，使用IndexedRowMatrix的构造函数将RDD转换为IndexedRowMatrix对象。接下来，使用multiply方法对IndexedRowMatrix进行乘法运算。最后，使用foreach方法遍历并打印结果矩阵的行。 IndexedRowMatrix提供了其他方法，如计算行数和列数、转换为RowMatrix等，可以根据具体需求进行使用。同样，需要注意在处理大规模矩阵时可能会受到内存限制，需要进行合理分区和计算资源管理。
中文源码
/**
* 表示[[org.apache.spark.mllib.linalg.distributed.IndexedRowMatrix]]的一行
*/
@Since("1.0.0")
case class IndexedRow(index: Long, vector: Vector)
/**
* 表示带有索引行的行导向[[org.apache.spark.mllib.linalg.distributed.DistributedMatrix]]。
*
* @param rows 该矩阵的索引行
* @param nRows 行数。非正数表示未知，行数将由最大行索引加一来确定。
* @param nCols 列数。非正数表示未知，列数将由第一行的大小来确定。
*/
@Since("1.0.0")
class IndexedRowMatrix @Since("1.0.0") (
@Since("1.0.0") val rows: RDD[IndexedRow],
private var nRows: Long,
private var nCols: Int) extends DistributedMatrix {
/** 留出矩阵维度自动确定的替代构造函数。*/
@Since("1.0.0")
def this(rows: RDD[IndexedRow]) = this(rows, 0L, 0)
@Since("1.0.0")
override def numCols(): Long = {
if (nCols <= 0) {
// 如果`rows`为空，则调用`first`会抛出异常。
nCols = rows.first().vector.size
}
nCols
}
@Since("1.0.0")
override def numRows(): Long = {
if (nRows <= 0L) {
// 如果`rows`为空，则调用`reduce`会抛出异常。
nRows = rows.map(_.index).reduce(math.max) + 1L
}
nRows
}
/**
* 使用计算归一化点积的蛮力方法，计算该矩阵各列之间的余弦相似度。
*
* @return 一个n x n的稀疏上三角余弦相似度矩阵，表示该矩阵各列之间的相似度。
*/
@Since("1.6.0")
def columnSimilarities(): CoordinateMatrix = {
toRowMatrix().columnSimilarities()
}
/**
* 删除行索引并将此矩阵转换为
* [[org.apache.spark.mllib.linalg.distributed.RowMatrix]]。
*/
@Since("1.0.0")
def toRowMatrix(): RowMatrix = {
new RowMatrix(rows.map(_.vector), 0L, nCols)
}
/**
* 转换为BlockMatrix。使用1024 x 1024的块大小。
*/
@Since("1.3.0")
def toBlockMatrix(): BlockMatrix = {
toBlockMatrix(1024, 1024)
}
/**
* 转换为BlockMatrix。根据行的稀疏程度，块可以是稀疏或稠密的。
* @param rowsPerBlock 每个块的行数。底部边缘的块可能具有更小的值。必须是大于0的整数值。
* @param colsPerBlock 每个块的列数。右侧边缘的块可能具有更小的值。必须是大于0的整数值。
* @return 一个[[BlockMatrix]]
*/
@Since("1.3.0")
def toBlockMatrix(rowsPerBlock: Int, colsPerBlock: Int): BlockMatrix = {
require(rowsPerBlock > 0,
s"rowsPerBlock必须大于0。rowsPerBlock: $rowsPerBlock")
require(colsPerBlock > 0,
s"colsPerBlock必须大于0。colsPerBlock: $colsPerBlock")
val m = numRows()
val n = numCols()
// 由于块矩阵需要一个整数行索引
require(math.ceil(m.toDouble / rowsPerBlock) <= Int.MaxValue,
"行数除以rowsPerBlock不能超过最大整数。")
// 当m % rowsPerBlock != 0或n % colsPerBlock != 0时，余数计算才有意义
val remainderRowBlockIndex = m / rowsPerBlock
val remainderColBlockIndex = n / colsPerBlock
val remainderRowBlockSize = (m % rowsPerBlock).toInt
val remainderColBlockSize = (n % colsPerBlock).toInt
val numRowBlocks = math.ceil(m.toDouble / rowsPerBlock).toInt
val numColBlocks = math.ceil(n.toDouble / colsPerBlock).toInt
val blocks = rows.flatMap { ir: IndexedRow =>
val blockRow = ir.index / rowsPerBlock
val rowInBlock = ir.index % rowsPerBlock
ir.vector match {
case SparseVector(size, indices, values) =>
indices.zip(values).map { case (index, value) =>
val blockColumn = index / colsPerBlock
val columnInBlock = index % colsPerBlock
((blockRow.toInt, blockColumn.toInt), (rowInBlock.toInt, Array((value, columnInBlock))))
}
case DenseVector(values) =>
values.grouped(colsPerBlock)
.zipWithIndex
.map { case (values, blockColumn) =>
((blockRow.toInt, blockColumn), (rowInBlock.toInt, values.zipWithIndex))
}
}
}.groupByKey(GridPartitioner(numRowBlocks, numColBlocks, rows.getNumPartitions)).map {
case ((blockRow, blockColumn), itr) =>
val actualNumRows =
if (blockRow == remainderRowBlockIndex) remainderRowBlockSize else rowsPerBlock
val actualNumColumns =
if (blockColumn == remainderColBlockIndex) remainderColBlockSize else colsPerBlock
val arraySize = actualNumRows * actualNumColumns
val matrixAsArray = new Array[Double](arraySize)
var countForValues = 0
itr.foreach { case (rowWithinBlock, valuesWithColumns) =>
valuesWithColumns.foreach { case (value, columnWithinBlock) =>
matrixAsArray.update(columnWithinBlock * actualNumRows + rowWithinBlock, value)
countForValues += 1
}
}
val denseMatrix = new DenseMatrix(actualNumRows, actualNumColumns, matrixAsArray)
val finalMatrix = if (countForValues / arraySize.toDouble >= 0.1) {
denseMatrix
} else {
denseMatrix.toSparse
}
((blockRow, blockColumn), finalMatrix)
}
new BlockMatrix(blocks, rowsPerBlock, colsPerBlock, m, n)
}
/**
* 将此矩阵转换为
* [[org.apache.spark.mllib.linalg.distributed.CoordinateMatrix]]。
*/
@Since("1.3.0")
def toCoordinateMatrix(): CoordinateMatrix = {
val entries = rows.flatMap { row =>
val rowIndex = row.index
row.vector match {
case SparseVector(size, indices, values) =>
Iterator.tabulate(indices.length)(i => MatrixEntry(rowIndex, indices(i), values(i)))
case DenseVector(values) =>
Iterator.tabulate(values.length)(i => MatrixEntry(rowIndex, i, values(i)))
}
}
new CoordinateMatrix(entries, numRows(), numCols())
}
}
/**
* 计算此IndexedRowMatrix的奇异值分解。
* 将此矩阵记为A（m x n），将计算矩阵U，S，V使得A = U * S * V'。
*
* 此方法的成本和实现与[[org.apache.spark.mllib.linalg.distributed.RowMatrix]]中完全相同
* 在增加指标的情况下。
*
* 最多返回k个最大的非零奇异值和相关向量。
* 如果有k个这样的值，那么返回的维度将是：
*
* U是一个大小为m x k的[[org.apache.spark.mllib.linalg.distributed.IndexedRowMatrix]]，
* 满足 U'U = eye(k)，
* s是一个大小为k的向量，按降序保存奇异值，
* V是一个大小为n x k的本地矩阵，满足 V'V = eye(k)。
*
* @param k 要保留的奇异值的数量。如果有数值上的零奇异值，则可能返回少于k个。请参阅rCond。
* @param computeU 是否计算U
* @param rCond 倒数条件数。小于rCond * sigma(0)的所有奇异值都被视为零，其中sigma(0)是最大的奇异值。
* @return 奇异值分解（U，s，V）
*/
@Since("1.0.0")
def computeSVD(
k: Int,
computeU: Boolean = false,
rCond: Double = 1e-9): SingularValueDecomposition[IndexedRowMatrix, Matrix] = {
val n = numCols().toInt
require(k > 0 && k <= n, s"请求了k个奇异值，但得到的k=$k和numCols=$n。")
val indices = rows.map(_.index)
val svd = toRowMatrix().computeSVD(k, computeU, rCond)
val U = if (computeU) {
val indexedRows = indices.zip(svd.U.rows).map { case (i, v) =>
IndexedRow(i, v)
}
new IndexedRowMatrix(indexedRows, nRows, svd.U.numCols().toInt)
} else {
null
}
SingularValueDecomposition(U, svd.s, svd.V)
}
/**
* 将此矩阵与右侧的本地矩阵相乘。
*
* @param B 行数必须与此矩阵的列数相同的本地矩阵
* @return 表示乘积的IndexedRowMatrix，保留分区
*/
@Since("1.0.0")
def multiply(B: Matrix): IndexedRowMatrix = {
val mat = toRowMatrix().multiply(B)
val indexedRows = rows.map(_.index).zip(mat.rows).map { case (i, v) =>
IndexedRow(i, v)
}
new IndexedRowMatrix(indexedRows, nRows, B.numCols)
}
/**
* 计算Gramian矩阵`A^T A`。
*
* @note 无法在列数超过65535的矩阵上计算。
*/
@Since("1.0.0")
def computeGramianMatrix(): Matrix = {
toRowMatrix().computeGramianMatrix()
}
private[mllib] override def toBreeze(): BDM[Double] = {
val m = numRows().toInt
val n = numCols().toInt
val mat = BDM.zeros[Double](m, n)
rows.collect().foreach { case IndexedRow(rowIndex, vector) =>
val i = rowIndex.toInt
vector.foreachActive { case (j, v) =>
mat(i, j) = v
}
}
mat
}
}
BlockMatrix
Spark的BlockMatrix是用于存储和处理分布式稠密矩阵的数据结构。它将稠密矩阵划分为一组方块，并以分布式方式存储和操作这些方块。
原理
BlockMatrix将整个稠密矩阵分割为多个方块，每个方块存储为一个DenseMatrix对象。在分布式环境中，矩阵方块被分布在不同的分区中，并通过RDD进行处理。这使得可以并行地对矩阵进行操作。 BlockMatrix的数据存储方式是一个键值对RDD，其中键是表示方块位置的元组(blockRowIndex, blockColIndex)，值是方块对应的DenseMatrix对象。
用法示例
下面是使用BlockMatrix的简单示例代码：
object BlockMatrixTest{
def main(args: Array[String]): Unit = {
import org.apache.spark.mllib.linalg.distributed.{CoordinateMatrix, MatrixEntry}
val spark = SparkSession.builder()
.appName("BlockMatrixTest")
.master("local")
.getOrCreate()
val sc = spark.sparkContext
import org.apache.spark.mllib.linalg._
import org.apache.spark.mllib.linalg.distributed.{BlockMatrix, CoordinateMatrix, MatrixEntry}
// 创建一个CoordinateMatrix对象
val entries: RDD[MatrixEntry] = sc.parallelize(Seq(
MatrixEntry(0, 0, 1.0),
MatrixEntry(0, 1, 2.0),
MatrixEntry(1, 0, 3.0),
MatrixEntry(1, 1, 4.0)
))
val coordMatrix: CoordinateMatrix = new CoordinateMatrix(entries)
// 将CoordinateMatrix转换为BlockMatrix
val blockSize = 2
val blockMatrix: BlockMatrix = coordMatrix.toBlockMatrix(rowsPerBlock = blockSize, colsPerBlock = blockSize)
// 进行矩阵乘法
val result: BlockMatrix = blockMatrix.multiply(blockMatrix)
// 显示结果矩阵
result.toLocalMatrix().toArray.foreach(println)
}
}
//7.0
//15.0
//10.0
//22.0
在示例中，首先创建了一个CoordinateMatrix对象，该对象表示一个2x2的稀疏矩阵。然后，使用toBlockMatrix方法将CoordinateMatrix转换为BlockMatrix，并指定每个方块的大小为2。然后，使用multipy方法对BlockMatrix进行乘法运算。最后，使用toLocalMatrix方法将结果BlockMatrix转换为本地矩阵，并以数组的形式进行打印。 注意，BlockMatrix的大小和分区可以根据具体情况进行调整，以便适应不同的任务需求。
中文源码
/**
* 使用规则的网格来对坐标进行划分的网格划分器。
*
* @param rows 行数。
* @param cols 列数。
* @param rowsPerPart 每个分区的行数，下边缘的可能较少。
* @param colsPerPart 每个分区的列数，右边缘的可能较少。
*/
private[mllib] class GridPartitioner(
val rows: Int,
val cols: Int,
val rowsPerPart: Int,
val colsPerPart: Int) extends Partitioner {
require(rows > 0)
require(cols > 0)
require(rowsPerPart > 0)
require(colsPerPart > 0)
private val rowPartitions = math.ceil(rows * 1.0 / rowsPerPart).toInt
private val colPartitions = math.ceil(cols * 1.0 / colsPerPart).toInt
override val numPartitions: Int = rowPartitions * colPartitions
/**
* 返回输入坐标所属分区的索引。
*
* @param key 分区的ID i（通过此方法计算对于坐标（i，j）在`simulateMultiply`中，坐标（i，j）或元组（i，j，k），其中k是在乘法中使用的内部索引。k在计算分区时被忽略。
* @return 坐标所属分区的索引。
*/
override def getPartition(key: Any): Int = {
key match {
case i: Int => i
case (i: Int, j: Int) =>
getPartitionId(i, j)
case (i: Int, j: Int, _: Int) =>
getPartitionId(i, j)
case _ =>
throw new IllegalArgumentException(s"无法识别的键: $key.")
}
}
/** 将子矩阵分区为具有相邻子矩阵的块。 */
private def getPartitionId(i: Int, j: Int): Int = {
require(0 <= i && i < rows, s"行索引 $i 超出范围 [0, $rows).")
require(0 <= j && j < cols, s"列索引 $j 超出范围 [0, $cols).")
i / rowsPerPart + j / colsPerPart * rowPartitions
}
override def equals(obj: Any): Boolean = {
obj match {
case r: GridPartitioner =>
(this.rows == r.rows) && (this.cols == r.cols) &&
(this.rowsPerPart == r.rowsPerPart) && (this.colsPerPart == r.colsPerPart)
case _ =>
false
}
}
override def hashCode: Int = {
com.google.common.base.Objects.hashCode(
rows: java.lang.Integer,
cols: java.lang.Integer,
rowsPerPart: java.lang.Integer,
colsPerPart: java.lang.Integer)
}
}
private[mllib] object GridPartitioner {
/** 创建一个新的[[GridPartitioner]]实例。 */
def apply(rows: Int, cols: Int, rowsPerPart: Int, colsPerPart: Int): GridPartitioner = {
new GridPartitioner(rows, cols, rowsPerPart, colsPerPart)
}
/** 使用输入的建议分区数创建一个新的[[GridPartitioner]]实例。 */
def apply(rows: Int, cols: Int, suggestedNumPartitions: Int): GridPartitioner = {
require(suggestedNumPartitions > 0)
val scale = 1.0 / math.sqrt(suggestedNumPartitions)
val rowsPerPart = math.round(math.max(scale * rows, 1.0)).toInt
val colsPerPart = math.round(math.max(scale * cols, 1.0)).toInt
new GridPartitioner(rows, cols, rowsPerPart, colsPerPart)
}
}
/**
* 以本地矩阵块的形式表示分布式矩阵。
*
* @param blocks 形成此分布式矩阵的子矩阵块的RDD(((blockRowIndex, blockColIndex), sub-matrix))。如果存在多个具有相同索引的块，则像add和multiply等操作的结果将是不可预测的。
* @param rowsPerBlock 构成每个块的行数。最终行组成的块不需要具有给定的行数
* @param colsPerBlock 构成每个块的列数。最终列组成的块不需要具有给定的列数
* @param nRows 此矩阵的行数。如果提供的值小于或等于零，则在调用`numRows`时将计算行数。
* @param nCols 此矩阵的列数。如果提供的值小于或等于零，则在调用`numCols`时将计算列数。
*/
@Since("1.3.0")
class BlockMatrix @Since("1.3.0") (
@Since("1.3.0") val blocks: RDD[((Int, Int), Matrix)],
@Since("1.3.0") val rowsPerBlock: Int,
@Since("1.3.0") val colsPerBlock: Int,
private var nRows: Long,
private var nCols: Long) extends DistributedMatrix with Logging {
private type MatrixBlock = ((Int, Int), Matrix) // ((blockRowIndex, blockColIndex), sub-matrix)
/**
* BlockMatrix的另一种构造函数，不需要输入行数和列数。
*
* @param blocks 形成此分布式矩阵的子矩阵块的RDD(((blockRowIndex, blockColIndex), sub-matrix))。如果存在多个具有相同索引的块，则像add和multiply等操作的结果将是不可预测的。
* @param rowsPerBlock 构成每个块的行数。最终行组成的块不需要具有给定的行数
* @param colsPerBlock 构成每个块的列数。最终列组成的块不需要具有给定的列数
*/
@Since("1.3.0")
def this(
blocks: RDD[((Int, Int), Matrix)],
rowsPerBlock: Int,
colsPerBlock: Int) = {
this(blocks, rowsPerBlock, colsPerBlock, 0L, 0L)
}
@Since("1.3.0")
override def numRows(): Long = {
if (nRows <= 0L) estimateDim()
nRows
}
@Since("1.3.0")
override def numCols(): Long = {
if (nCols <= 0L) estimateDim()
nCols
}
@Since("1.3.0")
val numRowBlocks = math.ceil(numRows() * 1.0 / rowsPerBlock).toInt
@Since("1.3.0")
val numColBlocks = math.ceil(numCols() * 1.0 / colsPerBlock).toInt
private[mllib] def createPartitioner(): GridPartitioner =
GridPartitioner(numRowBlocks, numColBlocks, suggestedNumPartitions = blocks.partitions.length)
private lazy val blockInfo = blocks.mapValues(block => (block.numRows, block.numCols)).cache()
/** 估计矩阵的维度。 */
private def estimateDim(): Unit = {
val (rows, cols) = blockInfo.map { case ((blockRowIndex, blockColIndex), (m, n)) =>
(blockRowIndex.toLong * rowsPerBlock + m,
blockColIndex.toLong * colsPerBlock + n)
}.reduce { (x0, x1) =>
(math.max(x0._1, x1._1), math.max(x0._2, x1._2))
}
if (nRows <= 0L) nRows = rows
assert(rows <= nRows, s"行数 $rows 大于声明的 $nRows。")
if (nCols <= 0L) nCols = cols
assert(cols <= nCols, s"列数 $cols 大于声明的 $nCols。")
}
}
/**
* 验证块矩阵信息与矩阵数据（`blocks`）是否一致，如果发现任何错误则抛出异常。
*/
@Since("1.3.0")
def validate(): Unit = {
logDebug("验证块矩阵...")
// 检查矩阵是否大于声明的维度
estimateDim()
logDebug("块矩阵维度正确...")
// 检查是否有相同索引的多个MatrixBlocks。
blockInfo.countByKey().foreach { case (key, cnt) =>
if (cnt > 1) {
throw new SparkException(s"找到多个具有索引$key的MatrixBlocks，请删除具有重复索引的块。")
}
}
logDebug("MatrixBlock索引正确...")
// 检查每个MatrixBlock（除边缘）是否具有rowsPerBlock x colsPerBlock的维度
// 第一个元组是索引，第二个元组是MatrixBlock的维度
val dimensionMsg = s"维度与rowsPerBlock: $rowsPerBlock和colsPerBlock: $colsPerBlock不同。" +
s"右边和底边的块可以具有较小的维度。您可以使用repartition方法来解决此问题。"
blockInfo.foreach { case ((blockRowIndex, blockColIndex), (m, n)) =>
if ((blockRowIndex < numRowBlocks - 1 && m != rowsPerBlock) ||
(blockRowIndex == numRowBlocks - 1 && (m <= 0 || m > rowsPerBlock))) {
throw new SparkException(s"($blockRowIndex, $blockColIndex)处的MatrixBlock具有" + dimensionMsg)
}
if ((blockColIndex < numColBlocks - 1 && n != colsPerBlock) ||
(blockColIndex == numColBlocks - 1 && (n <= 0 || n > colsPerBlock))) {
throw new SparkException(s"($blockRowIndex, $blockColIndex)处的MatrixBlock具有" + dimensionMsg)
}
}
logDebug("MatrixBlock维度正确...")
logDebug("块矩阵有效！")
}
/** 缓存底层RDD。 */
@Since("1.3.0")
def cache(): this.type = {
blocks.cache()
this
}
/** 使用指定的存储级别持久化底层RDD。 */
@Since("1.3.0")
def persist(storageLevel: StorageLevel): this.type = {
blocks.persist(storageLevel)
this
}
/** 转换为CoordinateMatrix。 */
@Since("1.3.0")
def toCoordinateMatrix(): CoordinateMatrix = {
val entryRDD = blocks.flatMap { case ((blockRowIndex, blockColIndex), mat) =>
val rowStart = blockRowIndex.toLong * rowsPerBlock
val colStart = blockColIndex.toLong * colsPerBlock
val entryValues = new ArrayBuffer[MatrixEntry]()
mat.foreachActive { (i, j, v) =>
if (v != 0.0) entryValues += new MatrixEntry(rowStart + i, colStart + j, v)
}
entryValues
}
new CoordinateMatrix(entryRDD, numRows(), numCols())
}
/** 转换为IndexedRowMatrix。列数必须在整数范围内。 */
@Since("1.3.0")
def toIndexedRowMatrix(): IndexedRowMatrix = {
val cols = numCols().toInt
require(cols < Int.MaxValue, s"列数应小于Int.MaxValue($cols)。")
val rows = blocks.flatMap { case ((blockRowIdx, blockColIdx), mat) =>
mat.rowIter.zipWithIndex.map {
case (vector, rowIdx) =>
blockRowIdx * rowsPerBlock + rowIdx -> ((blockColIdx, vector.asBreeze))
}
}.groupByKey().map { case (rowIdx, vectors) =>
val numberNonZeroPerRow = vectors.map(_._2.activeSize).sum.toDouble / cols.toDouble
val wholeVector = if (numberNonZeroPerRow <= 0.1) { // 1/10th nnz的稀疏向量
BSV.zeros[Double](cols)
} else {
BDV.zeros[Double](cols)
}
vectors.foreach { case (blockColIdx: Int, vec: BV[_]) =>
val offset = colsPerBlock * blockColIdx
wholeVector(offset until Math.min(cols, offset + colsPerBlock)) := vec
}
new IndexedRow(rowIdx, Vectors.fromBreeze(wholeVector))
}
new IndexedRowMatrix(rows)
}
/**
* 在驱动程序中将分布式矩阵收集为DenseMatrix。
*/
@Since("1.3.0")
def toLocalMatrix(): Matrix = {
require(numRows() < Int.MaxValue, "这个矩阵的行数应小于Int.MaxValue。" +
s"当前行数：${numRows()}")
require(numCols() < Int.MaxValue, "这个矩阵的列数应小于Int.MaxValue。" +
s"当前列数：${numCols()}")
require(numRows() * numCols() < Int.MaxValue, "值数组的长度必须小于Int.MaxValue。" +
s"当前行数 * 列数：${numRows() * numCols()}")
val m = numRows().toInt
val n = numCols().toInt
val mem = m * n / 125000
if (mem > 500) logWarning(s"存储这个矩阵将需要 $mem MB 的内存!")
val localBlocks = blocks.collect()
val values = new Array[Double](m * n)
localBlocks.foreach { case ((blockRowIndex, blockColIndex), submat) =>
val rowOffset = blockRowIndex * rowsPerBlock
val colOffset = blockColIndex * colsPerBlock
submat.foreachActive { (i, j, v) =>
val indexOffset = (j + colOffset) * m + rowOffset + i
values(indexOffset) = v
}
}
new DenseMatrix(m, n, values)
}
/**
* 转置这个`BlockMatrix`。返回一个共享相同底层数据的新的`BlockMatrix`实例。这是一个惰性操作。
*/
@Since("1.3.0")
def transpose: BlockMatrix = {
val transposedBlocks = blocks.map { case ((blockRowIndex, blockColIndex), mat) =>
((blockColIndex, blockRowIndex), mat.transpose)
}
new BlockMatrix(transposedBlocks, colsPerBlock, rowsPerBlock, nCols, nRows)
}
/** 收集数据并组装为本地稠密breeze矩阵（仅用于测试）。 */
private[mllib] def toBreeze(): BDM[Double] = {
val localMat = toLocalMatrix()
new BDM[Double](localMat.numRows, localMat.numCols, localMat.toArray)
}
/**
* 给定与`this`和`other`具有兼容维度和兼容块维度的矩阵，对它们的相应块应用二元函数。
*
* @param other `binMap`指定的运算符的第二个BlockMatrix参数
* @param binMap 一个接受两个breeze矩阵并返回一个breeze矩阵的函数
* @return 一个[[BlockMatrix]]，其块是`this`和`other`的指定二元映射在块上的结果
* 注意：`blockMap`仅适用于`add`和`subtract`方法，并不支持诸如(a, b) => -a + b的运算符
* TODO: 优化使用零矩阵的存储效率。
*/
private[mllib] def blockMap(
other: BlockMatrix,
binMap: (BM[Double], BM[Double]) => BM[Double]): BlockMatrix = {
require(numRows() == other.numRows(), "两个矩阵必须具有相同的行数. " +
s"A.numRows: ${numRows()}, B.numRows: ${other.numRows()}")
require(numCols() == other.numCols(), "两个矩阵必须具有相同的列数. " +
s"A.numCols: ${numCols()}, B.numCols: ${other.numCols()}")
if (rowsPerBlock == other.rowsPerBlock && colsPerBlock == other.colsPerBlock) {
val newBlocks = blocks.cogroup(other.blocks, createPartitioner())
.map { case ((blockRowIndex, blockColIndex), (a, b)) =>
if (a.size > 1 || b.size > 1) {
throw new SparkException("具有以下索引的MatrixBlocks有多个: " +
s"($blockRowIndex, $blockColIndex)。请删除它们。")
}
if (a.isEmpty) {
val zeroBlock = BM.zeros[Double](b.head.numRows, b.head.numCols)
val result = binMap(zeroBlock, b.head.asBreeze)
new MatrixBlock((blockRowIndex, blockColIndex), Matrices.fromBreeze(result))
} else if (b.isEmpty) {
new MatrixBlock((blockRowIndex, blockColIndex), a.head)
} else {
val result = binMap(a.head.asBreeze, b.head.asBreeze)
new MatrixBlock((blockRowIndex, blockColIndex), Matrices.fromBreeze(result))
}
}
new BlockMatrix(newBlocks, rowsPerBlock, colsPerBlock, numRows(), numCols())
} else {
throw new SparkException("无法处理具有不同块维度的矩阵")
}
}
/**
* 将给定的块矩阵`other`添加到`this`块矩阵：`this + other`。
* 矩阵必须具有相同的大小和匹配的`rowsPerBlock`和`colsPerBlock`值。
* 如果要添加的块之一是`SparseMatrix`的实例，即使它被添加到`DenseMatrix`中，结果子矩阵也将是`SparseMatrix`。
* 如果两个稠密矩阵相加，输出将仍为`DenseMatrix`。
*/
@Since("1.3.0")
def add(other: BlockMatrix): BlockMatrix =
blockMap(other, (x: BM[Double], y: BM[Double]) => x + y)
/**
* 从`this`块矩阵中减去给定的块矩阵`other`：`this - other`。
* 矩阵必须具有相同的大小和匹配的`rowsPerBlock`和`colsPerBlock`值。
* 如果要减去的块之一是`SparseMatrix`的实例，即使它从`DenseMatrix`中减去，结果子矩阵也将是`SparseMatrix`。
* 如果两个稠密矩阵相减，输出将仍为`DenseMatrix`。
*/
@Since("2.0.0")
def subtract(other: BlockMatrix): BlockMatrix =
blockMap(other, (x: BM[Double], y: BM[Double]) => x - y)
/** 块(i,j) --> 目标分区的集合 */
private type BlockDestinations = Map[(Int, Int), Set[Int]]
/**
* 模拟只使用块索引进行乘法，以便在实际洗牌矩阵时降低通信成本。
* 此矩阵的`colsPerBlock`必须等于`other`的`rowsPerBlock`。
* 暴露给测试。
*
* @param other 要乘的BlockMatrix
* @param partitioner 将用于结果矩阵`C = A * B`的分区器
* @return 一个[[BlockDestinations]]的元组。第一个元素是我们需要洗牌`this`的每个块的分区集合的Map，
* 第二个元素是`other`的Map。
*/
private[distributed] def simulateMultiply(
other: BlockMatrix,
partitioner: GridPartitioner,
midDimSplitNum: Int): (BlockDestinations, BlockDestinations) = {
val leftMatrix = blockInfo.keys.collect()
val rightMatrix = other.blockInfo.keys.collect()
val rightCounterpartsHelper = rightMatrix.groupBy(_._1).mapValues(_.map(_._2))
val leftDestinations = leftMatrix.map { case (rowIndex, colIndex) =>
val rightCounterparts = rightCounterpartsHelper.getOrElse(colIndex, Array.empty[Int])
val partitions = rightCounterparts.map(b => partitioner.getPartition((rowIndex, b)))
val midDimSplitIndex = colIndex % midDimSplitNum
((rowIndex, colIndex),
partitions.toSet.map((pid: Int) => pid * midDimSplitNum + midDimSplitIndex))
}.toMap
val leftCounterpartsHelper = leftMatrix.groupBy(_._2).mapValues(_.map(_._1))
val rightDestinations = rightMatrix.map { case (rowIndex, colIndex) =>
val leftCounterparts = leftCounterpartsHelper.getOrElse(rowIndex, Array.empty[Int])
val partitions = leftCounterparts.map(b => partitioner.getPartition((b, colIndex)))
val midDimSplitIndex = rowIndex % midDimSplitNum
((rowIndex, colIndex),
partitions.toSet.map((pid: Int) => pid * midDimSplitNum + midDimSplitIndex))
}.toMap
(leftDestinations, rightDestinations)
}
/**
* 将此[[BlockMatrix]]与`other`，另一个[[BlockMatrix]]进行左乘。此矩阵的`colsPerBlock`必须等于`other`的`rowsPerBlock`。
* 如果`other`包含`SparseMatrix`，则必须将其转换为`DenseMatrix`。输出的[[BlockMatrix]]将仅包含`DenseMatrix`的块。
* 这可能会导致一些性能问题，直到添加了两个稀疏矩阵相乘的支持。
*
* @note 乘法的行为在1.6.0中发生了变化。`multiply`在存在重复索引块时曾经抛出错误。现在，具有重复索引的块将会相加。
*/
@Since("1.3.0")
def multiply(other: BlockMatrix): BlockMatrix = {
multiply(other, 1)
}
/**
* 将此[[BlockMatrix]]与`other`，另一个[[BlockMatrix]]进行左乘。此矩阵的`colsPerBlock`必须等于`other`的`rowsPerBlock`。
* 如果`other`包含`SparseMatrix`，则必须将其转换为`DenseMatrix`。输出的[[BlockMatrix]]将仅包含`DenseMatrix`的块。
* 这可能会导致一些性能问题，直到添加了两个稀疏矩阵相乘的支持。
* 具有重复索引的块将会相加。
*
* @param other 矩阵 `B` 在 `A * B = C` 中
* @param numMidDimSplits 在进行乘法时，沿中间维度切分的切分数。
* 例如，当将大小为 `m x n` 的矩阵 `A` 与大小为 `n x k` 的矩阵 `B` 相乘时，
* 此参数配置了在分组矩阵时使用的并行性。并行性将从 `m x k` 增加到 `m x k x numMidDimSplits`，
* 在某些情况下还可以减少总的洗牌数据。
*/
@Since("2.2.0")
def multiply(other: BlockMatrix, numMidDimSplits: Int): BlockMatrix = {
require(numCols() == other.numRows(), "A的列数和B的行数必须相等。A.numCols: ${numCols()}, B.numRows: ${other.numRows()}。" +
"如果您认为它们应该相等，请在初始化时明确设置A和B的维度。")
require(numMidDimSplits > 0, "numMidDimSplits应该是一个正整数。")
if (colsPerBlock == other.rowsPerBlock) {
val resultPartitioner = GridPartitioner(numRowBlocks, other.numColBlocks,
math.max(blocks.partitions.length, other.blocks.partitions.length))
val (leftDestinations, rightDestinations) = simulateMultiply(other, resultPartitioner, numMidDimSplits)
// A的每个块必须与B的列中对应的块相乘
val flatA = blocks.flatMap { case ((blockRowIndex, blockColIndex), block) =>
val destinations = leftDestinations.getOrElse((blockRowIndex, blockColIndex), Set.empty)
destinations.map(j => (j, (blockRowIndex, blockColIndex, block)))
}
// B的每个块必须与A的每一行中对应的块相乘
val flatB = other.blocks.flatMap { case ((blockRowIndex, blockColIndex), block) =>
val destinations = rightDestinations.getOrElse((blockRowIndex, blockColIndex), Set.empty)
destinations.map(j => (j, (blockRowIndex, blockColIndex, block)))
}
val intermediatePartitioner = new Partitioner {
override def numPartitions: Int = resultPartitioner.numPartitions * numMidDimSplits
override def getPartition(key: Any): Int = key.asInstanceOf[Int]
}
val newBlocks = flatA.cogroup(flatB, intermediatePartitioner).flatMap { case (pId, (a, b)) =>
a.flatMap { case (leftRowIndex, leftColIndex, leftBlock) =>
b.filter(_._1 == leftColIndex).map { case (rightRowIndex, rightColIndex, rightBlock) =>
val C = rightBlock match {
case dense: DenseMatrix => leftBlock.multiply(dense)
case sparse: SparseMatrix => leftBlock.multiply(sparse.toDense)
case _ => throw new SparkException(s"未识别的矩阵类型 ${rightBlock.getClass}。")
}
((leftRowIndex, rightColIndex), C.asBreeze)
}
}
}.reduceByKey(resultPartitioner, (a, b) => a + b).mapValues(Matrices.fromBreeze)
// TODO: 尝试使用aggregateByKey而不是reduceByKey以摆脱中间矩阵
new BlockMatrix(newBlocks, rowsPerBlock, other.colsPerBlock, numRows(), other.numCols())
} else {
throw new SparkException("A的colsPerBlock与B的rowsPerBlock不匹配。" +
s"A.colsPerBlock: $colsPerBlock, B.rowsPerBlock: ${other.rowsPerBlock}")
}
}
2308
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
