layout: doc-page title: Distributed Stochastic Singular Value Decomposition
Mahout has a distributed implementation of Stochastic Singular Value Decomposition [1] using the parallelization strategy comprehensively defined in Nathan Halko's dissertation “Randomized methods for computing low-rank approximations of matrices” [2].
Given an (m\times n) matrix (\mathbf{A}), a target rank (k\in\mathbb{N}{1}) , an oversampling parameter (p\in\mathbb{N}{1}), and the number of additional power iterations (q\in\mathbb{N}_{0}), this procedure computes an (m\times\left(k+p\right)) SVD (\mathbf{A\approx U}\boldsymbol{\Sigma}\mathbf{V}^{\top}):
Create seed for random (n\times\left(k+p\right)) matrix (\boldsymbol{\Omega}). The seed defines matrix (\mathbf{\Omega}) using Gaussian unit vectors per one of suggestions in [Halko, Martinsson, Tropp].
(\mathbf{Y=A\boldsymbol{\Omega}},,\mathbf{Y}\in\mathbb{R}^{m\times\left(k+p\right)})
Column-orthonormalize (\mathbf{Y}\rightarrow\mathbf{Q}) by computing thin decomposition (\mathbf{Y}=\mathbf{Q}\mathbf{R}). Also, (\mathbf{Q}\in\mathbb{R}^{m\times\left(k+p\right)},,\mathbf{R}\in\mathbb{R}^{\left(k+p\right)\times\left(k+p\right)}); denoted as (\mathbf{Q}=\mbox{qr}\left(\mathbf{Y}\right).\mathbf{Q})
(\mathbf{B}_{0}=\mathbf{Q}^{\top}\mathbf{A}:,,\mathbf{B}\in\mathbb{R}^{\left(k+p\right)\times n}).
If (q>0) repeat: for (i=1..q): (\mathbf{B}{i}^{\top}=\mathbf{A}^{\top}\mbox{qr}\left(\mathbf{A}\mathbf{B}{i-1}^{\top}\right).\mathbf{Q}) (power iterations step).
Compute Eigensolution of a small Hermitian (\mathbf{B}{q}\mathbf{B}{q}^{\top}=\mathbf{\hat{U}}\boldsymbol{\Lambda}\mathbf{\hat{U}}^{\top}), (\mathbf{B}{q}\mathbf{B}{q}^{\top}\in\mathbb{R}^{\left(k+p\right)\times\left(k+p\right)}).
Singular values (\mathbf{\boldsymbol{\Sigma}}=\boldsymbol{\Lambda}^{0.5}), or, in other words, (s_{i}=\sqrt{\sigma_{i}}).
If needed, compute (\mathbf{U}=\mathbf{Q}\hat{\mathbf{U}}).
If needed, compute (\mathbf{V}=\mathbf{B}_{q}^{\top}\hat{\mathbf{U}}\boldsymbol{\Sigma}^{-1}). Another way is (\mathbf{V}=\mathbf{A}^{\top}\mathbf{U}\boldsymbol{\Sigma}^{-1}).
Mahout dssvd(...)
is implemented in the mahout math-scala
algebraic optimizer which translates Mahout's R-like linear algebra operators into a physical plan for both Spark and H2O distributed engines.
def dssvd[K: ClassTag](drmA: DrmLike[K], k: Int, p: Int = 15, q: Int = 0): (DrmLike[K], DrmLike[Int], Vector) = { val drmAcp = drmA.checkpoint() val m = drmAcp.nrow val n = drmAcp.ncol assert(k <= (m min n), "k cannot be greater than smaller of m, n.") val pfxed = safeToNonNegInt((m min n) - k min p) // Actual decomposition rank val r = k + pfxed // We represent Omega by its seed. val omegaSeed = RandomUtils.getRandom().nextInt() // Compute Y = A*Omega. var drmY = drmAcp.mapBlock(ncol = r) { case (keys, blockA) => val blockY = blockA %*% Matrices.symmetricUniformView(n, r, omegaSeed) keys -> blockY } var drmQ = dqrThin(drmY.checkpoint())._1 // Checkpoint Q if last iteration if (q == 0) drmQ = drmQ.checkpoint() var drmBt = drmAcp.t %*% drmQ // Checkpoint B' if last iteration if (q == 0) drmBt = drmBt.checkpoint() for (i <- 0 until q) { drmY = drmAcp %*% drmBt drmQ = dqrThin(drmY.checkpoint())._1 // Checkpoint Q if last iteration if (i == q - 1) drmQ = drmQ.checkpoint() drmBt = drmAcp.t %*% drmQ // Checkpoint B' if last iteration if (i == q - 1) drmBt = drmBt.checkpoint() } val (inCoreUHat, d) = eigen(drmBt.t %*% drmBt) val s = d.sqrt // Since neither drmU nor drmV are actually computed until actually used // we don't need the flags instructing compute (or not compute) either of the U,V outputs val drmU = drmQ %*% inCoreUHat val drmV = drmBt %*% (inCoreUHat %*%: diagv(1 /: s)) (drmU(::, 0 until k), drmV(::, 0 until k), s(0 until k)) }
Note: As a side effect of checkpointing, U and V values are returned as logical operators (i.e. they are neither checkpointed nor computed). Therefore there is no physical work actually done to compute (\mathbf{U}) or (\mathbf{V}) until they are used in a subsequent expression.
The scala dssvd(...)
method can easily be called in any Spark or H2O application built with the math-scala
library and the corresponding Spark
or H2O
engine module as follows:
import org.apache.mahout.math._ import decompositions._ import drm._ val(drmU, drmV, s) = dssvd(drma, k = 40, q = 1)
[1]: Mahout Scala and Mahout Spark Bindings for Linear Algebra Subroutines
[2]: Randomized methods for computing low-rank approximations of matrices