Distributed Stochastic Singular Value Decomposition

Intro

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.

Modified SSVD Algorithm

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}\):

1. 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].

2. \(\mathbf{Y=A\boldsymbol{\Omega}},\,\mathbf{Y}\in\mathbb{R}^{m\times\left(k+p\right)}\)

3. 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}\)

4. \(\mathbf{B}_{0}=\mathbf{Q}^{\top}\mathbf{A}:\,\,\mathbf{B}\in\mathbb{R}^{\left(k+p\right)\times n}\).

5. 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).

6. 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)}\).

7. Singular values \(\mathbf{\boldsymbol{\Sigma}}=\boldsymbol{\Lambda}^{0.5}\), or, in other words, \(s_{i}=\sqrt{\sigma_{i}}\).

8. If needed, compute \(\mathbf{U}=\mathbf{Q}\hat{\mathbf{U}}\).

9. 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}\).

Implementation

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.

Usage

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)

References

approximations of matrices](http://amath.colorado.edu/faculty/martinss/Pubs/2012_halko_dissertation.pdf)