GPCCA for matrix with absorbing states

[Marius1311]

Hi @msmdev , I have a large, reducible matrix with 5 recurrent and a number of transient states. I run GPCCA on it, both with method brandts as well as krylov_schur. In both cases, I get the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-96-a7b1ed0126a1> in <module>
----> 1 gp.compute_metastable_states(n_states=5, method='brandts')
      2 gp.plot_metastable_states()

~/Projects/cellrank/cellrank/tools/estimators/_gppca.py in compute_metastable_states(self, n_states, initial_distribution, use_min_chi, method, which, n_cells, cluster_key, en_cutoff, p_thresh)
    371             start = logg.info("Computing metastable states")
    372 
--> 373             gpcca = gpcca.optimize(m=n_states)
    374 
    375             # when `n_cells!=None` and the overlap is high, we're skipping some metastable states

~/Projects/msmtools/msmtools/analysis/dense/gpcca.py in optimize(self, m, return_extra)
   1335         # Calculate Schur matrix R and Schur vector matrix X, if not adequately given.
   1336 
-> 1337         self._do_schur_helper(max(m_list))
   1338 
   1339         # Initialize lists to collect results.

~/Projects/msmtools/msmtools/analysis/dense/gpcca.py in _do_schur_helper(self, m)
   1137                     self.X, self.R = _do_schur(self.P, self.eta, m, self.z, self.method)
   1138             else:
-> 1139                 self.X, self.R = _do_schur(self.P, self.eta, m, self.z, self.method)
   1140 
   1141     def minChi(self, m_min, m_max):

~/Projects/msmtools/msmtools/analysis/dense/gpcca.py in _do_schur(P, eta, m, z, method, tol_krylov)
    300         # TODO @Marius: I'd decrease the tolerance, getting error here for 3+ MS states
    301         print(X.conj().T.dot(np.diag(eta)).dot(X))
--> 302         raise ValueError("Schur vectors appear to not be D-orthogonal.")
    303     # Raise, if X doesn't fullfill the invariant subspace condition!
    304 

ValueError: Schur vectors appear to not be D-orthogonal.

The deviations from the unit matrix are quite large (~1e-1), so this is not a tolerance problem. Would you expect the current implementation to run on reducible matrices? I remember you had some handling for this case in PCCA.

[msmdev]

What are you using as eta here? The uniform distribution or something else? Differently asked: does X.T.dot(X) behave well?

[msmdev]

Another question: You say that you have "5 recurrent and a number of transient states". This is a little bit confusing to me, since a state is either transient or recurrent (to be perfectly precise: a state is either transient, positive recurrent or null recurrent). Since you say that your matrix is huge, I'm confused that you seem to only have a few states at all (5 recurrent plus a few transient ones). Shouldn't every state of your many (thousands?) states be either recurrent or transient? Also: How did you identify (technically) that your matrix is reducible and how did you characterize the states as recurrent or transient?

[msmdev]

Could you mb. provide me with the matrix?

[Marius1311]

What are you using as eta here? The uniform distribution or something else? Differently asked: does X.T.dot(X) behave well?

I'm using the uniform distribution.

[msmdev]

What are you using as eta here? The uniform distribution or something else? Differently asked: does X.T.dot(X) behave well?

I'm using the uniform distribution.

Hm, ok... thus it should work - there is no theoretical justification for your problem. Reducible matrices have a Schur decomposition too. You will have multiple eigenvalues equal 1, but this is no problem. From the transient states several eigenvalues at zero will result, thus the matrix won't have full rank, but this also shouldn't be a problem...

[Marius1311]

Another question: You say that you have "5 recurrent and a number of transient states". This is a little bit confusing to me, since a state is either transient or recurrent (to be perfectly precise: a state is either transient, positive recurrent or null recurrent). Since you say that your matrix is huge, I'm confused that you seem to only have a few states at all (5 recurrent plus a few transient ones). Shouldn't every state of your many (thousands?) states be either recurrent or transient? Also: How did you identify (technically) that your matrix is reducible and how did you characterize the states as recurrent or transient?

Sorry, I could have been more precise here. I meant that I have 5 recurrent classes, some of which only contain one state, i.e. are absorbing, see Fig 1 below. The remaining states are distributed across transient classes, I have 602 different transient classes, most of which only contain a single state. The largest transient class contains ~1.5k states. To compute recurrent and transient classes, I use this function in CellRank: https://github.com/theislab/cellrank/blob/040e758b3c1736d64c8ce9e29b2cc2ed198d0916/cellrank/tools/_utils.py#L694

I will send you the matrix via email.

Fig 1

[Marius1311]

Matrix send!

[msmdev]

Update from Matlab side:

The problem was not reproducible in Matlab: When sorting the top ten eigenvalues up and afterwards checking X.T*X for e.g. the 5 top Schurvectors (associated to 5 degenerated ones), the result was perfect orthogonality... subspace was in the range of 1e-15 too... Even clustering in e.g. 5 states worked perfectly - the vectors constituted a perfect simplex. Chi's looked fine too.

Conclusion:

Let's do the merge of request 3 (just change the tolerance to 2*eps here) and then try it again: The problem must be somewhere inside our code or in the scipy Schur routines... Something you could try: perform sorting of Q, R with Brandts once, afterwards sort the already sorted Q, R again with Brandts and check for changes... Have you tested the same matrix using Krylov? Any problems there?

[Marius1311]

Hi @msmdev , thanks for the update! Will look into the re-sorting with Brandts

I tried this both with krylov and brands, see above, same result.