Python scipy.sparse.linalg.eigs() Examples

def scaled_laplacian(W):
    '''
    Normalized graph Laplacian

    Parameters
    ----------
    W: np.ndarray, adjacency matrix,
       shape is (num_of_vertices, num_of_vertices)

    Returns
    ----------
    np.ndarray, shape is (num_of_vertices, num_of_vertices)

    '''

    num_of_vertices = W.shape[0]
    d = np.sum(W, axis=1)
    L = np.diag(d) - W
    for i in range(num_of_vertices):
        for j in range(num_of_vertices):
            if (d[i] > 0) and (d[j] > 0):
                L[i, j] = L[i, j] / np.sqrt(d[i] * d[j])
    # lambda_max \approx 2.0, the largest eigenvalues of L.
    lambda_max = eigs(L, k=1, which='LR')[0][0].real
    return 2 * L / lambda_max - np.identity(num_of_vertices) 

def test_power_method(self):
        """test for inverse power iteration method"""

        # solve eigenvalue problem in TT format
        evp.power_method(self.operator_tt, self.initial_tt, operator_gevp=self.operator_gevp)
        eigenvalue_tt, eigenvector_tt = evp.power_method(self.operator_tt, self.initial_tt)

        # solve eigenvalue problem in matrix format
        eigenvalue_mat, eigenvector_mat = splin.eigs(self.operator_mat, k=1)

        # compute relative error between exact and approximate eigenvalues
        rel_err_val = np.abs(eigenvalue_mat - eigenvalue_tt) / np.abs(eigenvalue_mat)

        # compute relative error between exact and approximate eigenvectors
        norm_1 = np.linalg.norm(eigenvector_mat + eigenvector_tt.matricize()[:, None])
        norm_2 = np.linalg.norm(eigenvector_mat - eigenvector_tt.matricize()[:, None])
        rel_err_vec = np.amin([norm_1, norm_2]) / np.linalg.norm(eigenvector_mat)

        # check if relative errors are smaller than tolerance
        self.assertLess(rel_err_val, self.tol_eigval)
        self.assertLess(rel_err_vec, self.tol_eigvec) 

def learn_embedding(self, graph=None, edge_f=None,
                        is_weighted=False, no_python=False):
        if not graph and not edge_f:
            raise Exception('graph/edge_f needed')
        if not graph:
            graph = graph_util.loadGraphFromEdgeListTxt(edge_f)
        graph = graph.to_undirected()
        t1 = time()
        L_sym = nx.normalized_laplacian_matrix(graph)

        w, v = lg.eigs(L_sym, k=self._d + 1, which='SM')
        idx = np.argsort(w) # sort eigenvalues
        w = w[idx]
        v = v[:, idx]
        t2 = time()
        self._X = v[:, 1:]

        p_d_p_t = np.dot(v, np.dot(np.diag(w), v.T))
        eig_err = np.linalg.norm(p_d_p_t - L_sym)
        print('Laplacian matrix recon. error (low rank): %f' % eig_err)
        return self._X.real, (t2 - t1) 

def get_pagerank_with_teleportation_from_transition_matrix(rw_transition, rw_transition_t, rho):
    number_of_nodes = rw_transition.shape[0]

    # Set up the random walk with teleportation matrix.
    non_teleportation = 1-rho
    mv = lambda l, v: non_teleportation*l.dot(v) + (rho/number_of_nodes)*np.ones_like(v)
    teleport = lambda vec: mv(rw_transition_t, vec)

    rw_transition_operator = spla.LinearOperator(rw_transition.shape, matvec=teleport, dtype=np.float64)

    # Calculate stationary distribution.
    try:
        eigenvalue, stationary_distribution = spla.eigs(rw_transition_operator,
                                                        k=1,
                                                        which='LM',
                                                        return_eigenvectors=True)
    except spla.ArpackNoConvergence as e:
        print("ARPACK has not converged.")
        eigenvalue = e.eigenvalues
        stationary_distribution = e.eigenvectors

    stationary_distribution = stationary_distribution.flatten().real/stationary_distribution.sum()

    return stationary_distribution 

