api.mc

ctmc_saddlepoint(D0, D1, t, k, method, pi0)

[P, LOGP, THETA, INFO] = CTMC_SADDLEPOINT(D0, D1, T, K, METHOD, PI0)

Saddlepoint approximation of Pr{N(t)=k}, the probability that the counting process of the Markovian arrival process (D0,D1) records exactly k events in (0,t], starting from the phase distribution PI0. Throughout, K is the number of phases and k the event count.

The counting generating function of a MAP is the matrix exponential

sum_{k>=0} P(k,t) z^k = exp(t*(D0 + z*D1)),

where P(k,t) is the (i,j) probability of k events and phase j at time t. Its cumulant generating function is therefore

eta(theta) = spectral abscissa of A(theta) = D0 + exp(theta)*D1,

the Perron root of an irreducible Metzler matrix, which is real, simple and strictly convex in theta with eta(0)=0 and eta’(0)=lambda, the arrival rate. Inverting the generating function by the method of steepest descent gives Daniels (1954),

Pr{N(t)=k} ~ g(theta*) * exp(t*eta(theta*) - k*theta*) / sqrt(2*pi*t*eta’’(theta*)),

where the saddle theta* solves eta’(theta*) = k/t and g is the amplitude of the Perron projection, g(theta) = (pi0*v(theta))*(u(theta)*1) with u,v the left and right Perron vectors normalised by u*v = 1.

THE EXPANSION PARAMETER IS K2 = t*eta’’(theta*), THE VARIANCE OF THE COUNT, not its mean and not t. Measured across the process family, the two error laws are

err(‘daniels’) = 0.083 / K2, err(‘daniels2’) = 0.017 / K2^2,

with the constants flat to two digits over Erlang orders 1..8 and horizons 10..160. For a renewal Erlang(r) the count variance rate is lambda/r, so K2 = lambda*t/r and an Erlang-4 at t=50 is as accurate as a Poisson at t=12.5: low variability SHRINKS the parameter, it does not break the method. Below K2 = 5 the expansion is outside its regime and the call warns; INFO.K2 carries the value per point.

The approximation therefore IMPROVES with the horizon, which is where the block-by-block expansions of exp(t*X) lose their digits to cancellation, and it keeps relative accuracy on probabilities far below the double-precision floor when read through LOGP. Its cost is one eigenproblem of order K per Newton step and is independent of both t and k.

This is an asymptotic method, not a quadrature: use it for rare-event and large-deviation coefficients, where k/t is away from lambda or where the probability underflows. For the bulk of the transient distribution, i.e. every block k=0..N-1 at once at moderate t, uniformization (CTMC_UNIFORMIZATION, CTMC_FOXGLYNN) is both exact and faster.

@param D0 Generator of the phase process with the counted transitions: removed (K x K). May instead be a MAP cell {D0,D1}, in which case the remaining arguments shift left by one.
@param D1 Rates of the counted transitions (K x K, nonnegative). D0+D1 must: be an irreducible generator.

@param t Time horizon; scalar, or an array broadcast against the count @param k Event count, a nonnegative integer; scalar, or an array broadcast

against the horizon

@param method ‘daniels2’ (default) second-order saddlepoint, the Daniels: bracket plus the amplitude’s own curvature along the contour, error O(1/K2^2); ‘daniels’ is the first-order form with the Perron amplitude, error O(1/K2); ‘plain’ is the bare first-order form with the amplitude set to 1
@param pi0 Initial phase distribution (1 x K); the stationary distribution: of D0+D1 if empty or omitted

@return p Approximation of Pr{N(t)=k}, of the broadcast size of t and k @return logp Its natural logarithm, evaluated without forming P, so that it

stays accurate below the smallest positive double

@return theta The saddle theta*, -Inf where k=0, the one point at which it: runs off to minus infinity
@return info Struct with the per-point fields eta, deta, d2eta, d3eta,: d4eta, ampl (the Perron amplitude g), corr (the bracket multiplying the leading term), k2 (t*eta’’(theta*), the expansion parameter – read this to know how far the answer can be trusted), iter (Newton steps), exact (true where the value was computed exactly rather than approximated), and the scalar field lambda (the stationary event rate eta’(0))

Examples

MAP = {[-1 0.2 0.1; 0.05 -2 0.3; 0.1 0.1 -0.9], …: [0.4 0.2 0.1; 0.5 0.15 1.0; 0.3 0.2 0.2]};

p = ctmc_saddlepoint(MAP, 200, 66) % one coefficient [p, logp] = ctmc_saddlepoint(MAP, 1e6, 330298); % logp = -1.95e5 p = ctmc_saddlepoint(MAP, 100, 120:140) % a stretch of counts [p, ~, ~, info] = ctmc_saddlepoint(MAP, 40, 34); % info.k2 = 36.06

Measured relative error on that MAP, against Pr{N(t)=k} read off the block chain by dense expm (independently cross-checked against uniformization summed from n=0):

t k K2 Pr{N(t)=k} daniels daniels2 plain

20 26 27.2 1.109e-02 3.1e-03 6.0e-06 5.0e-03

