Error in calculating eigenvalues with Schur decomposition

[hattiechung]

I get the following error at the step of identifying terminal states. This is after calculating the velocity vectors.

cr.tl.terminal_states(ad_inf, cluster_key="annotation", weight_connectivities=0.2) #, n_states=8)

100%
2249/2249 [00:01<00:00, 2009.24cell/s]
100%
2249/2249 [00:01<00:00, 1999.66cell/s]
Mat Object: 1 MPI process
  type: seqdense
9.9999999999999489e-01 1.4114090614410588e-03 1.1568412248044814e-02 -2.9786989288699676e-03 8.9986905021165534e-03 
0.0000000000000000e+00 9.9780745509596869e-01 1.9245281401477242e-03 -9.4907745539946108e-04 -4.2135097291487239e-03 
0.0000000000000000e+00 0.0000000000000000e+00 9.9160368866063386e-01 -4.3116686139400328e-03 1.7720437331391087e-03 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.8229576348708580e-01 2.0405866724520394e-03 
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.6517123047558040e-01 
---------------------------------------------------------------------------
Error                                     Traceback (most recent call last)
/tmp/ipykernel_420933/3529991258.py in <module>
----> 1 cr.tl.terminal_states(ad_inf, cluster_key="annotation", weight_connectivities=0.2) #, n_states=8)

~/.local/lib/python3.7/site-packages/cellrank/tl/_utils.py in wrapper(wrapped, instance, args, kwargs)
   1641                 category=DeprecationWarning,
   1642             )
-> 1643         return wrapped(*args, **kwargs)
   1644 
   1645     return wrapper

~/.local/lib/python3.7/site-packages/cellrank/tl/_init_term_states.py in terminal_states(adata, estimator, mode, n_states, cluster_key, key, force_recompute, show_plots, copy, return_estimator, fit_kwargs, **kwargs)
    276         return_estimator=return_estimator,
    277         fit_kwargs=fit_kwargs,
--> 278         **kwargs,
    279     )

~/.local/lib/python3.7/site-packages/cellrank/tl/_init_term_states.py in _initial_terminal(adata, estimator, backward, mode, backward_mode, n_states, cluster_key, key, force_recompute, show_plots, copy, return_estimator, fit_kwargs, **kwargs)
    176         cluster_key=cluster_key,
    177         compute_absorption_probabilities=False,
--> 178         **fit_kwargs,
    179     )
    180 

~/.local/lib/python3.7/site-packages/cellrank/tl/_init_term_states.py in _fit(estim, n_lineages, keys, cluster_key, compute_absorption_probabilities, **kwargs)
     97 
     98         if n_lineages > 1:
---> 99             estim.compute_schur(n_lineages, method=kwargs.pop("method", "krylov"))
    100 
    101         try:

~/.local/lib/python3.7/site-packages/cellrank/tl/estimators/mixins/decomposition/_schur.py in compute_schur(self, n_components, initial_distribution, method, which, alpha)
    163 
    164         try:
--> 165             self._gpcca._do_schur_helper(n_components)
    166         except ValueError as e:
    167             if "will split complex conjugate eigenvalues" not in str(e):

~/.local/lib/python3.7/site-packages/pygpcca/_gpcca.py in _do_schur_helper(self, m)
    852         else:
    853             self._p_X, self._p_R, self._p_eigenvalues = _do_schur(
--> 854                 self._P, eta=self._eta, m=m, z=self._z, method=self._method
    855             )
    856 

~/.local/lib/python3.7/site-packages/pygpcca/_gpcca.py in _do_schur(P, eta, m, z, method, tol_krylov)
    252 
    253     # Make a Schur decomposition of P_bar and sort the Schur vectors (and form).
--> 254     R, Q, eigenvalues = sorted_schur(P_bar, m, z, method, tol_krylov=tol_krylov)  # Pbar!!!
    255 
    256     # Orthonormalize the sorted Schur vectors Q via modified Gram-Schmidt-orthonormalization,

