splearn/decomposition/truncated_svd.py
def svd_em(blocked_rdd, k, maxiter=20, tol=1e-6, compute_u=True, seed=None):
    """
    Calculate the SVD of a blocked RDD using an expectation maximization
    algorithm (from Roweis, NIPS, 1997) that avoids explicitly
    computing the covariance matrix, returning only the leading k
    singular vectors. Assumes n rows and d columns, does not require
    d^2 to fit into memory on a single machine.
    Parameters
    ----------
    blocked_rdd : ArrayRDD
        ArrayRDD with data points in numpy array blocks
    k : Int
        Number of singular vectors to return
    maxiter : Int, optional, default = 20
        Number of iterations to perform
    tol : Double, optional, default = 1e-5
        Tolerance for stopping iterative updates
    seed : Int, optional, default = None
        Seed for random number generator for initializing subspace
    Returns
    ----------
    u : RDD of blocks
        Left eigenvectors
    s : numpy array
        Singular values
    v : numpy array
        Right eigenvectors
    """

    n, m = blocked_rdd.shape[:2]
    sc = blocked_rdd._rdd.context

    def outerprod(x):
        return x.T.dot(x)

    # global run_sum

    # def accumsum(x):
    #     global run_sum
    #     run_sum += x

    # class MatrixAccum(AccumulatorParam):

    #     def zero(self, value):
    #         return np.zeros(np.shape(value))

    #     def addInPlace(self, val1, val2):
    #         val1 += val2
    #         return val1

    if seed is not None:
        rng = np.random.RandomState(seed)
        c = rng.randn(k, m)
    else:
        c = np.random.randn(k, m)

    iter = 0
    error = 100

    # iteratively update subspace using expectation maximization
    # e-step: x = (cc')^-1 c y
    # m-step: c = y x' (xx')^-1
    while (iter < maxiter) & (error > tol):
        c_old = c

        # pre compute (cc')^-1 c
        c_inv = np.dot(c.T, ln.inv(np.dot(c, c.T)))
        premult1 = sc.broadcast(c_inv)

        # compute (xx')^-1 through a map reduce
        xx = blocked_rdd.map(lambda x: outerprod(safe_sparse_dot(x, premult1.value))) \
                        .treeReduce(add)

        # compute (xx')^-1 using an accumulator
        # run_sum = sc.accumulator(np.zeros((k, k)), MatrixAccum())
        # blocked_rdd.map(lambda x: outerprod(safe_sparse_dot(x, premult1.value))) \
        #            .foreachPartition(lambda l: accumsum(sum(l)))
        # xx = run_sum.value
        xx_inv = ln.inv(xx)

        # pre compute (cc')^-1 c (xx')^-1
        premult2 = blocked_rdd.context.broadcast(np.dot(c_inv, xx_inv))

        # compute the new c through a map reduce
        c = blocked_rdd.map(lambda x: safe_sparse_dot(x.T, safe_sparse_dot(x, premult2.value))) \
                       .treeReduce(add)

        # compute the new c using an accumulator
        # run_sum = sc.accumulator(np.zeros((m, k)), MatrixAccum())
        # blocked_rdd.map(lambda x: safe_sparse_dot(x.T, safe_sparse_dot(x, premult2.value))) \
        #            .foreachPartition(lambda l: accumsum(sum(l)))
        # c = run_sum.value
        c = c.T

        error = np.sum((c - c_old) ** 2)
        iter += 1

    # project data into subspace spanned by columns of c
    # use standard eigendecomposition to recover an orthonormal basis
    c = ln.orth(c.T).T
    cov = blocked_rdd.map(lambda x: safe_sparse_dot(x, c.T)) \
                     .map(lambda x: outerprod(x)) \
                     .treeReduce(add)
    w, v = ln.eig(cov / n)
    w = np.real(w)
    v = np.real(v)
    inds = np.argsort(w)[::-1]
    s = np.sqrt(w[inds[0:k]]) * np.sqrt(n)
    v = np.dot(v[:, inds[0:k]].T, c)
    if compute_u:
        v_broadcasted = blocked_rdd.context.broadcast(v)
        u = blocked_rdd.map(
            lambda x: safe_sparse_dot(x, v_broadcasted.value.T) / s)
        return u, s, v
    else:
        return s, v