100 132 137.8 9.749e-07 6.1e-04 2.4e-07 2.7e-03 200 66 71.4 5.354e-19 1.1e-03 6.2e-07 1.1e-02 500 165 178.6 1.316e-44 4.6e-04 1.0e-07 1.1e-02 500 660 688.8 2.984e-25 1.2e-04 9.6e-09 2.2e-03

The three columns are the three regimes of the expansion: ‘plain’ carries an amplitude error that does not vanish, ‘daniels’ is O(1/K2), ‘daniels2’ is O(1/K2^2). Over the point processes and horizons of LINE’s own transient examples (MAP, MMPP2, Erlang, APH and Exp sources at t = 5..1000; 25 cases, 222 sampled counts at k/lambda*t in [0.15, 2.5]), ‘daniels2’ was never worse than ‘daniels’ at any point, and better by one to three orders almost everywhere. Where both fail is small K2, and there they fail together: the second-order correction goes nonpositive, the call falls back to first order and warns.

ATTRIBUTION. The first-order form is Daniels (1954). The AMPLITUDE g and the whole ‘daniels2’ bracket are NOT a rederivation: they are Jensen, “Saddlepoint Expansions for Sums of Markov Dependent Variables on a Continuous State Space”, Probab. Th. Rel. Fields 89, 1991, Eq. (4.4) with the coefficients on p.191. His gamma_0(s) = (sum_i c_i)(sum_i r_i P(Y_0=i)) is exactly g, under his own normalisation sum_i r_i c_i = 1; expanding his alpha_0 + (1/n){-alpha_3/2 + alpha_4/8 - 5*alpha_5/24} reproduces g*(1 + lam4/8 - 5*lam3^2/24) - g’’/(2*K2) + g’*K3/(2*K2^2) term for term, and his Theorem 4.1 gives the O(n^-2) error this file measures as 0.017/K2^2. Jensen works with discrete-n sums over a Markov chain; the continuous-time MAP counting process here is that result transcribed, n -> t and the kernel eigenvalue -> the Perron root of D0+exp(theta)*D1. He also states the reason the amplitude cannot be dropped (p.183): the large-deviation limit discards any factor of order one, but a local approximation must carry the projection onto the eigenspace, not just the maximal eigenvalue.

mexify_mc

@brief MATLAB Coder script to generate MEX functions for mc_ module.

This script generates MEX (MATLAB Executable) versions of Markov chain functions for improved performance.

Skipped functions (Coder-incompatible):: ctmc_courtois, ctmc_kms - cell array args (MS), dtmc_solve dependency ctmc_takahashi - cell array args, dtmc_solve, ctmc_randomization ctmc_multi - cell array args, ctmc_solve_reducible ctmc_simulate - exprnd() not supported ctmc_uniformization - speye(), sparse() not supported ctmc_transient - function handles (@ctmc_transientode) ctmc_ssg, ctmc_ssg_reachability - State class OOP methods dtmc_solve_reducible - stronglyconncomp() not supported ctmc_solve_reducible - depends on dtmc_solve_reducible ctmc_solve_reducible_blkdecomp - stronglyconncomp() not supported

See also CODER, CODER.CONFIG, CODER.TYPEOF, CODEGEN.

ctmc_takahashi(Q, MS, numSteps): CTMC_TAKAHASHI - Takahashi’s aggregation-disaggregation method [p,p_1,pcourt,Qperm,eps,epsMAX] = CTMC_TAKAHASHI(Q,MS,numSteps) – Input Q : infinitesimal generator matrix MS : cell array where MS{i} is the set of rows of Q in macrostate i numSteps: number of iterative steps – Output p : estimated steady-state probability vector p_1 : pcourt : steady-state probability vector estimated by ctmc_courtois Qperm : permuted Q matrix w.r.t MS as returned by ctmc_courtois eps : NCD index as returned by ctmc_courtois epsMAX : max acceptable value for eps (otherwise Q is not NCD) – Remarks * The initial approximate solutions is obtained by calling CTMC_COURTOIS(Q,MS) * No convergence stop criterion is currently implented

ctmc_state_space_logsize(sn, options)

LOGNSTATES = CTMC_STATE_SPACE_LOGSIZE(SN, OPTIONS)

Worst-case log-size of the CTMC state space induced by SN. The estimate is the product of four factors, summed in log space:

job placements: stars-and-bars C(n_k+M-1,M-1) per class, over the stations that do NOT keep an ordered buffer, with open classes truncated at the cutoff;

buffer orderings: a station outside the share family keeps the CLASS SEQUENCE of the jobs it holds, so with K>1 classes each occupancy vector is as many states as its sequences. Bounded by the number of sequences of length at most sum(n_k) over K symbols, once per such station;

service phases: sn.phasessz raised to the number of jobs that can be in service concurrently at each station;

routing state: one pointer over the outgoing links per (node,class) using RROBIN or WRROBIN.

This is the quantity fed to CTMC_MEMORY_GATE. It is exposed separately so that SolverAUTO can screen CTMC out of a candidate ranking without building the chain; see SolverCTMC.isStateSpaceTractable.