def _method_2(data, num_pcs=None):
    """Compute OPCA when num_observations <= num_dimensions."""
    data = np.nan_to_num(data - nanmean(data, axis=0))
    T = data.shape[0]
    tmp = np.dot(data, data.T)
    corr_offset = np.zeros(tmp.shape)
    corr_offset[1:] = tmp[:-1]
    corr_offset[:-1] += tmp[1:]
    if num_pcs is None:
        eivals, eivects = eig(corr_offset)
    else:
        eivals, eivects = eigs(corr_offset, num_pcs, which='LR')
    eivals = np.real(eivals)
    eivects = np.real(eivects)
    idx = np.argsort(-eivals)  # sort the eigenvectors and eigenvalues
    eivals = old_div(eivals[idx], (2. * (T - 1)))
    eivects = eivects[:, idx]
    transformed_eivects = np.dot(data.T, eivects)
    for i in range(transformed_eivects.shape[1]):  # normalize the eigenvectors
        transformed_eivects[:, i] /= np.linalg.norm(transformed_eivects[:, i])
    return eivals, transformed_eivects, np.dot(data, transformed_eivects) 

def ldos0d_wf(h,e=0.0,delta=0.01,num_wf = 10,robust=False,tol=0):
  """Calculates the local density of states of a hamiltonian and
     writes it in file, using arpack"""
  if h.dimensionality==0:  # only for 0d
    intra = csc_matrix(h.intra) # matrix
  else: raise # not implemented...
  if robust: # go to the imaginary axis for stability
    eig,eigvec = slg.eigs(intra,k=int(num_wf),which="LM",
                        sigma=e+1j*delta,tol=tol) 
    eig = eig.real # real part only
  else: # Hermitic Hamiltonian
    eig,eigvec = slg.eigsh(intra,k=int(num_wf),which="LM",sigma=e,tol=tol) 
  d = np.array([0.0 for i in range(intra.shape[0])]) # initialize
  for (v,ie) in zip(eigvec.transpose(),eig): # loop over wavefunctions
    v2 = (np.conjugate(v)*v).real # square of wavefunction
    fac = delta/((e-ie)**2 + delta**2) # factor to create a delta
    d += fac*v2 # add contribution
#  d /= num_wf # normalize
  d /= np.pi # normalize
  d = spatial_dos(h,d) # resum if necessary
  g = h.geometry  # store geometry
  write_ldos(g.x,g.y,d,z=g.z) # write in file 

def ldos_arpack(intra,num_wf=10,robust=False,tol=0,e=0.0,delta=0.01):
  """Use arpack to calculate hte local density of states at a certain energy"""
  if robust: # go to the imaginary axis for stability
    eig,eigvec = slg.eigs(intra,k=int(num_wf),which="LM",
                        sigma=e+1j*delta,tol=tol) 
    eig = eig.real # real part only
  else: # Hermitic Hamiltonian
    eig,eigvec = slg.eigsh(intra,k=int(num_wf),which="LM",sigma=e,tol=tol) 
  d = np.array([0.0 for i in range(intra.shape[0])]) # initialize
  for (v,ie) in zip(eigvec.transpose(),eig): # loop over wavefunctions
    v2 = (np.conjugate(v)*v).real # square of wavefunction
    fac = delta/((e-ie)**2 + delta**2) # factor to create a delta
    d += fac*v2 # add contribution
#  d /= num_wf # normalize
  d /= np.pi # normalize
  return d 

