r"""Compute eigenvectors of reversible transition matrix.

    Parameters
    ----------
    T : (M, M) scipy.sparse matrix
        Transition matrix (stochastic matrix)
    k : int
        Number of eigenvalues to compute
    right : bool, optional
        If True compute right eigenvectors, left eigenvectors otherwise
    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
    -------
    eigvec : (M, k) ndarray
        k-eigenvectors of T

def eigenvectors_rev(T, k, right=True, ncv=None, mu=None):
    r"""Compute eigenvectors of reversible transition matrix.

    Parameters
    ----------
    T : (M, M) scipy.sparse matrix
        Transition matrix (stochastic matrix)
    k : int
        Number of eigenvalues to compute
    right : bool, optional
        If True compute right eigenvectors, left eigenvectors otherwise
    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
    -------
    eigvec : (M, k) ndarray
        k-eigenvectors of T

    """
    if mu is None:
        mu = stationary_distribution(T)
    """ symmetrize T """
    smu = np.sqrt(mu)
    D = diags(smu, 0)
    Dinv = diags(1.0/smu, 0)
    S = (D.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 eigenvectors"""
    ind = np.argsort(np.abs(val))[::-1]
    eigvec = eigvec[:, ind]
    if right:
        return eigvec / smu[:, np.newaxis]
    else:
        return eigvec * smu[:, np.newaxis]