ctmc_multi(Q, MS, MSS): CTMC_MULTI - Multigrid aggregation-disaggregation method CTMC_Multigrid is the basic one step implementation of standard method of Multigrid. The complete implementation (e.g multi-level disaggregation) requires repetitive coarsening using the base method. For illustrative purpose, we only use two levels.

ctmc_kms(Q, MS, numSteps): CTMC_KMS - Koury-McAllister-Stewart aggregation-disaggregation method [p,p_1,Qperm,eps,epsMAX,pcourt] = CTMC_KMS(Q,MS,numSteps) – Input Q : infinitesimal generator matrix MS : cell array where MS{i} is the set of rows of Q in macrostate i numSteps: number of iterative steps – Output p : estimated steady-state probability vector p_1 Qperm : permuted Q matrix w.r.t MS as returned by ctmc_courtois eps : NCD index as returned by ctmc_courtois epsMAX : max acceptable value for eps (otherwise Q is not NCD) pcourt : steady-state probability vector estimated by ctmc_courtois – Remarks * The initial approximate solutions is obtained by calling CTMC_COURTOIS(Q,MS) * No convergence stop criterion is currently implented

ctmc_foxglynn(pi0, Q, t, tol, maxiter)

[PI,LEFT,RIGHT,W]=CTMC_FOXGLYNN(PI0,Q,T,TOL,MAXITER)

Transient distribution of the CTMC by uniformization, with Poisson weights and truncation points obtained by the Fox-Glynn algorithm, with Jansen’s correction to the right tail estimate. Neither exp(-q*t) nor (q*t)^k/k! is ever formed, so the method is free of overflow and underflow and needs no horizon splitting. MAXITER caps the right truncation point; pass a nonpositive value (or omit it) to leave it uncapped.

ctmc_courtois(Q, MS, q): CTMC_COURTOIS - Courtois decomposition [p,Qperm,Qdec,eps,epsMAX,P,B,C,q] = CTMC_COURTOIS(Q,MS) – Input Q : infinitesimal generator matrix MS : cell array where MS{i} is the set of rows of Q in macrostate i q : (optional) randomization coefficient – Output p : approximate steady-state probability vector Qperm : Q reordered according to macrostates Qdec : infinitesimal generator for the macrostates P : probability matrix obtained from Qperm with randomization B : part of P not modelled by decomposition eps : nearly-complete decomposability (NCD) index epsMAX : max acceptable value for eps (otherwise Q is not NCD) q : randomization coefficient

smp_passage_time(P, hlst, pi0, target, tset, options)

[F, f, OUT] = SMP_PASSAGE_TIME(P, HLST, PI0, TARGET, TSET, OPTIONS)

Cumulative distribution and density of the first passage time into the target state set for a semi-Markov chain, by inverting SMP_PASSAGE_LST through api/lti.

There is no matrix-exponential route here: a semi-Markov chain has no generator to exponentiate, which is exactly the case uniformization does not reach and the transform does. This is the point the paper makes for preferring transform inversion over uniformization.

OPTIONS.lti_method defaults to ‘euler’ RATHER THAN ‘weeks’. Semi-Markov passage densities are the case Sec. 4.2 singles out as slow-converging for a Laguerre series: a kernel with a deterministic or discontinuous holding time gives a density whose derivatives jump, and LAPLACE_WEEKS_SCALING then refuses by name rather than returning noise.

smp_passage_lst(P, hlst, pi0, target, s)

L = SMP_PASSAGE_LST(P, HLST, PI0, TARGET, S)

Laplace-Stieltjes transform of the first passage time into the target state set for a semi-Markov chain, at the (possibly complex) points S.

Reference: P. G. Harrison and W. J. Knottenbelt, “Passage Time Distributions in Large Markov Chains”, 2002, Eqs. 4-5:

L_i(s) = sum_{k not in B} r*_ik(s) L_k(s) + sum_{k in B} r*_ik(s)

so (I - R*_AA(s)) L_A(s) = R*_AB(s) 1, one linear system per value of s.

HLST is either

(nstates x 1) cell of handles h*_i(s) the sojourn in i depends only on i,: so r*_ik(s) = P(i,k) h*_i(s) and the complex numbers stay on the DIAGONAL of the system (Eq. 5)
(nstates x nstates) cell of r*_ik(s) the full Markov-renewal kernel; the: coefficients are then complex functions of s throughout, which is the harder case the paper flags

Distribution objects supply their own transform: MARKOVIAN.evalLST gives the closed form pie (sI-D0)^{-1} (-D0) e for the phase-type family, so hlst{i} = @(s) dist.evalLST(s) is the intended way to build these.

ctmc_passage_moments(Q, pi0, target, nmax)

[MALL, M] = CTMC_PASSAGE_MOMENTS(Q, PI0, TARGET, NMAX)

Moments of order 1..NMAX of the first passage time into the target state set. MALL is (nstates x NMAX): row i is the moment vector for a passage STARTED IN STATE i, zero on target states and Inf where the target cannot be reached. M is the (1 x NMAX) moment vector for the initial law PI0.

This is Eq. 3 of P. G. Harrison and W. J. Knottenbelt, “Passage Time Distributions in Large Markov Chains”, 2002,