def getLeadingEigenvector(A, normalized=True, lanczosVecs = 15, maxiter = 1000):
        """Compute normalized leading eigenvector of a given matrix A.

        @param A: sparse matrix for which leading eigenvector will be computed
        @param normalized: wheter or not to normalize. Default is C{True}
        @param lanczosVecs: number of Lanczos vectors to be used in the approximate
            calculation of eigenvectors and eigenvalues. This maps to the ncv parameter 
            of scipy's underlying function eigs. 
        @param maxiter: scaling factor for the number of iterations to be used in the 
            approximate calculation of eigenvectors and eigenvalues. The number of iterations 
            passed to scipy's underlying eigs function will be n*maxiter where n is the 
            number of rows/columns of the Laplacian matrix.
        """

        if _sparse.issparse(A) == False:
            raise TypeError("A must be a sparse matrix")

        # NOTE: ncv sets additional auxiliary eigenvectors that are computed
        # NOTE: in order to be more confident to find the one with the largest
        # NOTE: magnitude, see https://github.com/scipy/scipy/issues/4987
        w, pi = _sla.eigs( A, k=1, which="LM", ncv=lanczosVecs, maxiter=maxiter)
        pi = pi.reshape(pi.size,)
        if normalized:
            pi /= sum(pi)
        return pi 

def learn_embedding(self, graph=None, edge_f=None,
                        is_weighted=False, no_python=False):
        if not graph and not edge_f:
            raise Exception('graph/edge_f needed')
        if not graph:
            graph = graph_util.loadGraphFromEdgeListTxt(edge_f)
        graph = graph.to_undirected()
        t1 = time()
        L_sym = nx.normalized_laplacian_matrix(graph)
        
        try:
            w, v = lg.eigs(L_sym, k=self._d + 1, which='SM')
            t2 = time()
            self._X = v[:, 1:]

            p_d_p_t = np.dot(v, np.dot(np.diag(w), v.T))
            eig_err = np.linalg.norm(p_d_p_t - L_sym)
            print ('Laplacian matrix recon. error (low rank): %f' % eig_err)
            return self._X, (t2 - t1)
        except:
            print ('SVD did not converge. Assigning random emebdding')
            self._X = np.random.randn(L_sym.shape[0], self._d)
            t2 = time()
            return self._X, (t2 - t1) 

def check_stability(self):
        """
        Check that the weight matrix is stable

        :return:
        """
        if self.K < 100:
            eigs = np.linalg.eigvals(self.weight_model.W_effective)
            maxeig = np.amax(np.real(eigs))
        else:
            from scipy.sparse.linalg import eigs
            maxeig = eigs(self.weight_model.W_effective, k=1)[0]

        print("Max eigenvalue: ", maxeig)
        if maxeig < 1.0:
            return True
        else:
            return False 

def check_stability(self, verbose=False):
        """
        Check that the weight matrix is stable

        :return:
        """
        if self.K < 100:
            eigs = np.linalg.eigvals(self.weight_model.W_effective)
            maxeig = np.amax(np.real(eigs))
        else:
            from scipy.sparse.linalg import eigs
            maxeig = eigs(self.weight_model.W_effective, k=1)[0]

        if verbose:
            print("Max eigenvalue: ", maxeig)

        return maxeig < 1.0 

def get_scipy_eig_vec(self):
    """Computes scipy estimate of min eigenvalue for matrix M.

    Returns:
      eig_vec: Minimum absolute eigen value
      eig_val: Corresponding eigen vector
    """
    if not self.params['has_conv']:
      matrix_m = self.sess.run(self.dual_object.matrix_m)
      min_eig_vec_val, estimated_eigen_vector = eigs(matrix_m, k=1, which='SR',
                                                     tol=1E-4)
      min_eig_vec_val = np.reshape(np.real(min_eig_vec_val), [1, 1])
      return np.reshape(estimated_eigen_vector, [-1, 1]), min_eig_vec_val
    else:
      dim = self.dual_object.matrix_m_dimension
      input_vector = tf.placeholder(tf.float32, shape=(dim, 1))
      output_vector = self.dual_object.get_psd_product(input_vector)

      def np_vector_prod_fn(np_vector):
        np_vector = np.reshape(np_vector, [-1, 1])
        output_np_vector = self.sess.run(output_vector, feed_dict={input_vector:np_vector})
        return output_np_vector
      linear_operator = LinearOperator((dim, dim), matvec=np_vector_prod_fn)
      # Performing shift invert scipy operation when eig val estimate is available
      min_eig_vec_val, estimated_eigen_vector = eigs(linear_operator,
                                                     k=1, which='SR', tol=1E-4)
      min_eig_vec_val = np.reshape(np.real(min_eig_vec_val), [1, 1])
      return np.reshape(estimated_eigen_vector, [-1, 1]), min_eig_vec_val 

