ctmc_takahashi.h File Reference

Takahashi's aggregation-disaggregation for a nearly completely decomposable CTMC. More...

#include <cstddef>
#include <vector>
#include "line/api/mc/ctmc_bicgstab.h"
#include "line/api/mc/ctmc_courtois.h"
#include "line/api/mc/ctmc_kms.h"
#include "line/api/mc/ctmc_randomization.h"
#include "line/api/mc/dtmc_solve.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_takahashi.h:

Go to the source code of this file.

Classes
struct	line::mc::TakahashiResult< T >

Namespaces
namespace	line
namespace	line::mc

Functions
template<class T>
TakahashiResult< T >	line::mc::ctmc_takahashi (const Matrix< T > &Q, const std::vector< std::vector< std::size_t > > &MS, std::size_t numSteps, double massTol=1e-14)
	Takahashi's aggregation-disaggregation for a nearly completely decomposable CTMC.

Detailed Description

Takahashi's aggregation-disaggregation for a nearly completely decomposable CTMC.

Templated port of matlab/src/api/mc/ctmc_takahashi.m and jar/src/main/java/jline/api/mc/Ctmc_takahashi.java. Like KMS it starts from the Courtois approximation and alternates aggregation with disaggregation, but the disaggregation step is a per-macro-state fixed point rather than a Gauss-Seidel sweep: for macro-state I,

(I - P_II') x_I = b_I, b_I(i) = sum_{K != I} gamma_K GI(K, i),

with gamma the aggregated macro-state distribution and GI(K, .) the conditional one-step flow out of macro-state K.

Everything here is in the ORIGINAL state ordering: P is the uniformization of the UNPERMUTED Q and each macro-state is addressed by the indices MS[I] carries. Contiguous block offsets would silently solve for a different partition of the same block sizes whenever a macro-state is not a contiguous range, which is the defect the MATLAB version records in its own comments.

The uniformization rate is passed explicitly as (21/20) max|Q|, matching the rate ctmc_courtois derives. The reference used to leave it to the ctmc_randomization default, max|Q| + rand, which made the whole iteration irreproducible for no benefit: the aggregation and disaggregation equations are homogeneous in P - I = Q/q, so the fixed point does not depend on q.

REFERENCE DEFECT, repaired here. MATLAB guards the aggregation against an empty macro-state (S > 1e-14) but then divides by that same S unguarded when building GI, so a macro-state carrying no mass poisons the whole iterate with NaN. Here a macro-state whose current mass is below the threshold contributes nothing to either G or GI, which is the limit of the expression as S -> 0.

GATED ON TRANSCENDENTAL ARITHMETIC: seeded by ctmc_courtois, whose epsMAX is an eigenvalue modulus, and its large-system block solve is GMRES.

Definition in file ctmc_takahashi.h.