-q_ii M_i(n) = sum_{k not in target} q_ik M_k(n) + n M_i(n-1)

i.e. (-S) M(n) = n M(n-1) with M(0) = 1, solved once per order: NMAX sparse solves and no transform inversion at all. The equivalent closed form is n! alpha (-S)^{-n} 1, which is NOT how it is evaluated here – forming the inverse of the sub-generator destroys the sparsity the recursion preserves.

ctmc_passage_lst(Q, pi0, target, s)

L = CTMC_PASSAGE_LST(Q, PI0, TARGET, S)

Laplace-Stieltjes transform of the first passage time from PI0 into the target state set, evaluated at the (possibly complex) points S.

L(s) = alpha (sI - S)^{-1} s0 + atom

which is Eqs. 1-2 of P. G. Harrison and W. J. Knottenbelt, “Passage Time Distributions in Large Markov Chains”, 2002: one linear system per value of s, of the size of the non-target block.

ONE SPARSE SOLVE PER s, NOT PER (s,t) PAIR. The saving over a dense matrix exponential is that the solves are sparse, so this route reaches chains a dense expm cannot hold. It is NOT a saving in the number of time points: every Abate-Whitt inverter places its nodes at s = beta/t, so a grid of T points costs T*|beta| solves. On a small chain ctmc_passage_time’s default ‘expm’ route is faster.

ctmc_hitting_time(Q, targetStates)

H = CTMC_HITTING_TIME(Q, TARGETSTATES)

Mean time to reach any state in TARGETSTATES from each state of a CTMC with generator Q. Target states have zero hitting time; a state that cannot reach the set has an infinite one.

Continuous-time twin of DTMC_HITTING_TIME, and the first-moment special case of CTMC_PASSAGE_MOMENTS: (-S) h = 1 on the non-target block, where DTMC_HITTING_TIME solves (I - P_NT) h = 1. TARGETSTATES is 1-based here and 0-based in the Java, Python and C++ twins.

dtmc_transient(P, pi0, steps)

PI_T = DTMC_TRANSIENT(P, PI0, STEPS)

Transient distribution of a DTMC with transition matrix P, i.e. the rows PI(k) = PI0*P^k for k = 0,…,STEPS. Twin of the Python api.mc.dtmc_transient.

dtmc_solve_reducible(P, pin, options)

[PI,PIS,PI0,SCC,ISREC,PL,PIL] = DTMC_SOLVE_REDUCIBLE(P, PIN, OPTIONS)

Estimate limiting distribution for a DTMC P that may have reducible components

Input: P: dtmc transition matrix pin: initial vector, to be set to [] if not available options: struct where options.tol sets the tolerance

Output: pi: - For an ergodic DTMC, this is the unique limiting distribution.

For a reducible DTMC:

if there is a single transient SCC then this is assumed to be the starting state with uniform probability across its states.

if there are multiple transient SCCs then pi is the average of the limiting distributions in pis

pis: limiting distribution given initialization pi0 in a single SCC pi0: start vector of the lumped DTMC for each row of pis. For ergodic

DTMCs, this is the empty vector.

scc: mapping of state of P to SCCs isrec: element i is true if SCC i is recurrent or false otherwise Pl: lumped DTMC where each SCC is replaced by a single state pil: limiting distribution of the lumped DTMC

WARNING: the script does not consider explicitly periodic SCCs

dtmc_hitting_time(P, targetStates)

H = DTMC_HITTING_TIME(P, TARGETSTATES)

Mean number of steps to reach any state in TARGETSTATES from each state of a DTMC with transition matrix P. Target states have zero hitting time; the others solve (I - P_NT) h_NT = 1 over the non-target block. Twin of the Python api.mc.dtmc_hitting_time; TARGETSTATES is 1-based here and 0-based there, as elsewhere between the two codebases.

ctmc_isfeasible(Q, tol)

BOOL = CTMC_ISFEASIBLE(Q, TOL)

True when Q is a valid infinitesimal generator: square, non-negative off-diagonal entries, non-positive diagonal, and zero row sums, each up to TOL (default 1e-10). Twin of the Python api.mc.ctmc_isfeasible; note that dtmc_isfeasible instead returns a precision level rather than a flag.

ctmc_transient_sens(Q, dQ, pi0, t0, t1)

[DPI, PI, T] = CTMC_TRANSIENT_SENS(Q, DQ, PI0, T0, T1)

Sensitivity of the transient distribution of a CTMC to a scalar parameter theta, given the generator Q and its derivative DQ = dQ/dtheta.

Differentiating the forward equations d pi(t)/dt = pi(t) Q with respect to theta, and assuming the initial vector does not depend on theta, gives

d/dt (dpi(t)/dtheta) = (dpi(t)/dtheta) Q + pi(t) (dQ/dtheta), dpi(0)/dtheta = 0,

i.e. Trivedi and Bobbio (2017), Eq. (9.82). The state and its sensitivity are integrated as one augmented system of size 2n, since the sensitivity equation is driven by pi(t) and the two cannot be advanced separately.