def compute_certificate(self, current_step, feed_dictionary):
    """ Function to compute the certificate based either current value
    or dual variables loaded from dual folder """
    feed_dict = feed_dictionary.copy()
    nu = feed_dict[self.nu]
    second_term = self.make_m_psd(nu, feed_dict)
    tf.logging.info('Nu after modifying: ' + str(second_term))
    feed_dict.update({self.nu: second_term})
    computed_certificate = self.sess.run(self.unconstrained_objective, feed_dict=feed_dict)

    tf.logging.info('Inner step: %d, current value of certificate: %f',
                    current_step, computed_certificate)

    # Sometimes due to either overflow or instability in inverses,
    # the returned certificate is large and negative -- keeping a check
    if LOWER_CERT_BOUND < computed_certificate < 0:
      _, min_eig_val_m = self.get_lanczos_eig(feed_dict=feed_dict)
      tf.logging.info('min eig val from lanczos: ' + str(min_eig_val_m))
      input_vector_m = tf.placeholder(tf.float32, shape=(self.matrix_m_dimension, 1))
      output_vector_m = self.get_psd_product(input_vector_m)

      def np_vector_prod_fn_m(np_vector):
        np_vector = np.reshape(np_vector, [-1, 1])
        feed_dict.update({input_vector_m:np_vector})
        output_np_vector = self.sess.run(output_vector_m, feed_dict=feed_dict)
        return output_np_vector
      linear_operator_m = LinearOperator((self.matrix_m_dimension,
                                          self.matrix_m_dimension),
                                         matvec=np_vector_prod_fn_m)
      # Performing shift invert scipy operation when eig val estimate is available
      min_eig_val_m_scipy, _ = eigs(linear_operator_m, k=1, which='SR', tol=TOL)

      tf.logging.info('min eig val m from scipy: ' + str(min_eig_val_m_scipy))

      if min_eig_val_m - TOL > 0:
        tf.logging.info('Found certificate of robustness!')
        return True

    return False 

def _make_diffusion_coords(self, L):
        evals, evecs = spsl.eigs(L, k=(self.n_evecs+1), which='LR')
        ix = evals.argsort()[::-1][1:]
        evals = np.real(evals[ix])
        evecs = np.real(evecs[:, ix])
        dmap = np.dot(evecs, np.diag(np.sqrt(-1. / evals)))
        return dmap, evecs, evals 

