Python scipy.sparse.linalg.lobpcg() Examples

def solveEigenspace(x, m = 2):
	k = numpy.random.uniform(.5, 1.0, (len(neighborhoods),m))

	for ix in range(20):
		t0 = time.time()
		res = lobpcg(x,k,largest=False,maxiter=50)

		k = res[1]

	return res 

def _get_fiedler_func(method):
    """Return a function that solves the Fiedler eigenvalue problem.
    """
    match = _tracemin_method.match(method)
    if match:
        method = match.group(1)
        def find_fiedler(L, x, normalized, tol):
            q = 2 if method == 'pcg' else min(4, L.shape[0] - 1)
            X = asmatrix(normal(size=(q, L.shape[0]))).T
            sigma, X = _tracemin_fiedler(L, X, normalized, tol, method)
            return sigma[0], X[:, 0]
    elif method == 'lanczos' or method == 'lobpcg':
        def find_fiedler(L, x, normalized, tol):
            L = csc_matrix(L, dtype=float)
            n = L.shape[0]
            if normalized:
                D = spdiags(1. / sqrt(L.diagonal()), [0], n, n, format='csc')
                L = D * L * D
            if method == 'lanczos' or n < 10:
                # Avoid LOBPCG when n < 10 due to
                # https://github.com/scipy/scipy/issues/3592
                # https://github.com/scipy/scipy/pull/3594
                sigma, X = eigsh(L, 2, which='SM', tol=tol,
                                 return_eigenvectors=True)
                return sigma[1], X[:, 1]
            else:
                X = asarray(asmatrix(x).T)
                M = spdiags(1. / L.diagonal(), [0], n, n)
                Y = ones(n)
                if normalized:
                    Y /= D.diagonal()
                sigma, X = lobpcg(L, X, M=M, Y=asmatrix(Y).T, tol=tol,
                                  maxiter=n, largest=False)
                return sigma[0], X[:, 0]
    else:
        raise nx.NetworkXError("unknown method '%s'." % method)

    return find_fiedler 

def _get_fiedler_func(method):
    """Returns a function that solves the Fiedler eigenvalue problem.
    """
    if method == "tracemin":  # old style keyword <v2.1
        method = "tracemin_pcg"
    if method in ("tracemin_pcg", "tracemin_chol", "tracemin_lu"):
        def find_fiedler(L, x, normalized, tol, seed):
            q = 1 if method == 'tracemin_pcg' else min(4, L.shape[0] - 1)
            X = asmatrix(seed.normal(size=(q, L.shape[0]))).T
            sigma, X = _tracemin_fiedler(L, X, normalized, tol, method)
            return sigma[0], X[:, 0]
    elif method == 'lanczos' or method == 'lobpcg':
        def find_fiedler(L, x, normalized, tol, seed):
            L = csc_matrix(L, dtype=float)
            n = L.shape[0]
            if normalized:
                D = spdiags(1. / sqrt(L.diagonal()), [0], n, n, format='csc')
                L = D * L * D
            if method == 'lanczos' or n < 10:
                # Avoid LOBPCG when n < 10 due to
                # https://github.com/scipy/scipy/issues/3592
                # https://github.com/scipy/scipy/pull/3594
                sigma, X = eigsh(L, 2, which='SM', tol=tol,
                                 return_eigenvectors=True)
                return sigma[1], X[:, 1]
            else:
                X = asarray(asmatrix(x).T)
                M = spdiags(1. / L.diagonal(), [0], n, n)
                Y = ones(n)
                if normalized:
                    Y /= D.diagonal()
                sigma, X = lobpcg(L, X, M=M, Y=asmatrix(Y).T, tol=tol,
                                  maxiter=n, largest=False)
                return sigma[0], X[:, 0]
    else:
        raise nx.NetworkXError("unknown method '%s'." % method)

    return find_fiedler 