@param Q Generator matrix (n x n) @param dQ Derivative of the generator with respect to theta (n x n) @param pi0 Initial distribution (1 x n); uniform if empty @param t0 Initial time; 0 if omitted @param t1 Final time @return dpi Sensitivity of the distribution at each time point (length(T) x n) @return pi Distribution at each time point (length(T) x n) @return t Column vector of time points

ctmc_sens(Q, dQ, pi)

DPI = CTMC_SENS(Q, DQ, PI)

Sensitivity of the steady-state distribution of a CTMC to a scalar parameter theta, given the generator Q, its derivative DQ = dQ/dtheta, and the steady-state vector PI.

Differentiating the balance equations pi*Q = 0 and pi*e = 1 with respect to theta gives the linear system

(dpi/dtheta) * Q = -pi * (dQ/dtheta), sum_i dpi_i/dtheta = 0,

i.e. Trivedi and Bobbio (2017), Eq. (9.81). The system has the same coefficient matrix as the steady-state solve itself, so obtaining a sensitivity costs one extra solve against a matrix that is already assembled. The normalization replaces one column of the singular Q, exactly as in the steady-state solve.

@param Q Generator matrix (n x n) @param dQ Derivative of the generator with respect to theta (n x n) @param pi Steady-state distribution (1 x n); computed if omitted @return dpi Derivative of the steady-state distribution (1 x n)

lineGetAvailableMemory()

LINEGETAVAILABLEMEMORY Portable available-physical-memory probe (bytes).

Routes through the JVM OperatingSystemMXBean that MATLAB ships on every platform, so the same code path returns a valid figure on Windows, macOS and Linux without parsing /proc. Falls back to a conservative constant if the bean or method is unavailable (e.g. MATLAB started with -nojvm).

Mirrors jline.solvers.ctmc.MemoryGuard.getAvailableMemoryBytes (JAR) and line_solver … memory_guard.get_available_memory_bytes (Python native).

ctmc_uniformization(pi0, Q, t, tol, maxiter)

[PI,KMAX]=CTMC_UNIFORMIZATION(PI0,Q,T,TOL,MAXITER)

MAXITER caps the Poisson series truncation depth; pass a nonpositive value (or omit it) to size it adaptively as max(100, q*t+10*sqrt(q*t)+20).

ctmc_timeaverage(pi0, Q, t, tol, maxiter)

[PITIMEAVG,PIEXIT,KMAX]=CTMC_TIMEAVERAGE(PI0,Q,T,TOL,MAXITER)

Time-averaged transient distribution of a CTMC with generator Q over [0,T], starting from the (arbitrary) initial distribution PI0, via Jensen’s uniformization. Companion of CTMC_UNIFORMIZATION, which returns only the endpoint PI0*exp(Q*T); this function additionally returns the time average

PITIMEAVG = PI0 * (1/T) * int_0^T exp(Q*tau) d(tau)

as well as the endpoint PIEXIT = PI0*exp(Q*T) (computed from the same series). Both are obtained without forming any dense matrix exponential.

Uniformization: with q = 1.1*max|diag(Q)| and P = I + Q/q (row-stochastic),: PI0*exp(Q*T) = sum_j w_j(qT) * (PI0*P^j) PI0*int_0^T exp(Q*tau) = (1/q) * sum_j (1 - W_j(qT)) * (PI0*P^j)

where w_j and W_j are the Poisson(qT) PMF and CDF. The time average divides the integral by T (equivalently the integral sum by q*T).

ctmc_ssg(sn, options)

ctmc_simulate(Q, pi0, n)

ctmc_memory_gate(logNstates, force, verbose, safetyFraction)

CTMC_MEMORY_GATE Hardware-aware, profiling-calibrated CTMC memory pre-gate.

[ok,msg] = ctmc_memory_gate(logNstates, force, verbose, safetyFraction)

Decides whether a CTMC steady-state solve of a state space of worst-case size exp(logNstates) is safe on the current host. The budget is a fraction of the memory actually available (see lineGetAvailableMemory); the per-state cost is calibrated by profiling the local sparse LU once and caching the fitted power law per machine in tempdir.

ok is false only when the predicted peak footprint exceeds the budget and force is not set. This replaces the historical hard-coded threshold.

ctmc_ssg_reachability(sn, options)

ctmc_solve_reducible(Q, pi0, options)

ctmc_rand(n): Q=CTMC_RAND(N)

ctmc_makeinfgen(Q): Q=CTMC_MAKEINFGEN(Q)

ctmc_stochcomp(Q, I): [S,Q11,Q12,Q21,Q22,T] = CTMC_STOCHCOMP(Q,I) Copyright (c) 2012-2026, Imperial College London All rights reserved.

ctmc_solve_reducible_blkdecomp(Q, pin, options)

[PI,PIS,PI0,SCC,ISREC] = CTMC_SOLVE_REDUCIBLE_BLKDECOMP(Q, PIN, OPTIONS)

Compute limiting distribution for a CTMC with reducible generator Q using direct block decomposition on the generator matrix.

Algorithm:

Decompose states into transient and recurrent classes via SCC
For transient states: solve sojourn * Q_tt = -p0_t for expected sojourn
Compute hitting probabilities: hit = sojourn * Q_ta + p0_r
For each recurrent class: solve pi_c * Q_cc = 0, scale by hitting prob