def test_tf_lanczos_smallest_eigval(self):
    tf_num_iter = tf.placeholder(dtype=tf.int32, shape=())
    tf_matrix = tf.placeholder(dtype=tf.float32)
    def _vector_prod_fn(x):
      return tf.matmul(tf_matrix, tf.reshape(x, [-1, 1]))

    min_eigen_fn = autograph.to_graph(utils.tf_lanczos_smallest_eigval)
    init_vec_ph = tf.placeholder(shape=(MATRIX_DIMENTION, 1), dtype=tf.float32)
    tf_eigval, tf_eigvec = min_eigen_fn(
        _vector_prod_fn, MATRIX_DIMENTION, init_vec_ph, tf_num_iter, dtype=tf.float32)
    eigvec = np.zeros((MATRIX_DIMENTION, 1), dtype=np.float32)

    with self.test_session() as sess:
      # run this test for a few random matrices
      for _ in range(NUM_RANDOM_MATRICES):
        matrix = np.random.random((MATRIX_DIMENTION, MATRIX_DIMENTION))
        matrix = matrix + matrix.T  # symmetrizing matrix
        eigval, eigvec = sess.run(
            [tf_eigval, tf_eigvec],
            feed_dict={tf_num_iter: NUM_LZS_ITERATIONS, tf_matrix: matrix, init_vec_ph: eigvec})

        scipy_min_eigval, scipy_min_eigvec = eigs(
            matrix, k=1, which='SR')
        scipy_min_eigval = np.real(scipy_min_eigval)
        scipy_min_eigvec = np.real(scipy_min_eigvec)
        scipy_min_eigvec = scipy_min_eigvec / np.linalg.norm(scipy_min_eigvec)

        np.testing.assert_almost_equal(eigval, scipy_min_eigval, decimal=3)
        np.testing.assert_almost_equal(np.linalg.norm(eigvec), 1.0, decimal=3)
        abs_dot_prod = abs(np.dot(eigvec.flatten(), scipy_min_eigvec.flatten()))
        np.testing.assert_almost_equal(abs_dot_prod, 1.0, decimal=3) 

def solver_eigs(A, Neigs, guess_value=0, guess_vector=None, timing=False):
    # solves for the eigenmodes of A

    if timing:
        start = time()
    (values, vectors) = spl.eigs(A, k=Neigs, sigma=guess_value, v0=guess_vector, which='LM')
    if timing:
        end = time()
        print('Elapsed time for eigs() is %.4f secs' % (end - start))
    return (values, vectors) 

def estimate(self, preference_matrix):
        super()._check_matrix(preference_matrix)
        width = preference_matrix.shape[0]

        _, vectors = eigs(preference_matrix, k=(width-2), sigma=width, which='LM', v0=np.ones(width))

        real_vector = np.real([vec for vec in np.transpose(vectors) if not np.all(np.imag(vec))][:1])
        sum_vector = np.sum(real_vector)

        self._evaluate_consistency(preference_matrix)

        return np.around(real_vector, decimals=3)[0] / sum_vector 

def keigs(T, k, P, take_diagonal=0):
        """ return k largest magnitude eigenvalues for the matrix T.
        :param T: Matrix to find eigen values/vectors of
        :param k: number of eigen values/vectors to return
        :param P: in the case of symmetric normalizations,
                  this is the NxN diagonal matrix which relates the nonsymmetric
                  version to the symmetric form via conjugation
        :param take_diagonal: if 1, returns the eigenvalues as a vector rather than as a
                              diagonal matrix.
        """
        D, V = eigs(T, k, tol=1e-4, maxiter=1000)
        D = np.real(D)
        V = np.real(V)
        inds = np.argsort(D)[::-1]
        D = D[inds]
        V = V[:, inds]
        if P is not None:
            V = P.dot(V)

        # Normalize
        for i in range(V.shape[1]):
            V[:, i] = V[:, i] / norm(V[:, i])
        V = np.round(V, 10)

        if take_diagonal == 0:
            D = np.diag(D)

        return V, D 

def eigenvector(network):
    """Calculates eigenvector centralities of nodes.

    Parameters
    ----------
    network

    Returns
    -------
    dict

    """
    assert isinstance(network, Network), \
        "network must be an instance of Network"
    adj_mat = network.adjacency_matrix(weighted=False, transposed=True).asfptype()
    # calculate leading eigenvector of A
    _, v = sla.eigs(adj_mat, k=1, which="LM")

    v = v.reshape(v.size, )

    eigen_vec_cent = dict(zip(network.nodes, map(_np.abs, v)))
    S = sum(eigen_vec_cent.values())
    for v in eigen_vec_cent:
        eigen_vec_cent[v] /= S

    return eigen_vec_cent 

