ctmc_solve.h File Reference

Steady-state distribution of a continuous-time Markov chain. More...

#include <algorithm>
#include <cstddef>
#include <vector>
#include "line/api/mc/ctmc_bicgstab.h"
#include "line/api/mc/ctmc_gmres.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"

Include dependency graph for ctmc_solve.h:

Go to the source code of this file.

Namespaces
namespace	line
namespace	line::mc

Functions
template<class T>
Matrix< T >	line::mc::ctmc_makeinfgen (const Matrix< T > &Q)
	Set the diagonal so that every row sums to zero (ctmc_makeinfgen).
template<class T>
std::vector< T >	line::mc::ctmc_solve (const Matrix< T > &Qin)
	Steady-state distribution of a continuous-time Markov chain.

Detailed Description

Steady-state distribution of a continuous-time Markov chain.

Templated port of matlab/lib/kpctoolbox/mc/ctmc_solve.m (the numeric branch) and jar/src/main/java/jline/api/mc/Ctmc_solve.java. Solves pi Q = 0 with sum(pi) = 1 by replacing the last column of the generator with ones and solving Q' x = e_n, exactly as MATLAB's Qnnz' \ bnnz does.

Every step is a field operation, so the exact instantiation returns pi as a vector of rationals: for a generator with rational rates that is the true stationary distribution with no rounding whatsoever, which is what makes the exact backend worth its cost here (the double path of the JAR needs 50.9 s at n=800 where this LU needs tens of ms, see _kb/14-cpp-multiprecision).

Reducible generators are handled as MATLAB does: the weakly connected components are solved separately and renormalized. States that are ISOLATED (no transition in and none out) are trimmed first – an absorbing state is NOT isolated and is kept, being where the stationary mass ends up – and a generator that trims to nothing is an error rather than a plausible-looking uniform vector.

ABOVE GMRES_MIN_STATES THE LU IS ABANDONED FOR THE KRYLOV PATH, as in the other three codebases: restarted GMRES first, BiCGSTAB if it reports a nonzero flag, and the direct solve only if both do. The gate is compile-time as well as size-based – an exact instantiation always takes the LU, because the iteration stops on a residual tolerance and normalizes by a Euclidean norm, so there is nothing exact for it to converge to.

Definition in file ctmc_solve.h.