Input: Q: infinitesimal generator matrix pin: initial distribution vector, set to [] if not available options: struct where options.tol sets the tolerance

Output: pi: limiting distribution (1 x N)

For an ergodic CTMC, this is the unique limiting distribution.

For a reducible CTMC:

if PIN is supplied, this is the exact limiting distribution for that initial vector: each recurrent class (BSCC) carries its absorption probability from PIN and every transient state carries zero.

otherwise, if there is a single transient SCC then this is the limiting distribution when starting uniformly within it, and failing that pi is the weighted average of pis rows.

pis: limiting distribution given initialization in a single SCC (numSCC x N) pi0: starting distribution for each row of pis (numSCC x N), so PIN is the

vector the caller supplies and PI0 the ones this routine builds

scc: mapping of each state of Q to its SCC index isrec: element i is true if SCC i is recurrent

ctmc_gmres_multi(A, B, tol, restart, maxit)

[X,FLAG]=CTMC_GMRES_MULTI(A,B,TOL,RESTART,MAXIT)

Solve A*X=B for every column of B by restarted GMRES, reusing one ILUT factorization across all of them and starting each column from the previous solution. This is the shape of the stochastic complement, whose right-hand side is a whole block of the generator: refactorizing per column would cost more than the direct solve it replaces.

Preparation follows CTMC_GMRES: rows are equilibrated to unit max norm and the states reordered by reverse Cuthill-McKee, without which the unpivoted elimination overflows on a chain of a few thousand states.

FLAG is 0 only if every column converged. On any other value X must be discarded and the caller must fall back to the direct solve; returning a partial block would leave the fallback ambiguous.

ctmc_gmres(A, b, tol, restart, maxit, x0)

[X,FLAG,RELRES,ITER]=CTMC_GMRES(A,B,TOL,RESTART,MAXIT,X0)

Solve the sparse nonsymmetric linear system A*x=b by restarted GMRES with an ILUT right preconditioner, falling back to a Jacobi preconditioner when the incomplete factorization breaks down. This is the iterative counterpart of the direct sparse solve used by CTMC_SOLVE, intended for generators whose LU fill-in exceeds available memory.

Two preparation steps are not optional on a generator. Rows are equilibrated to unit max norm, so the O(1) normalization row does not mix with rows carrying rates of a different magnitude. The states are then reordered by reverse Cuthill-McKee: in the natural ordering of a birth-death chain the unpivoted elimination has growth factor (mu/lambda)^n, which overflows by a few thousand states, and a fill-reducing ordering rather than pivoting is what removes it.

A is the already-assembled coefficient matrix; no CTMC-specific processing is performed here, so the same kernel serves the stochastic complementation and aggregation kernels.

Defaults: TOL=1e-12, RESTART=min(n,50), MAXIT=ceil(n/RESTART), X0=ones(n,1)/n. TOL is a linear-solve residual and is therefore much tighter than the fixed-point tolerance options.iter_tol.

FLAG follows the MATLAB GMRES convention: 0 converged, 1 iteration limit reached, 2 preconditioner ill-conditioned, 3 stagnation. Callers must check it and fall back to the direct solve when it is nonzero.

ITER counts matrix-vector products with A, i.e. total inner iterations across all restart cycles, so that the four codebases report one comparable scalar.

The recurrence is in the private GMRES_ITERATE rather than the built-in GMRES, which applies the preconditioner on the LEFT; see that function for what that costs on a generator.

ctmc_fau(pi0, Q, t, epsilon, delta, maxsteps)

[PIT,INFO]=CTMC_FAU(PI0,Q,T,EPSILON,DELTA,MAXSTEPS)

Transient distribution of a CTMC at time T by fast adaptive uniformization.

ADAPTIVE UNIFORMIZATION (van Moorsel and Sanders, 1994). Ordinary uniformization fixes one rate q >= max_i |q_ii| over the WHOLE state space and mixes the powers of P = I + Q/q against a Poisson(q*t) law. Its cost is therefore set by the fastest state anywhere, including states that carry no probability at time t. Adaptive uniformization instead picks a rate per step from the states the iterate actually occupies,

Lambda_n >= max{ |q_ii| : i in supp(u^(n)) }, u^(n+1) = u^(n)(I + Q/Lambda_n),

which keeps every entry of u^(n+1) nonnegative. The subordinating process is then no longer Poisson but the pure birth process N(t) with rates Lambda_0, Lambda_1, …, and

pi(t) = sum_{n>=0} P{N(t)=n} u^(n).

FAST ADAPTIVE UNIFORMIZATION (Mateescu, Wolf, Didier and Henzinger, 2010) adds the second half: an entry of u^(n) below DELTA is dropped rather than propagated, so the support tracks the states of non-negligible occupancy instead of the reachable set. Dropping only removes nonnegative contributions, so PIT is a componentwise LOWER BOUND on the exact distribution. Nothing is renormalized anywhere, so the error is not estimated but MEASURED: the three approximations (the birth index truncated at K, the Poisson window of the weight computation, the DELTA threshold) each remove mass and none puts any back, whence

