Stat dist slepc (#22)

Use PETSc/SLEPc, if installed, to speed up computing the stationary distribution:

* Add SLEPc support for stat. dist

* Add test case matrices

* Add ground truth stationary distributions

* Add regression tests

* Add checks for stationary distribution

* Add tests for the stationary dist check

* Add release notes

Co-authored-by: Michal Klein <>
Co-authored-by: Marius Lange <marius.lange@t-online.de>

.gitignore

@@ -130,6 +130,7 @@ dmypy.json
 
 # Editors
 .idea
+*.DS_Store
 
 # docs
 docs/source/api/*rst

docs/source/index.rst

@@ -12,6 +12,7 @@
     api
     example
     acknowledgments
+    release_notes
     references
 
 .. |br| raw:: html

docs/source/release_notes.rst

@@ -0,0 +1,31 @@
+Release Notes
+=============
+
+.. role:: small
+
+Version 1.0
+-----------
+
+1.0.2 :small:`2021-03-26`
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. rubric:: Bugfixes
+
+- Fix not catching ``ArpackError`` when computing stationary distribution and mypy-related linting issues
+  `PR 21 <https://github.com/msmdev/pyGPCCA/pull/21>`_.
+
+.. rubric:: Improvements
+
+- Use PETSc/SLEPc, if installed, to speed up the computation of the stationary distribution
+  `PR 22 <https://github.com/msmdev/pyGPCCA/pull/22>`_.
+
+1.0.1 :small:`2021-02-01`
+~~~~~~~~~~~~~~~~~~~~~~~~~
+.. rubric:: General
+
+- Minor improvements/fixes in README and acknowledgements.
+
+1.0.0 :small:`2021-01-29`
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Initial release.

pygpcca/utils/_utils.py

@@ -17,7 +17,7 @@
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-from typing import List, Union
+from typing import List, Tuple, Union
 from functools import singledispatch
 
 from scipy.linalg import eig, lu_solve, lu_factor
@@ -217,13 +217,60 @@ def _sdd(P: np.ndarray) -> np.ndarray:
             mu = np.maximum(mu, 0.0)
             mu /= mu.sum()
 
+    # check whether this really is a stationary distribution
+    _is_stationary_distribution(P, mu)
+
     return mu
 
 
+def _eigs_slepc(P: spmatrix, k: int, which: "str" = "LR", tol: float = EPS) -> Tuple[np.ndarray, np.ndarray]:
+    from petsc4py import PETSc
+    from slepc4py import SLEPc
+
+    M = PETSc.Mat().create()
+    if not isspmatrix_csr(P):
+        P = csr_matrix(P)
+    M.createAIJ(size=P.shape, csr=(P.indptr, P.indices, P.data))
+
+    E = SLEPc.EPS()
+    E.create()
+    E.setOperators(M)
+    E.setDimensions(k)
+    E.setTolerances(tol=tol)
+    if which == "LR":
+        E.setWhichEigenpairs(E.Which.LARGEST_REAL)
+    elif which == "LM":
+        E.setWhichEigenpairs(E.Which.LARGEST_MAGNITUDE)
+    else:
+        raise NotImplementedError(f"`which={which}` is not implemented.")
+    E.solve()
+
+    nconv = E.getConverged()
+    if nconv < k:
+        raise ValueError(f"Requested `{k}` eigenvalues/vectors, but only `{nconv}` converged.")
+
+    xr, _ = M.getVecs()
+    xi, _ = M.getVecs()
+
+    eigenvalues, eigenvectors = [], []
+    for i in range(k):
+        # Get the i-th eigenvalue as computed by solve().
+        eigenvalues.append(E.getEigenpair(i, xr, xi))
+        if eigenvalues[-1].imag != 0.0:
+            eigenvectors.append([complex(r, i) for r, i in zip(xr.getArray(), xi.getArray())])
+        else:
+            eigenvectors.append(list(xr.getArray()))
+
+    return np.asarray(eigenvalues), np.asarray(eigenvectors).T
+
+
 @stationary_distribution.register(spmatrix)
 def _sds(P: spmatrix) -> np.ndarray:
     # get the top two eigenvalues and vecs so we can check for irreducibility
-    vals, vecs = eigs(P.transpose(), k=2, which="LR", ncv=None)
+    try:
+        vals, vecs = _eigs_slepc(P.T, k=2, which="LR")
+    except ImportError:
+        vals, vecs = eigs(P.T, k=2, which="LR", ncv=None)
 
     # check for irreducibility
     if np.allclose(vals, 1, rtol=1e2 * EPS, atol=1e2 * EPS):
@@ -244,9 +291,13 @@ def _sds(P: spmatrix) -> np.ndarray:
     if not (top_vec > -1e4 * EPS).all() and not (top_vec < 1e4 * EPS).all():
         raise ValueError("Top eigenvector has both positive and negative entries.")
     top_vec = np.abs(top_vec)
+    pi = top_vec / np.sum(top_vec)
+
+    # check whether this really is a stationary distribution
+    _is_stationary_distribution(P, pi)
 
     # normalize to 1 and return
-    return top_vec / np.sum(top_vec)
+    return pi
 
 
 def backward_iteration(A: np.ndarray, mu: float, x0: np.ndarray, tol: float = 1e-14, maxiter: int = 100) -> np.ndarray:
@@ -354,3 +405,24 @@ def stationary_distribution_from_eigenvector(P: np.ndarray) -> np.ndarray:
     nu = np.abs(L[:, 0])
 
     return nu / np.sum(nu)
+
+
+def _is_stationary_distribution(T: Union[np.ndarray, spmatrix], pi: np.ndarray) -> bool:
+
+    # check the shapes
+    if not T.shape[0] == T.shape[1] or not T.shape[0] == pi.shape[0]:
+        raise ValueError("Shape mismatch.")
+
+    # check for invariance
+    if not np.allclose(T.T.dot(pi), pi, rtol=1e4 * EPS, atol=1e4 * EPS):
+        raise ValueError("Stationary distribution is not invariant under the transition matrix.")
+
+    # check for positivity
+    if not (pi > -1e4 * EPS).all():
+        raise ValueError("Stationary distribution has negative elements.")
+
+    # check whether it sums to one
+    if not np.allclose(pi.sum(), 1, rtol=1e4 * EPS, atol=1e4 * EPS):
+        raise ValueError("Stationary distribution doe not sum to one.")
+
+    return True