~/.local/lib/python3.7/site-packages/pygpcca/_sorted_schur.py in sorted_schur(P, m, z, method, tol_krylov)
    381         R, Q, eigenvalues = sorted_brandts_schur(P=P, k=m, z=z)
    382     elif method == "krylov":
--> 383         R, Q, eigenvalues, _ = sorted_krylov_schur(P=P, k=m, z=z, tol=tol_krylov)
    384     else:
    385         raise ValueError(f"Unknown method `{method!r}`.")

~/.local/lib/python3.7/site-packages/pygpcca/_sorted_schur.py in sorted_krylov_schur(P, k, z, tol)
    286         eigenvalues.append(eigenval)
    287         # Computes the error (based on the residual norm) associated with the i-th computed eigenpair.
--> 288         eigenval_error = E.computeError(i)
    289         eigenvalues_error.append(eigenval_error)
    290 

SLEPc/EPS.pyx in slepc4py.SLEPc.EPS.computeError()

Error: error code 58
[0] EPSComputeError() at /tmp/pip-install-qjixudoz/slepc_b1c1741ad605407bac7221c6a94cb806/src/eps/interface/epssolve.c:710
[0] EPSGetEigenpair() at /tmp/pip-install-qjixudoz/slepc_b1c1741ad605407bac7221c6a94cb806/src/eps/interface/epssolve.c:406
[0] EPSGetEigenvector() at /tmp/pip-install-qjixudoz/slepc_b1c1741ad605407bac7221c6a94cb806/src/eps/interface/epssolve.c:502
[0] EPSComputeVectors() at /tmp/pip-install-qjixudoz/slepc_b1c1741ad605407bac7221c6a94cb806/src/eps/interface/epssolve.c:22
[0] EPSComputeVectors_Schur() at /tmp/pip-install-qjixudoz/slepc_b1c1741ad605407bac7221c6a94cb806/src/eps/interface/epsdefault.c:121
[0] DSVectors() at /tmp/pip-install-qjixudoz/slepc_b1c1741ad605407bac7221c6a94cb806/src/sys/classes/ds/interface/dsops.c:942
[0] DSVectors_NHEP() at /tmp/pip-install-qjixudoz/slepc_b1c1741ad605407bac7221c6a94cb806/src/sys/classes/ds/impls/nhep/dsnhep.c:239
[0] DSVectors_NHEP_Eigen_All() at /tmp/pip-install-qjixudoz/slepc_b1c1741ad605407bac7221c6a94cb806/src/sys/classes/ds/impls/nhep/dsnhep.c:179
[0] MatDenseGetArrayRead() at /tmp/pip-install-qjixudoz/petsc_4a14ae4aca574b31b612d08e29db0d32/src/mat/impls/dense/seq/dense.c:2264
[0] MatDenseGetArray_SeqDense() at /tmp/pip-install-qjixudoz/petsc_4a14ae4aca574b31b612d08e29db0d32/src/mat/impls/dense/seq/dense.c:2136
[0] Operation done in wrong order
[0] Need to call MatDenseRestoreSubMatrix() first

Versions:

cellrank==1.5.1 scanpy==1.9.1 anndata==0.8.0 python==3.7.12 pygpcca=1.0.3 scvelo==0.2.4 scipy==1.7.3

[michalk8]

Hi @hattiechung , I believe this comes with the new PETSc/SLEPc version (>=3.18), see related issue https://github.com/msmdev/pyGPCCA/pull/42. Could you please print the versions of the following packages as, e.g.:

conda list | egrep "petsc|slepc"

or

import petsc4py
import slepc4py

print(petsc4py.__version__, slepc4py.__version__, petsc)

[hattiechung]

Ah yes, both petsc4py and slepc4py versions are 3.18.0.

[hattiechung]

@michalk8 should I install the dev branch of cellrank?

[hattiechung]

I reinstalled pyGPCCA to pull the updated version using the following and now the issue is resolved! pip install git+https://github.com/msmdev/pyGPCCA.git