0 <= pi(t) - PIT componentwise, and |pi(t) - PIT|_1 = sum(PI0) - sum(PIT) = INFO.errorBound.

This is what replaces a blind population cutoff on an open model by a numerical one whose error is a returned quantity. It is a transient method: it produces no stationary distribution.

THE BIRTH WEIGHTS ARE COMPUTED EXACTLY, not quadratured. The rates Lambda_0..Lambda_K generate a bidiagonal generator B on the birth index 0..K plus one absorbing overflow index, and P{N(t)=n} is the transient distribution of that scalar chain, obtained by uniformizing it at Lstar = max_n Lambda_n and applying the shipped Fox-Glynn weights. Every entry of the uniformized bidiagonal kernel lies in [0,1], so there is no cancellation, and the mass reaching the overflow index IS the truncation error INFO.weightTail. The alternative found in most implementations, integrating the convolution b_n(s) = Lambda_{n-1} int_0^s b_{n-1}(v) exp(-Lambda_n (s-v)) dv on a time grid, carries a quadrature error that the reported bound would then have to absorb.

THE SWEEP RUNS TWICE, and this is the one real cost of the exact weights. b_n(t) needs the rates up to n, which are not known before the sweep ends, while u^(n) is needed after them; storing every iterate would cost K * |support| doubles. The sweep is therefore replayed, and since it is deterministic the second pass reproduces the first rate for rate and drop for drop. Two sparse sweeps still beat ordinary uniformization whenever Lstar is well below max_i |q_ii| or the occupied support is well below n.

STOPPING IS CERTIFIED, not heuristic. S_{K+1} = sum_{m<=K} Exp(Lambda_m) dominates an Erlang(K+1, Lstar) stochastically, so P{N(t) > K} = P{S_{K+1} <= t} <= P{Poisson(Lstar*t) >= K+1}, bounded above by the Chernoff exponent of that Poisson tail. The sweep stops at the first K meeting EPSILON, which also shows that this method never takes more steps than uniformization at the largest rate it visited.

– Input PI0 : initial distribution, 1xn Q : infinitesimal generator, nxn (sparse is used as given) T : time horizon, T >= 0 EPSILON : birth-process truncation tolerance (default 1e-6) DELTA : occupancy threshold below which a state is dropped (default 1e-12) MAXSTEPS : cap on birth steps; nonpositive for the default cap of 1e6

– Output PIT : 1xn defective distribution at time T, a lower bound on pi(T) INFO : struct with fields

steps number of birth steps K+1 actually taken lambdaMin smallest adaptive rate used lambdaMax largest adaptive rate used, the Lstar above uniformRate max_i |q_ii|, the rate ordinary uniformization would use weightTail mass reaching the overflow index, i.e. P{N(T) > K} weightWindow Poisson mass outside the Fox-Glynn window of the weights droppedMass probability removed by the DELTA threshold errorBound sum(PI0)-sum(PIT), which IS |pi(T)-PIT|_1 supportMax largest occupied support over the sweep supportFinal support at the last step truncated true if MAXSTEPS stopped the sweep absorbed true if the support emptied or became absorbing

ctmc_bicgstab_multi(A, B, tol, maxit)

[X,FLAG]=CTMC_BICGSTAB_MULTI(A,B,TOL,MAXIT)

Solve A*X=B for every column of B by stabilized biconjugate gradients, reusing one ILUT factorization across all of them and starting each column from the previous solution. This is the shape of the stochastic complement, whose right-hand side is a whole block of the generator: refactorizing per column would cost more than the direct solve it replaces.

Preparation follows CTMC_BICGSTAB: rows are equilibrated to unit max norm and the states reordered by reverse Cuthill-McKee, without which the unpivoted elimination overflows on a chain of a few thousand states. The iteration is the private BICGSTAB_ITERATE, the same one CTMC_BICGSTAB runs.

FLAG is 0 only if every column converged. On any other value X must be discarded and the caller must fall back to another solve; returning a partial block would leave the fallback ambiguous.

ctmc_bicgstab(A, b, tol, maxit, x0)

[X,FLAG,RELRES,ITER]=CTMC_BICGSTAB(A,B,TOL,MAXIT,X0)

Solve the sparse nonsymmetric linear system A*x=b by the stabilized biconjugate gradient method of van der Vorst (1992), with the same equilibration and reordering as CTMC_GMRES and a non-pivoting incomplete factorization (see the ILU note below, which is where the two kernels part company). It is the short-recurrence counterpart of that kernel: work and storage per iteration are constant rather than growing with the Krylov dimension, so it does not restart and does not lose the optimality that restarting costs GMRES. Where GMRES(m) stagnates because the useful subspace is wider than m, this converges; where it does not, GMRES(m) is the more robust of the two, hence the order in which CTMC_SOLVE tries them.

Preparation is not optional on a generator and is identical to CTMC_GMRES: rows are equilibrated to unit max norm, so the O(1) normalization row does not mix with rows carrying rates of a different magnitude, and the states are then reordered by reverse Cuthill-McKee, since in the natural ordering of a birth-death chain the unpivoted elimination has growth factor (mu/lambda)^n.

