def rdl_decomposition_rev(T, k, norm='reversible', ncv=None, mu=None):
    r"""Compute the decomposition into left and right eigenvectors.

    Parameters
    ----------
    T : sparse matrix
        Transition matrix
    k : int
        Number of eigenvector/eigenvalue pairs
    norm: {'standard', 'reversible'}
        standard: (L'R) = Id, L[:,0] is a probability distribution,
            the stationary distribution mu of T. Right eigenvectors
            R have a 2-norm of 1.
        reversible: R and L are related via L=L[:,0]*R.
    ncv : int (optional)
        The number of Lanczos vectors generated, `ncv` must be greater than k;
        it is recommended that ncv > 2*k
    mu : (M,) ndarray, optional
        Stationary distribution of T

    Returns
    -------
    R : (M, M) ndarray
        The normalized ("unit length") right eigenvectors, such that the
        column ``R[:,i]`` is the right eigenvector corresponding to the eigenvalue
        ``w[i]``, ``dot(T,R[:,i])``=``w[i]*R[:,i]``
    D : (M, M) ndarray
        A diagonal matrix containing the eigenvalues, each repeated
        according to its multiplicity
    L : (M, M) ndarray
        The normalized (with respect to `R`) left eigenvectors, such that the
        row ``L[i, :]`` is the left eigenvector corresponding to the eigenvalue
        ``w[i]``, ``dot(L[i, :], T)``=``w[i]*L[i, :]``

    """
    if mu is None:
        mu = stationary_distribution(T)
    """ symmetrize T """
    smu = np.sqrt(mu)
    Dpi = diags(smu, 0)
    Dinv = diags(1.0/smu, 0)
    S = (Dpi.dot(T)).dot(Dinv)
    """Compute eigenvalues, eigenvecs using a solver for
    symmetric/hermititan eigenproblems"""
    val, eigvec = scipy.sparse.linalg.eigsh(S, k=k, ncv=ncv, which='LM',
                                            return_eigenvectors=True)
    """Sort eigenvalues and eigenvectors"""
    ind = np.argsort(np.abs(val))[::-1]
    val = val[ind]
    eigvec = eigvec[:, ind]

    """Diagonal matrix of eigenvalues"""
    D = np.diag(val)

    """Right and left eigenvectors"""
    R = eigvec / smu[:, np.newaxis]
    L = eigvec * smu[:, np.newaxis]

    """Ensure that R[:,0] is positive and unity"""
    tmp = R[0, 0]
    R[:, 0] = R[:, 0] / tmp

    """Ensure that L[:, 0] is probability vector"""
    L[:, 0] = L[:, 0] *  tmp

    if norm == 'reversible':
        return R, D, L.T
    elif norm == 'standard':
        """Standard l2-norm of right eigenvectors"""
        w = np.diag(np.dot(R.T, R))
        sw = np.sqrt(w)
        """Don't change normalization of eigenvectors for dominant eigenvalue"""
        sw[0] = 1.0

        R = R / sw[np.newaxis, :]
        L = L * sw[np.newaxis, :]
        return R, D, L.T
    else:
        raise ValueError("Keyword 'norm' has to be either 'standard' or 'reversible'")
def rdl_decomposition_rev ( T , k , norm = 'reversible' , ncv = None , mu = None ) : if mu is None : mu = stationary_distribution ( T ) """ symmetrize T """ smu = np . sqrt ( mu ) Dpi = diags ( smu , 0 ) Dinv = diags ( 1.0 / smu , 0 ) S = ( Dpi . dot ( T ) ) . dot ( Dinv ) val , eigvec = scipy . sparse . linalg . eigsh ( S , k = k , ncv = ncv , which = 'LM' , return_eigenvectors = True ) ind = np . argsort ( np . abs ( val ) ) [ : : - 1 ] val = val [ ind ] eigvec = eigvec [ : , ind ] """Diagonal matrix of eigenvalues""" D = np . diag ( val ) """Right and left eigenvectors""" R = eigvec / smu [ : , np . newaxis ] L = eigvec * smu [ : , np . newaxis ] """Ensure that R[:,0] is positive and unity""" tmp = R [ 0 , 0 ] R [ : , 0 ] = R [ : , 0 ] / tmp """Ensure that L[:, 0] is probability vector""" L [ : , 0 ] = L [ : , 0 ] * tmp if norm == 'reversible' : return R , D , L . T elif norm == 'standard' : w = np . diag ( np . dot ( R . T , R ) ) sw = np . sqrt ( w ) sw [ 0 ] = 1.0 R = R / sw [ np . newaxis , : ] L = L * sw [ np . newaxis , : ] return R , D , L . T else : raise ValueError ( "Keyword 'norm' has to be either 'standard' or 'reversible'" )
r"""Compute the decomposition into left and right eigenvectors.

    Parameters
    ----------
    T : sparse matrix
        Transition matrix
    k : int
        Number of eigenvector/eigenvalue pairs
    norm: {'standard', 'reversible'}
        standard: (L'R) = Id, L[:,0] is a probability distribution,
            the stationary distribution mu of T. Right eigenvectors
            R have a 2-norm of 1.
        reversible: R and L are related via L=L[:,0]*R.
    ncv : int (optional)
        The number of Lanczos vectors generated, `ncv` must be greater than k;
        it is recommended that ncv > 2*k
    mu : (M,) ndarray, optional
        Stationary distribution of T

    Returns
    -------
    R : (M, M) ndarray
        The normalized ("unit length") right eigenvectors, such that the
        column ``R[:,i]`` is the right eigenvector corresponding to the eigenvalue
        ``w[i]``, ``dot(T,R[:,i])``=``w[i]*R[:,i]``
    D : (M, M) ndarray
        A diagonal matrix containing the eigenvalues, each repeated
        according to its multiplicity
    L : (M, M) ndarray
        The normalized (with respect to `R`) left eigenvectors, such that the
        row ``L[i, :]`` is the left eigenvector corresponding to the eigenvalue
        ``w[i]``, ``dot(L[i, :], T)``=``w[i]*L[i, :]``