def _get_fiedler_func(method):
    """Return a function that solves the Fiedler eigenvalue problem.
    """
    match = _tracemin_method.match(method)
    if match:
        method = match.group(1)

        def find_fiedler(L, x, normalized, tol):
            q = 2 if method == 'pcg' else min(4, L.shape[0] - 1)
            X = asmatrix(normal(size=(q, L.shape[0]))).T
            sigma, X = _tracemin_fiedler(L, X, normalized, tol, method)
            return sigma[0], X[:, 0]
    elif method == 'lanczos' or method == 'lobpcg':
        def find_fiedler(L, x, normalized, tol):
            L = csc_matrix(L, dtype=float)
            n = L.shape[0]
            if normalized:
                D = spdiags(1. / sqrt(L.diagonal()), [0], n, n, format='csc')
                L = D * L * D
            if method == 'lanczos' or n < 10:
                # Avoid LOBPCG when n < 10 due to
                # https://github.com/scipy/scipy/issues/3592
                # https://github.com/scipy/scipy/pull/3594
                sigma, X = eigsh(L, 2, which='SM', tol=tol,
                                 return_eigenvectors=True)
                return sigma[1], X[:, 1]
            else:
                X = asarray(asmatrix(x).T)
                M = spdiags(1. / L.diagonal(), [0], n, n)
                Y = ones(n)
                if normalized:
                    Y /= D.diagonal()
                sigma, X = lobpcg(L, X, M=M, Y=asmatrix(Y).T, tol=tol,
                                  maxiter=n, largest=False)
                return sigma[0], X[:, 0]
    else:
        raise nx.NetworkXError("unknown method '%s'." % method)

    return find_fiedler 

def solveSimply(AOrig,DOrig):
	AD = (AOrig+DOrig).tocsr()
	A = AOrig.tocsr()
	D = -1.0 / DOrig.diagonal()


	x = numpy.random.uniform(0, 1.0, (len(neighborhoods), 1))
	x = normalizeAndDemean(x)

	xSimple = x

	errors = []
	separation = []

	def sep(x):
		low = x[:len(neighborhoods)/2,0]
		high = x[len(neighborhoods)/2:,0]

		return abs((high.mean() - low.mean())/numpy.std(x))


	elapsed = 0.0
	elapsedSimple = 0.0
	for passIx in range(10):
		t0 = time.time()
		for ix in range(400):
			x = updateSingleLOBPCG(AD, x)
		elapsed += time.time() - t0

		t0 = time.time()
		for ix in range(100):
			xSimple = updateByAveragingAndStepping(A, D, xSimple)

		elapsedSimple += time.time() - t0

		print "LOBPCG-2 method: ", elapsed, " with error ", error(AD, x)
		print "Simple method:   ", elapsedSimple, " with error ", error(AD, xSimple)

	eigenT0 = time.time()
	eigenstyle = flip(solveEigenspace(AD,2)[1][:,1].reshape((len(neighborhoods),1)))
	eigenElapsed = time.time() - eigenT0
	print "Eigen: ", eigenElapsed, " with error ", error(AD, eigenstyle)


	plt.plot(flip(x),label="lobpcg-2")
	plt.plot(flip(xSimple),label="averaging")
	plt.plot(eigenstyle,label="eigenvector")
	plt.legend()
	plt.show() 

def check_eigen_solver(eigen_solver, solver_kwds, size=None, nvec=None):
    """Check that the selected eigensolver is valid

    Parameters
    ----------
    eigen_solver : string
        string value to validate
    size, nvec : int (optional)
        if both provided, use the specified problem size and number of vectors
        to determine the optimal method to use with eigen_solver='auto'

    Returns
    -------
    eigen_solver : string
        The eigen solver. This only differs from the input if
        eigen_solver == 'auto' and `size` is specified.
    """
    if eigen_solver in BAD_EIGEN_SOLVERS:
        raise ValueError(BAD_EIGEN_SOLVERS[eigen_solver])
    elif eigen_solver not in EIGEN_SOLVERS:
        raise ValueError("Unrecognized eigen_solver: '{0}'."
                         "Should be one of: {1}".format(eigen_solver,
                                                        EIGEN_SOLVERS))

    if size is not None and nvec is not None:
        # do some checks of the eigensolver
        if eigen_solver == 'lobpcg' and size < 5 * nvec + 1:
            warnings.warn("lobpcg does not perform well with small matrices or "
                          "with large numbers of vectors. Switching to 'dense'")
            eigen_solver = 'dense'
            solver_kwds = None

        elif eigen_solver == 'auto':
            if size > 200 and nvec < 10:
                if PYAMG_LOADED:
                    eigen_solver = 'amg'
                    solver_kwds = None
                else:
                    eigen_solver = 'arpack'
                    solver_kwds = None
            else:
                eigen_solver = 'dense'
                solver_kwds = None

    return eigen_solver, solver_kwds