variance laplacians

laplacians.py

def _laplacian_dense2(graph, normalization=None, return_diag=False,
                      max_iters=100, tol=1e-5):
    """Compute different graph Laplacians
    """

    dd = np.sum(graph, axis=0)
    dd_idx = np.nonzero(dd)[0]
    dd_inv = dd
    dd_inv[dd_idx] = np.power(dd_inv[dd_idx], -1.)

    if normalization == 'bimarkov':
        # BIMARKOV computes the bimarkov normalization function p(x) for the
        # nonnegative, symemtric  kernel K(x,y) using an iterative scheme.  The
        # function p(x) is the unique function s.t.
        #
        #    diag(1./sqrt(p)*K*diag(1./sqrt(p))
        #
        # is both row and column stochastic.  Note that a bimarkov kernel
        # necessarily has L^2 norm 1.

        # initialize output
        lap = graph.copy()
        N = graph.shape[0]
        dd = np.ones(N) # bimarkov normalization function
        for i in range(max_iters):
            S = np.sum(lap, axis=1)
            err = max(abs(1-max(S)), abs(1-min(S)))
            if err < tol:
                break
            dd = S * dd
            lap = _laplacian_dense2(lap, normalization='sym_markov')
    if normalization == 'sym_markov':
        D = np.diag(np.sqrt(dd_inv))
        lap = D.dot(graph.dot(D))
        lap = (lap.T + lap)/2  # iron out numerical wrinkles
    if normalization == 'markov':
        lap = np.diag(dd_inv).dot(graph)
    if normalization == 'sym_beltrami':
        D = np.diag(dd_inv)
        K = D.dot(graph.dot(D))
        lap = _laplacian_dense2(K, normalization='sym_markov')
    if normalization == 'beltrami':
        D = np.diag(dd_inv)
        K = D.dot(graph.dot(D))
        lap = _laplacian_dense2(K, normalization='markov')
    if normalization == 'sym_fokker_planck':
        D = np.diag(np.sqrt(dd_inv))
        K = D.dot(graph.dot(D))
        lap = _laplacian_dense2(K, normalization='sym_markov')
    if normalization == 'fokker_planck':
        D = np.diag(np.sqrt(dd_inv))
        K = D.dot(graph.dot(D))
        lap = _laplacian_dense2(K, normalization='markov')
    if return_diag:
        return lap, dd
    return lap