def algebraic_connectivity(network, lanczos_vectors=15, maxiter=20):
    """

    Parameters
    ----------
    network: Network
    lanczos_vectors: int
        number of Lanczos vectors to be used in the approximate calculation of
        eigenvectors and eigenvalues. This maps to the ncv parameter of scipy's underlying
        function eigs.
    maxiter: int
        scaling factor for the number of iterations to be used in the approximate
        calculation of eigenvectors and eigenvalues. The number of iterations passed to
        scipy's underlying eigs function will be n*maxiter where n is the number of
        rows/columns of the Laplacian matrix.

    Returns
    -------

    """
    assert isinstance(network, Network), \
        "network must be an instance of Network"
    Log.add('Calculating algebraic connectivity ... ', Severity.INFO)

    lapl_mat = network.laplacian_matrix()
    # NOTE: ncv sets additional auxiliary eigenvectors that are computed
    # NOTE: in order to be more confident to find the one with the largest
    # NOTE: magnitude, see https://github.com/scipy/scipy/issues/4987
    w = sla.eigs(lapl_mat, which="SM", k=2, ncv=lanczos_vectors,
                 return_eigenvectors=False, maxiter=maxiter)
    eigen_values_sorted = _np.sort(_np.absolute(w))

    Log.add('finished.', Severity.INFO)

    # TODO: result is unstable, it looks like it depends on a "warm start"
    # (i.e. run after other eigen velue calculations) see test_algebraic_connectivity
    # problems with order k=3

    return _np.abs(eigen_values_sorted[1]) 

def leading_eigenvector(A, normalized=True, lanczos_vecs=None, maxiter=None):
        """Compute normalized leading eigenvector of a given matrix A.

        Parameters
        ----------
        A:
            sparse matrix for which leading eigenvector will be computed
        normalized: bool
            whether or not to normalize, default is True
        lanczos_vecs: int
            number of Lanczos vectors to be used in the approximate
            calculation of eigenvectors and eigenvalues. This maps to the ncv parameter
            of scipy's underlying function eigs.
        maxiter: int
            scaling factor for the number of iterations to be used in the
            approximate calculation of eigenvectors and eigenvalues.

        Returns
        -------

        """
        if not _sparse.issparse(A):  # pragma: no cover
            raise TypeError("A must be a sparse matrix")

        # NOTE: ncv sets additional auxiliary eigenvectors that are computed
        # NOTE: in order to be more confident to find the one with the largest
        # NOTE: magnitude, see https://github.com/scipy/scipy/issues/4987
        if lanczos_vecs == None or maxiter == None:
            w, pi = _sla.eigs(A, k=1, which="LM")
        else:
            w, pi = _sla.eigs(A, k=1, which="LM", ncv=lanczos_vecs, maxiter=maxiter)
        pi = pi.reshape(pi.size, )
        if normalized:
            pi /= sum(pi)
        return pi 

def eigs(A, *args, **kwargs):
    return _eigs(A, *args, **kwargs) 

def test_als_eigh(self):
        """test for ALS with solver='eigh'"""

        # make problem symmetric
        self.operator_tt = 0.5 * (self.operator_tt + self.operator_tt.transpose())
        self.operator_mat = 0.5 * (self.operator_mat + self.operator_mat.transpose())

        # solve eigenvalue problem in TT format (number_ev=self.number_ev)
        eigenvalues_tt, eigenvectors_tt = evp.als(self.operator_tt, self.initial_tt, number_ev=self.number_ev,
                                                  repeats=10, solver='eigh')

        # solve eigenvalue problem in matrix format
        eigenvalues_mat, eigenvectors_mat = splin.eigs(self.operator_mat, k=self.number_ev)
        eigenvalues_mat = np.real(eigenvalues_mat)
        eigenvectors_mat = np.real(eigenvectors_mat)
        idx = eigenvalues_mat.argsort()[::-1]
        eigenvalues_mat = eigenvalues_mat[idx]
        eigenvectors_mat = eigenvectors_mat[:, idx]

        # compute relative error between exact and approximate eigenvalues
        rel_err_val = []
        for i in range(self.number_ev):
            rel_err_val.append(np.abs(eigenvalues_mat[i] - eigenvalues_tt[i]) / eigenvalues_mat[i])

        # compute relative error between exact and approximate eigenvectors
        rel_err_vec = []
        for i in range(self.number_ev):
            norm_1 = np.linalg.norm(eigenvectors_mat[:, i] + eigenvectors_tt[i].matricize())
            norm_2 = np.linalg.norm(eigenvectors_mat[:, i] - eigenvectors_tt[i].matricize())
            rel_err_vec.append(np.amin([norm_1, norm_2]) / np.linalg.norm(eigenvectors_mat[:, i]))

        # check if relative errors are smaller than tolerance
        for i in range(self.number_ev):
            self.assertLess(rel_err_val[i], self.tol_eigval)
            self.assertLess(rel_err_vec[i], self.tol_eigvec) 

