Convergence tweaks for reversible MLE and sampler

msmtools/estimation/api.py

@@ -832,12 +832,14 @@ def transition_matrix(C, reversible=False, mu=None, method='auto', **kwargs):
         stationary probabilities (:math:`x_i = \sum_k x_{ik}`). The relative stationary probability changes
         :math:`e_i = (x_i^{(1)} - x_i^{(2)})/(x_i^{(1)} + x_i^{(2)})` are used in order to track changes in small
         probabilities. The Euclidean norm of the change vector, :math:`|e_i|_2`, is compared to maxerr.
-    return_statdist : False : Boolean
+    return_statdist : bool, default=False
         Optional parameter with reversible = True.
         If set to true, the stationary distribution is also returned
-    return_conv : False : Boolean
+    return_conv : bool, default=False
         Optional parameter with reversible = True.
         If set to true, the likelihood history and the pi_change history is returned.
+    warn_not_converged : bool, default=True
+        Prints a warning if not converged.
 
     Returns
     -------

msmtools/estimation/dense/sampler_rev.pyx

@@ -96,9 +96,8 @@ class SamplerRev(object):
 
         """Set up initial state of the chain"""
         if P0 is None:
-            P0 = tmatrix(C, reversible=True)
-            # A = C + C.T
-            # V0 = A/A.sum()            
+            # only do a few iterations to get close to the MLE and suppress not converged warning
+            P0 = tmatrix(C, reversible=True, maxiter=100, warn_not_converged=False)
         pi0 = statdist(P0)
         V0 = pi0[:,np.newaxis] * P0            
 

msmtools/estimation/dense/transition_matrix.py

@@ -97,7 +97,7 @@ def __relative_error(x, y, norm=None):
 
 
 def estimate_transition_matrix_reversible(C, Xinit=None, maxiter=1000000, maxerr=1e-8,
-                                          return_statdist=False, return_conv=False):
+                                          return_statdist=False, return_conv=False, warn_not_converged=True):
     """
     iterative method for estimating a maximum likelihood reversible transition matrix
 
@@ -121,10 +121,12 @@ def estimate_transition_matrix_reversible(C, Xinit=None, maxiter=1000000, maxerr
         stationary probabilities (x_i = sum_k x_ik). The relative stationary probability changes
         e_i = (x_i^(1) - x_i^(2))/(x_i^(1) + x_i^(2)) are used in order to track changes in small
         probabilities. The Euclidean norm of the change vector, |e_i|_2, is compared to convtol.
-    return_statdist = False : Boolean
+    return_statdist : bool, default=False
         If set to true, the stationary distribution is also returned
-    return_conv = False : Boolean
+    return_conv : bool, default=False
         If set to true, the likelihood history and the pi_change history is returned.
+    warn_not_converged : bool, default=True
+        Prints a warning if not converged.
 
     Returns
     -------
@@ -193,7 +195,7 @@ def estimate_transition_matrix_reversible(C, Xinit=None, maxiter=1000000, maxerr
         i += 1
     # finalize and return
     T = X / xsum[:, np.newaxis]
-    if not converged:
+    if warn_not_converged and not converged:
         warnings.warn("Reversible transition matrix estimation didn't converge.", msmtools.util.exceptions.NotConvergedWarning)
     if (return_statdist and return_conv):
         return (T, xsum, lhist[0:i], diffs[0:i])
@@ -204,7 +206,8 @@ def estimate_transition_matrix_reversible(C, Xinit=None, maxiter=1000000, maxerr
     return T  # else just return T
 
 
-def transition_matrix_reversible_fixpi(Z, mu, maxerr=1e-10, maxiter=10000, return_iterations=False):
+def transition_matrix_reversible_fixpi(Z, mu, maxerr=1e-10, maxiter=10000, return_iterations=False,
+                                       warn_not_converged=True):
     r"""
     maximum likelihood transition matrix with fixed stationary distribution
 
@@ -223,6 +226,8 @@ def transition_matrix_reversible_fixpi(Z, mu, maxerr=1e-10, maxiter=10000, retur
         Will exit when reaching maxiter iterations without reaching convergence.
     return_iterations: bool (False)
         set true in order to return (T, it), where T is the transition matrix and it is the number of iterations needed
+    warn_not_converged : bool, default=True
+        Prints a warning if not converged.
 
     Returns
     -------
@@ -267,7 +272,7 @@ def transition_matrix_reversible_fixpi(Z, mu, maxerr=1e-10, maxiter=10000, retur
         # copy new to old l-vector
         l[:] = lnew[:]
         it += 1
-    if (not converged) and (it >= maxiter):
+    if warn_not_converged and (not converged) and (it >= maxiter):
         warnings.warn('NOT CONVERGED: 2-norm of Langrange multiplier vector is still ' +
                          str(err) + ' > ' + str(maxerr) + ' after ' + str(it) +
                          ' iterations. Increase maxiter or decrease maxerr',

msmtools/estimation/sparse/mle_trev.pyx

@@ -15,7 +15,7 @@ import msmtools.util.exceptions
 cdef extern from "_mle_trev.h":
   int _mle_trev_sparse(double * const T_data, const double * const CCt_data, const int * const i_indices, const int * const j_indices, const int len_CCt, const double * const sum_C, const int dim, const double maxerr, const int maxiter)
 
-def mle_trev(C, double maxerr = 1.0E-12, int maxiter = 1000000):
+def mle_trev(C, double maxerr = 1.0E-12, int maxiter = 1000000, warn_not_converged = True):
 
   assert maxerr > 0, 'maxerr must be positive'
   assert maxiter > 0, 'maxiter must be positive'
@@ -53,7 +53,7 @@ def mle_trev(C, double maxerr = 1.0E-12, int maxiter = 1000000):
     raise Exception('The update of the stationary distribution produced zero or NaN.')
   elif err == -3:
     raise Exception('Some row and corresponding column of C have zero counts.')
-  elif err == -5:
+  elif err == -5 and warn_not_converged:
     warnings.warn('Reversible transition matrix estimation didn\'t converge.', msmtools.util.exceptions.NotConvergedWarning)
 
   # T matrix has the same shape and positions of nonzero elements as CCt

msmtools/estimation/sparse/mle_trev_given_pi.pyx

@@ -27,7 +27,8 @@ def mle_trev_given_pi(
   mu,
   double maxerr = 1.0E-12,
   int maxiter = 1000000,
-  double eps = 0.0
+  double eps = 0.0,
+  warn_not_converged = True
   ):
 
   assert maxerr > 0, 'maxerr must be positive'
@@ -70,7 +71,7 @@ def mle_trev_given_pi(
     raise Exception('Some row and corresponding column of C have zero counts.')
   elif err == -4:
     raise Exception('Some element of pi is zero.')
-  elif err == -5:
+  elif err == -5 and warn_not_converged:
     warnings.warn('Reversible transition matrix estimation with fixed stationary distribution didn\'t converge.', msmtools.util.exceptions.NotConvergedWarning)
 
   # unnormalized T matrix has the same shape and positions of nonzero elements as the C matrix