A is the already-assembled coefficient matrix; no CTMC-specific processing is performed here, so the same kernel serves the stochastic complementation and aggregation kernels.

Defaults: TOL=1e-12, MAXIT=min(n,200), X0=ones(n,1)/n. TOL is a linear-solve residual and is therefore much tighter than the fixed-point tolerance options.iter_tol.

FLAG follows the MATLAB BICGSTAB convention: 0 converged, 1 iteration limit reached, 2 preconditioner ill-conditioned, 3 stagnation, 4 a scalar quantity became too small or too large to continue. Callers must check it and fall back to another solve when it is nonzero.

ITER counts matrix-vector products with A: two per complete iteration, and one when the iteration converges at its half step. Reporting the product count rather than the iteration count is what makes it comparable with the ITER of CTMC_GMRES and across the four codebases, whose iteration bookkeeping differs.

The recurrence is in the private BICGSTAB_ITERATE rather than the built-in BICGSTAB, which applies the preconditioner on the LEFT and breaks down on the systems a generator produces; see that function for the measurement.

smp_passage_moments(P, hmom, pi0, target, nmax)

[MALL, M] = SMP_PASSAGE_MOMENTS(P, HMOM, PI0, TARGET, NMAX)

Moments of order 1..NMAX of the first passage time into the target state set for a semi-Markov chain with embedded transition matrix P and holding-time moments HMOM. MALL is (nstates x NMAX), row i for a passage started in state i, zero on the target. M is the PI0-weighted moment vector.

Reference: P. G. Harrison and W. J. Knottenbelt, “Passage Time Distributions in Large Markov Chains”, 2002, Sec. 3.2.

HMOM selects which of the paper’s two recursions runs, and they are NOT the same computation:

(nstates x NMAX) matrix m_i(r), the holding time in i depends only on i.

Eq. 7 with the u_i(r) recurrence of Eq. 8,
u_i(r) = -sum_{j=1..r} C(r,j) m_i(j) u_i(r-j), u_i(0) = 1,

which are the derivatives at the origin of 1/h*_i(s). Cheaper: no per-pair moments.

(nstates x nstates) cell HMOM{i,k} = [m_ik(1) … m_ik(NMAX)], the r-th
moment of the holding time in i WHEN THE NEXT STATE IS k. Eq. 6, the full Markov-renewal kernel. m_ik(0) = P(i,k) is implied and must not be supplied.

Both solve one linear system of the size of the non-target block per order, with the lower-order terms already known – the iteration the paper describes. Unlike the Markov case, the n-th moment needs every moment from 1 to n, so NMAX cannot be raised for free.

ctmc_passage_time(Q, pi0, target, tset, options)

[F, f, OUT] = CTMC_PASSAGE_TIME(Q, PI0, TARGET, TSET, OPTIONS)

Cumulative distribution F and density f of the first passage time from the initial law PI0 into the target state set, on the grid TSET.

F(t) = 1 - alpha exp(St) 1 f(t) = alpha exp(St) s0

OPTIONS.method selects the route:

‘expm’ (default) exact, one matrix exponential reused along a uniform: grid, in the manner of the CTMC response-time getters
‘lt’ the transform of Eqs. 1-2 inverted through api/lti;: OPTIONS.lti_method picks the inverter (‘euler’ default)

‘lt’ exists for chains whose non-target block is too large for a dense exp(St), not because it needs fewer time points: see CTMC_PASSAGE_LST. On a small chain ‘expm’ is both faster and more accurate, which is why it is the default.

OUT carries .atom (the mass of PI0 already inside the target, i.e. F(0)), .alpha, .S, .s0 and .keep.

Reference: P. G. Harrison and W. J. Knottenbelt, “Passage Time Distributions in Large Markov Chains”, 2002.

ctmc_passage_ph(Q, pi0, target)

[ALPHA, S, S0, KEEP, ATOM] = CTMC_PASSAGE_PH(Q, PI0, TARGET)

Phase-type representation of the first passage time from the initial law PI0 into the target state set TARGET, in the CTMC with generator Q.

This is the primitive behind the whole ctmc_passage_* family. With A the complement of TARGET,

S = Q(A,A) sub-generator: the passage has not completed s0 = -S*1 (= Q(A,TARGET)*1) exit vector alpha = PI0(A) UNNORMALIZED, see below atom = sum(PI0(TARGET))

so that L(s) = alpha (sI - S)^{-1} s0 + atom and F(t) = 1 - alpha exp(St) 1. Compare Eqs. 1-2 of P. G. Harrison and W. J. Knottenbelt, “Passage Time Distributions in Large Markov Chains”, 2002, which write the same system as n linear equations with L_i = 1 on the target.

ALPHA IS DELIBERATELY NOT NORMALIZED. Its mass is 1 - ATOM; the missing mass is the ATOM AT ZERO carried by initial states already in the target set. A caller that normalizes alpha and forgets the atom reports F(0) = 0 for a passage that has already completed with probability ATOM.

PI0 may be empty, in which case the conditional stationary law on A is used. TARGET is 1-based here and 0-based in the Java, Python and C++ twins.

KEEP maps the rows of S back to state indices of Q.