def init_matrices(self):
        self.W_in = (2.0 * np.random.rand(self.N_x, self.num_dim) - 1.0) / (2.0 * self.scaleW_in)

        converged = False

        i = 0

        # repeat because could not converge to find eigenvalues
        while (not converged):
            i += 1

            # generate sparse, uniformly distributed weights
            self.W = sparse.rand(self.N_x, self.N_x, density=self.connect).todense()

            # ensure that the non-zero values are uniformly distributed
            self.W[np.where(self.W > 0)] -= 0.5

            try:
                # get the largest eigenvalue
                eig, _ = slinalg.eigs(self.W, k=1, which='LM')
                converged = True
            except:
                print('not converged ', i)
                continue

        # adjust the spectral radius
        self.W /= np.abs(eig) / self.rho 

def check_stability(self):
        """
        Check that the weight matrix is stable

        :return:
        """
        # Compute the effective weight matrix
        W_eff = self.weights.sum(axis=2)
        eigs = np.linalg.eigvals(W_eff)
        maxeig = np.amax(np.real(eigs))
        # print "Max eigenvalue: ", maxeig
        if maxeig < 1.0:
            return True
        else:
            return False 

def getSlowDownFactor(self, k=2, lanczosVecs = 15, maxiter = 1000):    
        """
        Returns a factor S that indicates how much slower (S>1) or faster (S<1)
        a diffusion process evolves in a k-order model of the path statistics
        compared to what is expected based on a first-order model. This value captures 
        the effect of order correlations of length k on a diffusion process which evolves 
        based on the observed paths.
        """

        assert k>1, 'Slow-down factor can only be calculated for orders larger than one'
    
        #NOTE to myself: most of the time goes for construction of the 2nd order
        #NOTE            null graph, then for the 2nd order null transition matrix
    
        gk = HigherOrderNetwork(self, k=k)
        gkn = HigherOrderNetwork(self, k=k, nullModel = True)
    
        Log.add('Calculating slow down factor ... ', Severity.INFO)

        # Build transition matrices
        Tk = gk.getTransitionMatrix()
        Tkn = gkn.getTransitionMatrix()
    
        # Compute eigenvector sequences
        # NOTE: ncv=13 sets additional auxiliary eigenvectors that are computed
        # NOTE: in order to be more confident to find the one with the largest
        # NOTE: magnitude, see
        # NOTE: https://github.com/scipy/scipy/issues/4987
        w2 = _sla.eigs(Tk, which="LM", k=2, ncv=lanczosVecs, return_eigenvectors=False, maxiter=maxiter)
        evals2_sorted = _np.sort(-_np.absolute(w2))

        w2n = _sla.eigs(Tkn, which="LM", k=2, ncv=lanczosVecs, return_eigenvectors=False, maxiter=maxiter)
        evals2n_sorted = _np.sort(-_np.absolute(w2n))

        Log.add('finished.', Severity.INFO)
    
        return _np.log(_np.abs(evals2n_sorted[1]))/_np.log(_np.abs(evals2_sorted[1])) 