def svd_em(blocked_rdd, k, maxiter=20, tol=1e-6, compute_u=True, seed=None):
    """
    Calculate the SVD of a blocked RDD using an expectation maximization
    algorithm (from Roweis, NIPS, 1997) that avoids explicitly
    computing the covariance matrix, returning only the leading k
    singular vectors. Assumes n rows and d columns, does not require
    d^2 to fit into memory on a single machine.
    Parameters
    ----------
    blocked_rdd : ArrayRDD
        ArrayRDD with data points in numpy array blocks
    k : Int
        Number of singular vectors to return
    maxiter : Int, optional, default = 20
        Number of iterations to perform
    tol : Double, optional, default = 1e-5
        Tolerance for stopping iterative updates
    seed : Int, optional, default = None
        Seed for random number generator for initializing subspace
    Returns
    ----------
    u : RDD of blocks
        Left eigenvectors
    s : numpy array
        Singular values
    v : numpy array
        Right eigenvectors
    """

    n, m = blocked_rdd.shape[:2]
    sc = blocked_rdd._rdd.context

    def outerprod(x):
        return x.T.dot(x)

    # global run_sum

    # def accumsum(x):
    #     global run_sum
    #     run_sum += x

    # class MatrixAccum(AccumulatorParam):

    #     def zero(self, value):
    #         return np.zeros(np.shape(value))

    #     def addInPlace(self, val1, val2):
    #         val1 += val2
    #         return val1

    if seed is not None:
        rng = np.random.RandomState(seed)
        c = rng.randn(k, m)
    else:
        c = np.random.randn(k, m)

    iter = 0
    error = 100

    # iteratively update subspace using expectation maximization
    # e-step: x = (cc')^-1 c y
    # m-step: c = y x' (xx')^-1
    while (iter < maxiter) & (error > tol):
        c_old = c

        # pre compute (cc')^-1 c
        c_inv = np.dot(c.T, ln.inv(np.dot(c, c.T)))
        premult1 = sc.broadcast(c_inv)

        # compute (xx')^-1 through a map reduce
        xx = blocked_rdd.map(lambda x: outerprod(safe_sparse_dot(x, premult1.value))) \
                        .treeReduce(add)

        # compute (xx')^-1 using an accumulator
        # run_sum = sc.accumulator(np.zeros((k, k)), MatrixAccum())
        # blocked_rdd.map(lambda x: outerprod(safe_sparse_dot(x, premult1.value))) \
        #            .foreachPartition(lambda l: accumsum(sum(l)))
        # xx = run_sum.value
        xx_inv = ln.inv(xx)

        # pre compute (cc')^-1 c (xx')^-1
        premult2 = blocked_rdd.context.broadcast(np.dot(c_inv, xx_inv))

        # compute the new c through a map reduce
        c = blocked_rdd.map(lambda x: safe_sparse_dot(x.T, safe_sparse_dot(x, premult2.value))) \
                       .treeReduce(add)

        # compute the new c using an accumulator
        # run_sum = sc.accumulator(np.zeros((m, k)), MatrixAccum())
        # blocked_rdd.map(lambda x: safe_sparse_dot(x.T, safe_sparse_dot(x, premult2.value))) \
        #            .foreachPartition(lambda l: accumsum(sum(l)))
        # c = run_sum.value
        c = c.T

        error = np.sum((c - c_old) ** 2)
        iter += 1

    # project data into subspace spanned by columns of c
    # use standard eigendecomposition to recover an orthonormal basis
    c = ln.orth(c.T).T
    cov = blocked_rdd.map(lambda x: safe_sparse_dot(x, c.T)) \
                     .map(lambda x: outerprod(x)) \
                     .treeReduce(add)
    w, v = ln.eig(cov / n)
    w = np.real(w)
    v = np.real(v)
    inds = np.argsort(w)[::-1]
    s = np.sqrt(w[inds[0:k]]) * np.sqrt(n)
    v = np.dot(v[:, inds[0:k]].T, c)
    if compute_u:
        v_broadcasted = blocked_rdd.context.broadcast(v)
        u = blocked_rdd.map(
            lambda x: safe_sparse_dot(x, v_broadcasted.value.T) / s)
        return u, s, v
    else:
        return s, v
Calculate the SVD of a blocked RDD using an expectation maximization
    algorithm (from Roweis, NIPS, 1997) that avoids explicitly
    computing the covariance matrix, returning only the leading k
    singular vectors. Assumes n rows and d columns, does not require
    d^2 to fit into memory on a single machine.
    Parameters
    ----------
    blocked_rdd : ArrayRDD
        ArrayRDD with data points in numpy array blocks
    k : Int
        Number of singular vectors to return
    maxiter : Int, optional, default = 20
        Number of iterations to perform
    tol : Double, optional, default = 1e-5
        Tolerance for stopping iterative updates
    seed : Int, optional, default = None
        Seed for random number generator for initializing subspace
    Returns
    ----------
    u : RDD of blocks
        Left eigenvectors
    s : numpy array
        Singular values
    v : numpy array
        Right eigenvectors