added slepc eigensolver as an option for megaman

.travis.yml

@@ -13,7 +13,7 @@ env:
     # Directory where tests are run from
     - TEST_DIR=/tmp/megaman
     - CONDA_CHANNEL="conda-forge"
-    - CONDA_DEPS="pip nose coverage cython scikit-learn flann h5py"
+    - CONDA_DEPS="pip nose coverage cython scikit-learn flann h5py slepc4py"
     - PIP_DEPS="coveralls"
   matrix:
     - EXTRA_DEPS="pyflann pyamg"
@@ -38,7 +38,7 @@ install:
 
 script:
 - mkdir -p $TEST_DIR
-- cd $TEST_DIR && nosetests -v --with-coverage --cover-package=megaman megaman
+- cd $TEST_DIR && nosetests -v --nocapture --with-coverage --cover-package=megaman megaman
 
 after_success:
 - coveralls

README.md

@@ -67,6 +67,7 @@ Optional requirements include
 - [pyflann](http://www.cs.ubc.ca/research/flann/) which offers another method of computing distance matrices (this is bundled with the FLANN source code)
 - [nose](https://nose.readthedocs.org/) for running the unit tests
 - [h5py](http://www.h5py.org) for reading testing .mat files
+- [slepc4py] (http://slepc.upv.es/), for parallel eigendecomposition.
 
 These requirements can be installed on Linux and MacOSX using the following conda command:
 

megaman/embedding/tests/test_lle.py

@@ -66,6 +66,9 @@ def test_lle_simple_grid():
     distance_matrix = G.compute_adjacency_matrix()
     N = lle.barycenter_graph(distance_matrix, X).todense()
     reconstruction_error = np.linalg.norm(np.dot(N, X) - X, 'fro')
+    print('*************')
+    print('Reconstruction error for {}: {}'.format('numpy',reconstruction_error))
+    print('*************')
     assert(reconstruction_error < tol)
     for eigen_solver in EIGEN_SOLVERS:
         clf = lle.LocallyLinearEmbedding(n_components = n_components, geom = G,
@@ -74,6 +77,9 @@ def test_lle_simple_grid():
         assert(clf.embedding_.shape[1] == n_components)
         reconstruction_error = np.linalg.norm(
         np.dot(N, clf.embedding_) - clf.embedding_, 'fro') ** 2
+        print('*************')
+        print('Reconstruction error for {}: {}'.format(eigen_solver,reconstruction_error))
+        print('*************')
         assert(reconstruction_error < tol)
 
 def test_lle_manifold():

megaman/utils/eigendecomp.py

@@ -1,5 +1,5 @@
 # LICENSE: Simplified BSD https://github.com/mmp2/megaman/blob/master/LICENSE
-
+from __future__ import absolute_import
 import warnings
 import numpy as np
 from scipy import sparse
@@ -10,7 +10,7 @@
 from .validation import check_array
 
 
-EIGEN_SOLVERS = ['auto', 'dense', 'arpack', 'lobpcg']
+EIGEN_SOLVERS = ['auto', 'dense', 'arpack', 'lobpcg', 'slepc']
 BAD_EIGEN_SOLVERS = {}
 AMG_KWDS = ['strength', 'aggregate', 'smooth', 'max_levels', 'max_coarse']
 
@@ -23,6 +23,14 @@
     BAD_EIGEN_SOLVERS['amg'] = """The eigen_solver was set to 'amg',
     but pyamg is not available. Please either
     install pyamg or use another method."""
+
+try:
+    from . import slepctools as slepc
+    SLEPC_LOADED = True
+except ImportError:
+    #You need to install slepc4py with conda:
+    #conda install -c conda-forge slepc4py 
+    SLEPC_LOADED = False
 
 
 def check_eigen_solver(eigen_solver, solver_kwds, size=None, nvec=None):
@@ -68,6 +76,10 @@ def check_eigen_solver(eigen_solver, solver_kwds, size=None, nvec=None):
             else:
                 eigen_solver = 'dense'
                 solver_kwds = None
+        elif eigen_solver == 'slepc':
+            if not SLEPC_LOADED:
+                eigen_solver = 'arpack'
+                solver_kwds = None
 
     return eigen_solver, solver_kwds
 
@@ -136,7 +148,6 @@ def eigen_decomposition(G, n_components=8, eigen_solver='auto',
         if G.getformat() is not 'csr':
             G.tocsr()
     G = G.astype(np.float)
-
     # Check for symmetry
     is_symmetric = _is_symmetric(G)
 
@@ -222,6 +233,9 @@ def eigen_decomposition(G, n_components=8, eigen_solver='auto',
             diffusion_map = diffusion_map[:, ::-1] # reverse order the vectors
         lambdas = lambdas[:n_components]
         diffusion_map = diffusion_map[:, :n_components]
+    elif eigen_solver == 'slepc':
+        if is_symmetric:
+            lambdas, diffusion_map = slepc.get_eigenpairs(G,npairs=n_components,largest=largest)        
     return (lambdas, diffusion_map)
 
 
@@ -273,7 +287,7 @@ def null_space(M, k, k_skip=1, eigen_solver='arpack',
         for symmetric problems and
         http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.linalg.eig.html#scipy.linalg.eig
         for non symmetric problems.
-    arpack sovler key words: see
+    arpack solver key words: see
         http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.eigsh.html
         for symmetric problems and http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.eigs.html#scipy.sparse.linalg.eigs
         for non symmetric problems.
@@ -348,5 +362,8 @@ def null_space(M, k, k_skip=1, eigen_solver='arpack',
             eigen_values = eigen_values[index]
             eigen_vectors = eigen_vectors[:, index]
             return eigen_vectors[:, k_skip:k+1], np.sum(eigen_values[k_skip:k+1])
+    elif eigen_solver == 'slepc':
+        eigen_values, eigen_vectors = slepc.get_eigenpairs(M,npairs=k+k_skip,largest=False) 
+        return eigen_vectors[:, k_skip:], np.sum(eigen_values[k_skip:])
     else:
         raise ValueError("Unrecognized eigen_solver '%s'" % eigen_solver)

megaman/utils/slepctools.py

@@ -0,0 +1,188 @@
+from __future__ import absolute_import, division, print_function, unicode_literals
+from petsc4py import PETSc
+from slepc4py import SLEPc
+import numpy as np
+from scipy.sparse import csr_matrix
+Print = PETSc.Sys.Print
+__author__ = 'Murat Keceli'
+__date__ = 'February 13, 2017'
+"""
+You can install slepc4py and all its dependencies with conda:
+conda install -c conda-forge slepc4py
+More info on PETSc and SLEPc:
+https://www.mcs.anl.gov/petsc/
+http://slepc.upv.es/
+http://slepc.upv.es/documentation/slepc.pdf
+"""
+def get_petsc_matrix(Acsr,comm=PETSc.COMM_WORLD):
+    """
+    Given a scipy csr matrix returns PETSC AIJ matrix on a given mpi communicator.
+    Parameters
+    ----------
+    Acsr: Array like, scipy CSR formated sparse matrix
+    comm: MPI communicator
+    Returns
+    ----------
+    PETSc AIJ Mat
+    """
+    A = PETSc.Mat().createAIJ(Acsr.shape[0])
+    A.setUp()
+    rstart, rend = A.getOwnershipRange()
+    return A.createAIJ(size=Acsr.shape[0],
+                       csr=(Acsr.indptr[rstart:rend+1] - Acsr.indptr[rstart],
+                            Acsr.indices[Acsr.indptr[rstart]:Acsr.indptr[rend]],
+                            Acsr.data[Acsr.indptr[rstart]:Acsr.indptr[rend]]),
+                       comm=comm) 
+
+def get_numpy_array(xr):
+    """
+    Convert a distributed PETSc Vec into a sequential numpy array.
+    Parameters:
+    -----------
+    xr: PETSc Vec
+    Returns:
+    --------
+    Array
+    """
+    xr_size = xr.getSize()
+    seqx = PETSc.Vec()
+    seqx.createSeq(xr_size,comm=PETSc.COMM_SELF)
+    seqx.setFromOptions()
+   # seqx.set(0.0)
+    fromIS = PETSc.IS().createGeneral(range(xr_size),comm=PETSc.COMM_SELF)
+    toIS = PETSc.IS().createGeneral(range(xr_size),comm=PETSc.COMM_SELF)
+    sctr=PETSc.Scatter().create(xr,fromIS,seqx,toIS)
+    sctr.begin(xr,seqx,addv=PETSc.InsertMode.INSERT,mode=PETSc.ScatterMode.FORWARD)
+    sctr.end(xr,seqx,addv=PETSc.InsertMode.INSERT,mode=PETSc.ScatterMode.FORWARD)
+    return seqx.getArray()
+
+
+def get_eigenpairs(A, npairs=8, largest=True, comm=PETSc.COMM_WORLD):
+    """
+    Parameters:
+    -----------
+    A: scipy CSR matrix, or numpy array or PETSc mat 
+    npairs: int, number of eigenpairs required
+    largest: boolean, True (False) if largest (smallest) magnitude eigenvalues are requried
+    comm: MPI communicator
+    Returns:
+    --------
+    evals: 1D array of npairs, eigenvalues from largest to smallest
+    evecs: 2D array of (n, npairs) where n is size of A matrix
+    """
+    matrixtype = str(type(A))
+    if 'scipy' in matrixtype:
+        n = A.shape[0]
+        A = get_petsc_matrix(A, comm)
+    elif 'petsc' in matrixtype:
+        n = A.getSize()
+    elif 'numpy' in matrixtype:
+        A = csr_matrix(A)
+        n = A.shape[0]
+        A = get_petsc_matrix(A, comm)
+    else:
+        Print('Matrix type {} is not compatible. Use scipy CSR or PETSc AIJ type.'.format(type(A)))
+    mpisize = comm.size
+    Print('Matrix size: {}'.format(n))
+    eps = SLEPc.EPS()
+    eps.create()
+    eps.setOperators(A)
+    eps.setProblemType(SLEPc.EPS.ProblemType.HEP)
+    if largest:
+        eps.setWhichEigenpairs(SLEPc.EPS.Which.LARGEST_REAL)
+    else:
+        eps.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)
+    eps.setDimensions(nev=npairs,ncv=SLEPc.DECIDE,mpd=SLEPc.DECIDE)
+    eps.setTolerances(tol=1.e-14)
+    eps.setFromOptions() #Use command line options
+    eps.setUp()
+    eps.solve()
+    nconv = eps.getConverged()
+    vr = A.getVecLeft()
+    vi = A.getVecLeft()
+    if nconv < npairs:
+        Print("{} eigenvalues required, {} converged.".format(npairs,nconv))
+        npairs = nconv
+    evals = np.zeros(npairs)
+    evecs = np.zeros((n,npairs))
+    for i in range(npairs):
+        k = eps.getEigenpair(i,vr,vi)
+        if abs(k.imag) > 1.e-10:
+            Print("Imaginary eigenvalue: {} + {}j".format(k.real,k.imag))
+            Print("Error: {}".format(eps.computeError(i)))
+        if mpisize > 1:
+            evecs[:,i] = get_numpy_array(vr)
+        else:
+            evecs[:,i] = vr.getArray()
+        evals[i] = k.real
+    return evals, evecs  
+
+
+def get_null_space(A, npairs, comm=PETSc.COMM_WORLD):
+    """
+    Returns 'npairs' eigenpairs with eigenvalues close to 0. Uses shift-and-invert.
+    Parameters:
+    -----------
+    A: scipy CSR matrix, or numpy array or PETSc mat 
+    npairs : integer
+        Number of eigenvalues/vectors to return
+    comm: MPI communicator
+
+    Returns:
+    --------
+    evals: 1D array of npairs, eigenvalues from largest to smallest
+    evecs: 2D array of (n, npairs) where n is size of A matrix
+    """
+    matrixtype = str(type(A))
+    if 'scipy' in matrixtype:
+        n = A.shape[0]
+        A = get_petsc_matrix(A, comm)
+    elif 'petsc' in matrixtype:
+        n = A.getSize()
+    elif 'numpy' in matrixtype:
+        A = csr_matrix(A)
+        n = A.shape[0]
+        A = get_petsc_matrix(A, comm)
+    else:
+        Print('Matrix type {} is not compatible. Use scipy CSR or PETSc AIJ type.'.format(type(A)))
+    mpisize = comm.size
+    Print('Matrix size: {}'.format(n))
+    eps = SLEPc.EPS()
+    eps.create()
+    eps.setOperators(A)
+    eps.setProblemType(SLEPc.EPS.ProblemType.HEP)
+    eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_REAL)
+    eps.setTarget(0.)
+    eps.setDimensions(nev=npairs,ncv=SLEPc.DECIDE,mpd=SLEPc.DECIDE)
+    st=eps.getST()
+    ksp=st.getKSP()
+    pc=ksp.getPC()
+    pc.setType(PETSc.PC.Type.CHOLESKY)
+    pc.setFactorShift(shift_type=A.FactorShiftType.POSITIVE_DEFINITE)
+    pc.setFactorSolverPackage('mumps')
+    st.setType(SLEPc.ST.Type.SINVERT)
+    eps.setTolerances(tol=1.e-14)
+    eps.setFromOptions() #Use command line options
+    eps.setUp()
+    eps.solve()
+    nconv = eps.getConverged()
+    vr = A.getVecLeft()
+    vi = A.getVecLeft()
+    if nconv < npairs:
+        Print("{} eigenvalues required, {} converged.".format(npairs,nconv))
+        npairs = nconv
+    evals = np.zeros(npairs)
+    evecs = np.zeros((n,npairs))
+    for i in range(npairs):
+        k = eps.getEigenpair(i,vr,vi)
+        if abs(k.imag) > 1.e-10:
+            Print("Imaginary eigenvalue: {} + {}j".format(k.real,k.imag))
+            Print("Error: {}".format(eps.computeError(i)))
+        if mpisize > 1:
+            evecs[:,i] = get_numpy_array(vr)
+        else:
+            evecs[:,i] = vr.getArray()
+        evals[i] = k.real
+    return evals, evecs 
+
+

megaman/utils/tests/test_eigendecomp.py

@@ -12,7 +12,8 @@
                     'dense':{'turbo':True, 'type':1},
                     'arpack':{'mode':'normal', 'tol':0, 'maxiter':None},
                     'lobpcg':{'maxiter':20, 'tol':None},
-                    'amg':{'maxiter':20, 'tol':None,'aggregate':'standard'}}
+                    'amg':{'maxiter':20, 'tol':None,'aggregate':'standard'},
+                    'slepc':None}
 
 def _check_with_col_sign_flipping(A, B, tol=0.0):
     """ Check array A and B are equal with possible sign flipping on