def getEigenValueGap(self, includeSubPaths=True, lanczosVecs = 15, maxiter = 20):
        """
        Returns the eigenvalue gap of the transition matrix.

        @param includeSubPaths: whether or not to include subpath statistics in the 
            calculation of transition probabilities.
        """
    
        #NOTE to myself: most of the time goes for construction of the 2nd order
        #NOTE            null graph, then for the 2nd order null transition matrix   
    
        Log.add('Calculating eigenvalue gap ... ', Severity.INFO)

        # Build transition matrices
        T = self.getTransitionMatrix(includeSubPaths)
    
        # Compute the two largest eigenvalues
        # NOTE: ncv sets additional auxiliary eigenvectors that are computed
        # NOTE: in order to be more confident to actually find the one with the largest
        # NOTE: magnitude, see https://github.com/scipy/scipy/issues/4987
        w2 = _sla.eigs(T, which="LM", k=2, ncv=lanczosVecs, return_eigenvectors=False, maxiter = maxiter)
        evals2_sorted = _np.sort(-_np.absolute(w2))
        
        Log.add('finished.', Severity.INFO)
    
        return _np.abs(evals2_sorted[1]) 

def est_CompGraph_norm(K, tol=1e-3, try_fast_norm=True):
    """Estimates operator norm for L = ||K||.

    Parameters
    ----------
    tol : float
        Accuracy of estimate if not trying for upper bound.
    try_fast_norm : bool
        Whether to try for a fast upper bound.

    Returns
    -------
    float
        Estimate of ||K||.
    """
    if try_fast_norm:
        output_mags = [NotImplemented]
        K.norm_bound(output_mags)
        if NotImplemented not in output_mags:
            return output_mags[0]

    input_data = np.zeros(K.input_size)
    output_data = np.zeros(K.output_size)

    def KtK(x):
        K.forward(x, output_data)
        K.adjoint(output_data, input_data)
        return input_data

    # Define linear operator
    A = LinearOperator((K.input_size, K.input_size),
                       KtK, KtK)

    Knorm = np.sqrt(eigs(A, k=1, M=None, sigma=None, which='LM', tol=tol)[0].real)
    return np.float(Knorm) 

def init_matrices(self):
		self.W_in = (2.0*np.random.rand(self.N_x,self.num_dim)-1.0)/(2.0*self.scaleW_in)

		converged = False

		i =0 

		# repeat because could not converge to find eigenvalues 
		while(not converged):
			i+=1

			# generate sparse, uniformly distributed weights
			self.W = sparse.rand(self.N_x,self.N_x,density=self.connect).todense()

			# ensure that the non-zero values are uniformly distributed 
			self.W[np.where(self.W>0)] -= 0.5

			try:
				# get the largest eigenvalue 
				eig, _ = slinalg.eigs(self.W,k=1,which='LM')
				converged = True
			except: 
				print('not converged ',i)
				continue

		# adjust the spectral radius
		self.W /= np.abs(eig)/self.rho 

def eigs(T, k=10, eps=1e-3, perc=None, random_state=None, v0=None):
    if random_state is not None:
        np.random.seed(random_state)
        v0 = np.random.rand(min(T.shape))
    try:
        # find k eigs with largest real part, and sort in descending order of eigenvals
        eigvals, eigvecs = linalg.eigs(T.T, k=k, which="LR", v0=v0)
        p = np.argsort(eigvals)[::-1]
        eigvals = eigvals.real[p]
        eigvecs = eigvecs.real[:, p]

        # select eigenvectors with eigenvalue of 1 - eps.
        idx = eigvals >= 1 - eps
        eigvals = eigvals[idx]
        eigvecs = np.absolute(eigvecs[:, idx])

        if perc is not None:
            lbs, ubs = np.percentile(eigvecs, perc, axis=0)
            eigvecs[eigvecs < lbs] = 0
            eigvecs = np.clip(eigvecs, 0, ubs)
            eigvecs /= eigvecs.max(0)

    except:
        eigvals, eigvecs = np.empty(0), np.zeros(shape=(T.shape[0], 0))

    return eigvals, eigvecs