Paper deep dive
Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)
Jianlin Li, Nick Guo, Peter Ye, Yizhou Zhang
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%
Last extracted: 8/4/2026, 4:17:51 AM
Summary
The paper introduces Tessa, a methodology for verifying finite-horizon Markov chains by casting probabilistic model checking as dense tensor computations. This approach leverages hardware accelerators like GPUs via off-the-shelf compiler toolchains (e.g., JAX/XLA), offering significant speedups over traditional sparse matrix or symbolic methods for dense transition dynamics.
Entities (10)
Relation Signals (9)
Tessa → implements → Tensor Probabilistic Model Checking
confidence 95% · We implement our approach in a tool called Tessa.
Tessa → uses → Dense Tensors
confidence 95% · Our insight is to cast probabilistic model checking of Markov chains as computations over dense tensors.
Markov Chain → verifiedby → Probabilistic Model Checking
confidence 92% · Probabilistic model checking is a formal verification technique for verifying quantitative properties of systems exhibiting stochastic behavior. Given a mathematical model (such as a Markov chain)...
Probabilistic Model Checking → verifies → Step-bounded Reachability
confidence 92% · verifying Markov chains with respect to step-bounded reachability probabilities.
Tessa → targets → GPU
confidence 90% · offload them to GPU devices for massively parallel execution.
Tessa → uses → JAX
confidence 88% · compiling DTMC models to dense tensor computations (i.e., array programs in an array programming language like JAX [9]).
Tessa → uses → XLA
confidence 88% · we can use an off-the-shelf tensor compiler (XLA [2], in this case) to optimize the generated tensor code
Tessa → outperforms → PRISM
confidence 85% · Empirical evaluation shows that Tessa unlocks massive speedups over state-of-the-art methods... compared to Storm
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We reexamine the problem of verifying Markov chains with respect to step-bounded reachability probabilities. Prevailing approaches rely on encoding the state-transition matrix using either explicit or symbolic representations. While these approaches are effective for sparse transition dynamics, they scale less favorably in the dense regime. Our insight is to cast probabilistic model checking of Markov chains as computations over dense tensors. This methodology enables the use of off-the-shelf compiler toolchains for optimized execution of these tensor computations on hardware accelerators. We prove the soundness of the methodology of mapping probabilistic model checking to tensor computations. We implement our approach in a tool called Tessa . Empirical evaluation shows that Tessa unlocks massive speedups over state-of-theart methods on selected benchmarks from the literature.
Tags
Links
- Source: https://arxiv.org/abs/2608.00374v1
- Canonical: https://arxiv.org/abs/2608.00374v1
Trouble viewing inline? Open PDF directly →
Full Text
70,968 characters extracted from source content.
Expand or collapse full text
Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version) Jianlin Li , Nick Guo , Peter Ye , and Yizhou Zhang ( B ) David R. Cheriton School of Computer Science University of Waterloo, Canada jianlin.li, nick.guo, p2ye, yizhou@uwaterloo.ca Abstract.We reexamine the problem of verifying Markov chains with respect to step-bounded reachability probabilities. Prevailing approaches rely on encoding the state-transition matrix using either explicit or sym- bolic representations. While these approaches are effective for sparse tran- sition dynamics, they scale less favorably in the dense regime. Our insight is to cast probabilistic model checking of Markov chains as computations over dense tensors. This methodology enables the use of off-the-shelf compiler toolchains for optimized execution of these tensor computations on hardware accelerators. We prove the soundness of the methodology of mapping probabilistic model checking to tensor compu- tations. We implement our approach in a tool called Tessa. Empirical evaluation shows that Tessa unlocks massive speedups over state-of-the- art methods on selected benchmarks from the literature. 1 Introduction Probabilistic model checking is a formal verification technique for verifying quan- titative properties of systems exhibiting stochastic behavior. Given a mathemat- ical model (such as a Markov chain) representing the behavior of a system over time, an algorithmic procedure determines whether the system satisfies proper- ties of interest (such as step-bounded reachability). As an example, consider the problem of scheduling an urgent faculty meeting (Figure 1) involvingNprofessors by collecting their availability via a Doodle poll. At any time step, each professoriis in one of three states:Away,Doodling, orDone. A professorAwayfrom their inbox may, with probabilityp i , notice the Doodle poll email and transition toDoodlingto answer the Doodle poll. However, theDoodlingstate is fragile: with probabilityq i , they successfully complete the poll and reachDone; but with probability1−q i , they are interrupted by a new email or a knock on their door, forcing them back to beAwaywithout hitting submit. The dynamics capture the intermittent nature of professors’ attention. Each Markov chain in Figure1(a)depicts the behavior of a professor. The quantitative property to verify is the probability that allNprofessors have reachedDonewithin a horizon ofHtime steps. ©2026 The Authors. Technical report extending the authors’ CAV 2026 paper [31] with an appendix. arXiv:2608.00374v1 [cs.LO] 1 Aug 2026 2Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang Away s 1 = 0 Doodling s 1 = 1 Done s 1 = 2 1−p 1 p 1 1−q 1 q 1 1 ‖·‖ Away s N = 0 Doodling s N = 1 Done s N = 2 1−p N p N 1−q N q N 1 (a)Nparallel Markov chains capturing the stochastic behavior ofNprofessors. moduleP 1 vars 1 : [0..2]init0; [a] (s 1 = 0)→p 1 : (s 1 : = 1) + (1−p 1 ) : (s 1 : = 0); [a] (s 1 = 1)→q 1 : (s 1 : = 2) + (1−q 1 ) : (s 1 : = 0); [a] (s 1 = 2)→1.0 : (s 1 : = 2); module·//modulesP 2 ·P N goals 1 = 2∧·∧s N = 2 (b) A model ofNprofessors (in PRISM-like syntax). 68101214 0 100 200 300 N 0200400 0 100 200 300 H Storm (MTBDD) Dice Storm (Sparse) Tessa time (s) (c) Scaling plots. Figure 1. Running example: collecting availability fromNprofessors. The target property is the probability that allNprofessors reachDonewithin horizonH. The program in Figure1(b)models this system of stochastic processes, sim- ilarly to how one would structure it in the PRISM [1] and JANI [10] modeling languages. Each moduleP i models the behavior of professori: it declares a state variables i with domain0,1,2(indicatingAway,Doodling, andDone, respectively), and threecommandsthat updates i according to the transition probabilitiesp i andq i . Each command isguardedby a Boolean expression (e.g., s i = 0for the first command); a command is enabled only if its guard is satisfied in the current state. Thegoalclause of the global model specifies the target property that all professors are in stateDone. Despite the semantic simplicity of this model, the verification problem is in- herently computationally expensive due to state explosion. Since each process has 3 states, the global state-space size is3 N . It is unlikely that any model check- ing tool on a classical computer can avoid this exponential blowup. Nevertheless, we would like to push the boundary of tractable instances as far as possible, especially given the raw speed of modern hardware accelerators at our disposal. State-of-the-art methods largely fall into two categories: –Probabilistic model checkers, such as PRISM [29] and Storm [16], are the pre- vailing tools for such verification tasks. They represent the state-transition dynamics using either explicit (i.e., sparse matrices of size3 N ×3 N in theN- professor model) or symbolic (i.e., multi-terminal binary decision diagrams, or MTBDDs) representations. Both representations are highly optimized for CPU execution in the situation of sparse transition dynamics. However, they are not as effective in the dense regime, and their irregular memory access patterns make it challenging to map them efficiently to modern hardware ac- celerators. Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)3 –A recent development is the use of probabilistic inference for probabilistic model checking. Rubicon [19] compiles a DTMC into the Dice [18] proba- bilistic programming language; running the Dice program using Dice’s infer- ence engine—weighted model counting on binary decision diagrams (BDDs)— computes the step-bounded reachability probability. Still, BDD operations are pointer-rich, so they involve indirection that does not map well to devices like GPUs. We present Tessa, a new methodology for verifying step-bounded reachability properties of DTMCs, by compiling DTMC models to dense tensor computations (i.e., array programs in an array programming language like JAX [9]). As Figure1(c)shows, Tessa exhibits markedly better scalability than existing tools on theN-professor model. The two plots show running times for different values ofN(fixing the horizonH= 10) and different values ofH(fixingN= 12), respectively. The first plot shows that, while Tessa does not change the asymp- totic complexity of the problem (the state-space size still grows exponentially withN), it aggressively pushes down the constant factor of the exponential growth compared to Storm and Dice. The second plot shows that while all tools appear to scale linearly withH, Tessa has a significantly smaller slope—the curve is nearly flat. 1 A key reason for Tessa’s scalability is its representation of state distributions as dense tensors, and state-transition dynamics as transformations on these ten- sors. The name Tessa evokes a tesseract—a high-dimensional cube—reflecting our insight: rather than eagerly materializing the state-transition dynamics as a sparse matrix or a decision diagram, we treat the joint distribution ofNstate variables as an order-Ntensor, and the transition dynamics as a tensor program. A key advantage of this tensor representation is its amenability to optimiza- tion (by modern tensor compilers) and acceleration (by hardware such as GPUs). Because Tessa maps the verification problem to tensor operations, we can use an off-the-shelf tensor compiler (XLA [ 2], in this case) to optimize the generated tensor code 2 and offload them to GPU devices for massively parallel execution. In contrast, established methodologies target representations that induce irregu- lar memory accesses, which limit both the degree of parallelism and the ease of parallelization achievable. Targeting tensors also naturally supports parameter search for DTMCs with unknown transition probabilities. As Tessa generates differentiable tensor pro- grams, we can leverage JAX’s automatic differentiation to compute gradients of a rich selection of optimization objectives. Furthermore, the computationally in- tensive gradient calculation can be optimized by XLA and accelerated by GPUs. We believe the methodology embodied by Tessa usefully complements the toolbox of probabilistic model checkers. It is not a silver bullet; existing methods are already effective for sparse models. However, Tessa opens up a new regime of tractable models that were previously out of reach. 1 The comparison may make state-of-the-art tools look inefficient, but they are actually highly optimized for CPUs. The point is that a methodological shift unlocks massive speedups. 2 Tessa running times in Figure1(c)include compilation times of JAX and XLA. 4Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang Contributions.Our contributions are both theoretical and practical: –We base our development on a core language for specifying DTMCs (Section3). We show how to interpret programs in this core language as computations over tensors (Section4). This interpretation captures the essence of the implemen- tation of Tessa. We prove that this tensor interpretation is sound with respect to the verification of step-bounded reachability properties (Section5). –We implement Tessa and evaluate it on selected benchmarks from the liter- ature. Tessa employs a domain-specific compiler stack that automatically opti- mizes and parallelizes dense tensor computations for GPU execution (Section6). Experimental results highlight substantial performance gains unlocked by Tessa’s tensor-based approach (Section7). 2 Preliminaries We review discrete-time Markov chains (DTMCs) and the verification problem of step-bounded reachability properties. AMarkov chainis a tupleM= (S, ι, η,G): –Sis a finite set of states. –ι∈D(S)is the initial-state distribution. –η:S →D(S)is the transition function. –G ⊆ Sis the set of goal states. We useD(S)to denote the set of probability distributions overS. We write ι(s)to denote the probability of starting in states. We writeη(s, s ′ )to denote the probability of transitioning to states ′ from states. Thus, we have that ∑ s∈S ι(s) = 1and also that for any states∈ S, ∑ s ′ ∈S η(s, s ′ ) = 1. In this paper, we are interested in the verification of step-bounded reachabil- ity properties of the formPr M ( ♢ ≤n G ) , which denotes the probability of reaching a goal state withinnsteps. The probability of reaching a goal state withinnsteps starting from states is denotedPr M ( s⊨ ♢ ≤n G ) and is defined inductively as follows: Pr M ( s⊨ ♢ ≤0 G ) : = 1ifs∈ G 0otherwise (2.1) Pr M ( s⊨ ♢ ≤n+1 G ) : = 1ifs∈ G ∑ s ′ ∈S η(s, s ′ )·Pr M ( s ′ ⊨ ♢ ≤n G ) otherwise (2.2) The inductive case (2.2) is intuitive: the probability of reaching a goal state withinn+ 1steps is1if we are already in a goal state; otherwise, it is the weighted sum of the probabilities of reaching a goal state withinnsteps from each possible next state. ThenPr M ( ♢ ≤n G ) , the probability of reaching a goal state withinnsteps when starting from a state drawn from the initial-state distribution, is defined as Pr M ( ♢ ≤n G ) : = ∑ s∈S ι(s)·Pr M ( s⊨ ♢ ≤n G ) .(2.3) Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)5 We introduce another (perhaps less common) way to calculate step-bounded reachability probabilities. We writePr M ( μ ▶ ♢ ≤n G ) to denote the probability of reaching a goal state withinnsteps when starting from a state drawn from distributionμ∈D(S). It is defined inductively as follows: Pr M ( μ ▶ ♢ ≤0 G ) : = ∑ s∈G μ(s)(2.4) Pr M ( μ ▶ ♢ ≤n+1 G ) : = ( ∑ s∈G μ(s) ) + Pr M ( ∑ s /∈G μ(s)η(s) ▶ ♢ ≤n G ) (2.5) The inductive case (2.5) is intuitive: the probability of reaching a goal state withinn+ 1steps is the sum of (1) the probability of being in a goal state at the start, and (2) the probability of reaching a goal state withinnsteps after taking one step from a non-goal state. Notice that in the second term, the next- state distribution is given by ∑ s /∈G μ(s)η(s), which is the weighted sum of the transition distributions from all non-goal states. The correctness of this alternative way to calculate step-bounded reachability probabilities is established by Theorem2.1(and proven in the accompanying technical report [33]): Theorem 2.1.Pr M ( ♢ ≤n G ) = Pr M ( ι ▶ ♢ ≤n G ) . Its connection to tensor probabilistic model checking will become clear later, but for now we simply note that we usePr M ( μ ▶ ♢ ≤n G ) as a bridge to connect tensor-based reachability to the standard definition of reachability probabilities. 3 A Core Model-Specification Language Directly defining a DTMC via the four-tuple formalism in Section2is intractable for complex systems. The state space typically grows exponentially with the num- ber of variables, making such brute-force specification error-prone and difficult to maintain. As a result, probabilistic model checkers such as PRISM [29] and Storm [ 16] adopt high-level modeling languages that allow users to specify the system in a modular fashion. In this section, we define a high-level, modular specification language that captures the core aspects of a DTMC modeling language (e.g., thedtmcdialect of the PRISM language [ 1]). This allows describing a system concisely as a collection of interacting modules rather than a monolithic state-transition matrix. We call this core language PML. To place its tensor interpretation on a formal footing (Sections 4and5), we first define PML’s syntax in Section3.1and its Markov-chain semantics in Section3.2. 3.1 Syntax of PML Figure2presents the syntax of PML. 6Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang ModelM :: =modelm 1 ;·;m K goale Modulem :: =moduleDC DeclarationsD :: =d 1 ;·;d n Declarationd :: =varx:[0..n max ] initn init CommandsC :: =c 1 ;·;c n Commandc :: =[a]g→U Guardg :: =e Mixture of updatesU :: =θ 1 :u 1 +·+θ n :u n Updateu :: =x 1 : =e 1 ·x n : =e n Expressione :: =n|x|op(e 1 ,...,e n ) n∈Nx∈Variablesθ∈[0,1] Figure 2. PML syntax. A modelMcontainsKmodules. Each modulem k (1≤k≤K) declares a disjoint set of state variablesX k via a set of declarationsD k . The module further specifies its behavior via a set of commandsC k that update the state variables inX k . While each module declares and updates only its own state variables, it can read the variables of all the modules inM. A declarationddefines a state variablex, its domain (a bounded range of non-negative integers), and its initial value. A commandcconsists of an action labela, a guardg, and a mixtureU of updates. The action label allows commands from different modules to be synchronized. The guardgis a Boolean expression specifying a necessary (though not sufficient) condition for the command to be enabled. A command may also be unlabeled (written[]in PRISM), in which case it is assigned a globally unique action label. As a simplification, we assume that all unlabeled commands have already been assigned globally unique action labels. Each updateu i inUis associated with a probabilityθ i , where ∑ i θ i = 1. Each update further consists of deterministic assignments to a subset of the module’s declared state variables. An expressioneis either a constantn, a state variablex, or ann-ary operation over sub-expressionse 1 , . . . , e n . As a simplification, we treat Boolean expressions as integer expressions where0represents false and non-zero values represent true. The global modelMalso specifies the goal states via a Boolean expressione over the state variables of all modules. 3.2 Semantics of PML We now define the standard semantics of PML by interpreting a modelMas a DTMCM= (StateJMK,InitJMK,StepJMK,GoalJMK). LetX k be the set of state variables declared in modulem k (1≤k≤K). Let Vars(M) : = ⊎ K k=1 X k be the set of all state variables inM. Letdom(x)denote Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)7 the domain of variablex. ForX⊆Vars(M), the state space formed by the variables inXis State(X) : = (x∈X)→dom(x).(3.1) That is, a state is a mapping from each variable inXto a value in its domain. The state space of the entire modelMis then defined as StateJMK : =State(Vars(M)).(3.2) We will writeStateas a shorthand forStateJMKwhenMis clear from context. Notice thatState(X)is measurable since it is a finite set. The size ofState(X) is|State(X)|= ∏ x∈X |dom(x)|, which grows exponentially with the number of state variables inX. We writes(x)to denote the value of variablexin states. Lets 0 ∈Statebe the initial state where each variablex∈Vars(M)is initial- ized to the value specified in its declaration. Then the initial-state distribution is InitJMK : =δ s 0 (3.3) whereδ s is the delta distribution (point mass) ats. Letebe the goal expression ofM. Then the set of goal states is defined as GoalJMK : = s∈State ∣ ∣ ∣ JeK(s)6= 0 (3.4) whereJeKsis the interpretation of expressioneunder states. All that remains is to defineJeKandStepJMK. To defineStepJMK, we need to first define the semantics of each modulem k inM, which further requires defining the semantics of each command inm k and each update in a command. Interpreting expressions.JeK:State(X)→N, whereXis a superset of the variables appearing ine, is defined inductively as follows: JnKs : =nJxKs : =s(x) Je i Ks=n i fori= 1, . . . , lJopK=f Jop(e 1 ,· · ·, e l )Ks : =f(n 1 , . . . , n l ) Interpreting updates.Recall thatX k denotes the set of state variables de- clared inm k . LetY k be the set of variables each of which is either declared in m k or read inm k ’s commands. So we have thatX k ⊆Y k ⊆Vars(M). An updateuin modulem k is interpreted as a functionJuK:State(Y k )→ State(X k ). That is, fors∈State(Y k ), Jx 1 : =e 1 · · ·x n : =e n Ks : =s| X k [x 1 7→Je 1 Ks,· · ·, x n 7→Je n Ks]. (3.5) In (3.5), the notations| X k denotes the restriction ofsto the variables inX k , ands| X k [x 1 7→Je 1 Ks,· · ·, x n 7→Je n Ks]denotes the state obtained by updating s| X k with the assignments inu. 8Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang Interpreting mixtures of updates.A mixture of updatesUin modulem k is interpreted as a functionJUK:State(Y k )→D(State(X k )). That is, Jθ 1 :u 1 +· · ·+θ n :u n Ks : = n ∑ i=1 θ i ·δ Ju i Ks .(3.6) Interpreting modules.LetAction(M)be the modelM’s alphabet of action labels. LetAction(m k )be the set of action labels appearing inm k ’s commands. The modulem k is interpreted as a functionJm k K:Action(M)→State(Y k )→ D <∞ (State(X k )). Here,D <∞ (State(X k ))is the set of unnormalized distributions overState(X k ); these distributions have finite total mass, but the total mass is not necessarily1. Specifically, Jm k Ka s : = δ s| X k ifa6∈Action(m k ), ∑ [a]g→U∈Com(m k ,a) [JgKs6= 0]·JUKsotherwise. (3.7) In (3.7),Com(m k , a)is the set of commands in modulem k with action labela, and[P]is the indicator function that evaluates to1if predicatePis true and0 otherwise. Intuitively,Jm k Ka sdescribes the aggregate transition behavior of modulem k when the current state issand the chosen action isa. There are two cases. If ais not inm k ’s alphabet, thenm k stays put. Otherwise, a sum ranges over all commands with action labelawhose guard is satisfied bys. When exactly one such command[a]g→Uexists,Jm k Ka sis the distribution denoted by JUKs. When multiple commands with the action labelahave overlapping guards, Jm k Ka sis an unnormalized distribution whose total mass equals the number of enabled commands; PRISM allows this ambiguity and resolves it by a random choice among all the enabled command combinations (see ( 3.10)). If no command with action labelais enabled ins,Jm k Ka shas total mass0. Interpreting models.In a specification language like PRISM [1], the global behavior of a DTMC model is defined by the parallel execution of modules synchronized by action labels. A global transition step corresponds to the simul- taneous transitioning of all modules according to a common action label. LetEnJm k Ka s, wheres∈State(Y k ), count the number of commands with action labelathat are locally enabled in modulem k in states: EnJm k Ka s : = 1ifa6∈Action(m k ), ∑ [a]g→U∈Com(m k ,a) [JgKs6= 0]otherwise, (3.8) In casea /∈Action(m k ), it is considered that there is one enabled command—the implicit self-loop—soEnJm k Ka sis defined to be1. LetEnJMKa s, wheres∈State, count the total number of enabledcommand combinationsfor actionaacross all modules: EnJMKa s : = K ∏ k=1 EnJm k Ka s| Y k . (3.9) Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)9 Now we can define the DTMC transition functionStepJMK:State→D(State): StepJMKs : = ∑ a ⊗ K k=1 Jm k Ka s| Y k ∑ a EnJMKa s if ∑ a EnJMKa s >0, δ s otherwise. (3.10) The cases in (3.10) depend on whether there is at least one enabled command combination in states. –In PRISM, multiple commands—possibly with different action labels, or with the same label but overlapping guards—may be simultaneously enabled in a state for a module. This ambiguity is resolved by a uniformly random choice over all enabled command combinations across all modules. A command combination for actionaconsists of one enabled command from each module; the number of such command combinations isEnJMKa s= ∏ k EnJm k Ka s| Y k . The denominator ∑ a EnJMKa sin (3.10) is the total num- ber of enabled command combinations across all actions. The operator ⊗ in (3.10) denotes the product of (unnormalized) distributions over disjoint sets of variables. Since theKdisjoint sets of variablesX 1 , . . . , X K aggregate toVars(M), ⊗ K k=1 Jm k Ka s| Y k is an (unnormalized) distribution over the global state spaceState. Moreover, its total mass equalsEnJMKa s. –If no action is enabled ins, the DTMC loops back tos. While ⊗ K k=1 Jm k Ka s| Y k is in general an unnormalized probability distribu- tion,StepJMKsis a probability distribution. We prove this fact (Theorem 3.1) in the technical report [33]. Theorem 3.1.For anys∈State, StepJMKs∈D(State). 4 Probabilistic Model Checking as Tensor Computations In this section, we show how to cast the verification of PML models as compu- tations over tensors. Section4.1maps state distributions to tensors. Section4.2 maps PML models to tensor transformers. Section4.3maps the verification of step-bounded reachability properties to tensor computations. 4.1 Representing State Distributions as Tensors We use tensors to represent discrete probability distributions over DTMC states. A tensor is a multidimensional array generalizing scalars (order-0), vectors (order- 1), matrices (order-2), and so on. In this work, we use tensorsT∈R a 1 ×·×a N over the field of real numbersRand tensorsT∈N a 1 ×·×a N over the field of natural numbersN, whereNis theorder(a.k.a.rank) of the tensor anda k is the size of thek-th dimension. 10Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang The state space as an index space for tensors.A state spaceState(X) = (x∈X)→dom(x)serves as the index space for order-|X|tensors: –each dimension corresponds to a state variablex∈X, –and the size of each dimension is|dom(x)|. We write Tensor R ( State ( X )) to denote the set of R -valued order- | X | tensors with the index spaceState(X). That is,Tensor R (State(X)) : =R ∏ x∈X |dom(x)| . ForT∈Tensor R (State(X)),T[s]denotes the entry ofTat indexs∈State(X). Tensor R (State(X))is avector space(a.k.a.linear space) under element-wise addition and scalar multiplication. This|X|-dimensional vector space has a nat- ural set of basis vectors B s ∣ ∣ ∣ s∈State(X) , where each basis vectorB s ∈ Tensor R (State(X))is a tensor defined asB s [s ′ ] : = [s ′ =s]. State distributions as tensors.We define a functionΞ:D <∞ (State(X))→ Tensor R (State(X))mapping a distributionμ∈D <∞ (State(X))to a tensor: Ξ(μ) : = ∑ s∈State(X) μ(s)·B s . That is,Ξ(μ) =Tiff∀s∈State(X),T[s] =μ(s). For example, given a joint distributionμof two state variables (say, the states of professors in Figure1 whenN= 2),Ξmapsμto an order-2 tensorT∈R 3×3 : distribution μ s 1 s 2 probability 0 00.1 0 10.2 0 20.0 1 00.3 · · 0.10.20.0 0.30.00.1 0.00.10.2 s 1 =0 s 1 =1 s 1 =2 s 2 =0 s 2 =1 s 2 =2 tensor T Ξ(μ) =T Each dimension corresponds to a professor’s state. Both dimensions have size 3 (corresponding to statesAway,Doodling,Done). This tensor should not be con- fused with the state-transition matrix used by a probabilistic model checker like Storm [ 16], which would be of size9×9(since there are 9 joint states). Index tensors.For each variablex∈X, we define itsindex tensorI x∈X ∈ Tensor N (State(X))as aN-valued order-|X|tensor. For any indexs∈State(X), I x∈X [s] =s(x)is simply the value of that variable in states. 4.2 Tensor Semantics of PML Models Interpreting expressions.Section3.2interprets an expressioneas a function JeK:State(X)→N. We now defineLeM X ∈Tensor N (State(X))as follows, where the subscriptXindicates the shape of the tensor: LnM X : =n·1 X LxM X : =I x∈X Le i M X =T i fori= 1, . . . , lLopM X =f Lop(e 1 ,· · ·, e l )M X : =f(T 1 , . . . ,T l ) where1 X ∈Tensor N (State(X))is the all-ones tensor, andLopM X isJopKlifted to operate element-wise over tensors. Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)11 Interpreting updates.Section3.2interprets an updateuin modulem k as a functionJuK:State(Y k )→State(X k ). So givens∈State(Y k ),JuKsgives the value ofx∈X k after executinguin states. This reading suggests that we can de- fine the tensor interpretation ofuas a functionLuM Y k :X k →Tensor N (State(Y k )): Lx 1 : =e 1 · · ·x n : =e n M Y k x : = Le i M Y k ifx=x i for somei, I x∈Y k otherwise. (4.1) Interpreting mixtures of updates.Section3.2interprets a mixture of up- datesUin modulem k as a functionJUK:State(Y k )→D(State(X k )). We now interpretUas a tensorLUM∈Tensor R (State(Y k ]X k )), where]is the disjoint- union operator. The tensorLUMis of order|Y k |+|X k |. Specifically, Lθ 1 :u 1 +· · ·+θ n :u n M : = n ∑ i=1 θ i · ⊙ x∈X k [( Lu i M Y k x ) ⊗1 X k =1 Y k ⊗I x∈X k ] , (4.2) where ⊙ is the Hadamard (i.e., element-wise) product,⊗is the outer product, and=is overloaded to denote element-wise equality. Observe the correspondence between (3.6) and (4.2). The outer products in (4.2) are used to align the tensor dimensions. Intuitively, if we write(s, s ′ )as an index, wheres∈State(Y k )and s ′ ∈State(X k ), then we should haveLθ 1 :u 1 +· · ·+θ n :u n M[(s, s ′ )] = ∑ n i=1 θ i · [Ju i Ks=s ′ ]. That is, the tensor’s value at index(s, s ′ )is the probability of tran- sitioning to states ′ from stateswhen the mixture of updatesUis executed. Interpreting modules.Section3.2interprets a modulem k as a functionJm k K: Action(M)→State(Y k )→D <∞ (State(X k )). Accordingly, we define the tensor interpretationLm k M:Action(M)→Tensor R (State(Y k ]X k ))as follows: Lm k Ma : = ⊙ x∈X k [ I x∈Y k ⊗1 X k =1 Y k ⊗I x∈X k ] ifa6∈Action(m k ), ∑ [a]g→U∈Com(m k ,a) ([ LgM Y k 6= 0 ] ⊗1 X k ) LUMotherwise. (4.3) Observe the correspondence between (3.7) and (4.3). Intuitively,Lm k Mais a ten- sor whose value at index(s, s ′ )is the unnormalized probability of transitioning to states ′ from stateswhenm k is executed with actiona. When multiple com- mands have overlapping guards, the sum ∑ s ′ Lm k Ma[(s, s ′ )]equals the number of enabled commands in states; normalization to a probability distribution occurs at the model level rather than at the module level (see ( 4.6)). Interpreting models.Section3.2interprets a modelMas a functionStepJMK: State→D(State). Following the way that modules are interpreted as tensors, we could interpretMas a tensorLMM:Tensor R (State(Vars(M)]Vars(M))). How- ever, this tensor would be of order2|Vars(M)|, effectively materializing the full state-transition matrix. For space efficiency, we instead interpretMas a function 12Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang LMM:Tensor R (State)→Tensor R (State), representing the transition dynamics as a tensor transformer (i.e., code) rather than as a tensor (i.e., data). Corresponding toEnJm k Ka sandEnJMKa sin (3.8) and (3.9), we define tensorsEnLm k Ma∈Tensor N (State(Y k ))andEnLMMa∈Tensor N (State): EnLm k Ma : = 1 Y k ifa6∈Action(m k ), ∑ [a]g→U∈Com(m k ,a) [ LgM Y k 6= 0 ] otherwise, (4.4) EnLMMa[s] : = K ∏ k=1 EnLm k Ma [ s| Y k ] (4.5) DefineL : = ∑ a EnLMMaas the tensor counting the total number of enabled command combinations in each state. We defineLMMby specifying how it transforms a tensorT∈Tensor R (State) representing the current-state distribution into a tensor representing the next- state distribution. Specifically, for anys ′ ∈State,LMMTweights the probability of transitioning tos ′ from each statesby the probabilityT[s]of being ins: (LMMT)[s ′ ] : = ∑ s∈State T[s]· ∑ a ∏ K k=1 Lm k Ma [( s| Y k ,s ′ | X k )] L[s] ifL[s]>0, [s=s ′ ]otherwise. (4.6) The transition probability is given by a case analysis, similarly to (3.10). The branching control flow in (4.6) hinders parallelization over the index spaceState, however. Fortunately, we can encode the branching logic as pure tensor flow, through a sum of two terms. We redefineLMMas follows: LMMT : = ( ∑ a P a ) + ( T [ L=0 Vars(M) ]) (4.7) P a [s ′ ] : = ∑ s∈State ( T max ( L,1 Vars(M) ) ) [s]· K ∏ k=1 Lm k Ma [( s| Y k , s ′ | X k )] (4.8) The two terms in (4.7) correspond to the two branches of (4.6). –The first term ∑ a P a rearranges the first branch of ( 4.6): it pushes ∑ s inside ∑ a and weights the input tensorTby1 max(L,1), as required by the averaging semantics. In particular, whenL[s] = 0, the normalization factor max(L[s],1) = 1, whileLm k Ma [ (s| Y k ,·) ] is zero for some modulekfor every actiona, so the term vanishes. Notice thatP a as defined in ( 4.8) is a tensor contraction over dimensionss corresponding to the current state. The output dimensionss ′ correspond to the next state. The tensors in this contraction are the masked input tensor and theKmodule tensorsLm k Ma. This tensor contraction is the main work performed by the tensor transformerLMM. –The second termT [L=0]weights the input tensorTby the mask[L=0]. This mask is1exactly where no actions are enabled, leaving the probability mass in those states unchanged and zeroing out the mass in all other states. Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)13 111 111 110 mask∆ ¬e ¬(s 1 =2∧s 2 =2) T 00 T 01 T 02 T 10 T 11 T 12 T 20 T 21 T 22 current tensorT n+1steps remain = T 00 T 01 T 02 T 10 T 11 T 12 T 20 T 21 0 ∆ ¬e T T ′ 00 T ′ 01 T ′ 02 T ′ 10 T ′ 11 T ′ 12 T ′ 20 T ′ 21 T ′ 22 next tensorLMM(∆ ¬e T) nsteps remain LMM apply tensor transformer T 22 , goal state mass inT, is moved to the accumulated prob- ability of reaching the goal. Figure 3. Visualizing the computation in (4.10) of the reachability probability. The probability massT 22 at the goal state (s 1 = 2, s 2 = 2) is extracted and ac- cumulated to the reachability probability. The remaining mass is transformed by the model’s tensor-transformer interpretation to produce the next-state tensor. 4.3 Casting Probabilistic Model Checking as Tensor Computations With the tensor-transformer interpretation defined, we can now cast probabilistic model checking as tensor computations. In words, step-bounded reachability probabilities can be computed via repeated applications of the tensor transformer LMM. Specifically, for modelMwith goal expressione, we writePr M ( T ♢ ≤n e ) to denote the probability of reaching a state satisfying expressionewithinn steps when starting from a state drawn from the distribution represented by tensorT∈Tensor R (State). It is defined inductively as follows, where∆ e : = [ LeM Vars(M) 6=0 Vars(M) ] and∆ ¬e : = [ LeM Vars(M) =0 Vars(M) ] are mask tensors: Pr M ( T ♢ ≤0 e ) : =〈∆ e ,T〉(4.9) Pr M ( T♢ ≤n+1 e ) : =〈∆ e ,T〉+ Pr M ( LMM(∆ ¬e T)♢ ≤n e ) (4.10) The definition mirrors that ofPr M ( μ ▶ ♢ ≤n G ) in (2.4) and (2.5). In the base case (4.9), the probability of being in a goal state is given by the Frobenius inner product〈∆ e ,T〉= ∑ s∈State ∆ e [s]·T[s]. In the inductive case (4.10), the next- state tensor is given byLMM(∆ ¬e T), which is the result of applying the tensor transformerLMMto the current-state tensor masked by∆ ¬e . The mask tensor ∆ ¬e ensures that only the non-goal states in the current-state tensor contribute to the next-state tensor. Figure 3illustrates (4.10) for the 2-professor example. 4.4 Discussion The correctness ofPr M ( T♢ ≤n e ) with respect to the verification problem stated in Section2will be established in Section5. The close correspondence between the definitions ofPr M ( μ ▶ ♢ ≤n G ) andPr M ( T♢ ≤n e ) makes short work of proving the correctness of tensor probabilistic model checking. Two factors contribute to the efficiency of tensor probabilistic model checking. One factor is that the computations in ( 4.7), (4.8), (4.9), and (4.10) avoid mate- rializing the full3 N ×3 N transition matrix of the DTMC (using theN-professor example for concreteness). Instead, the transition matrix is implicitly encoded 14Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang as tensor computations (i.e., code rather than data) that transform the order-N tensors. A more important factor contributing to the efficiency is that the tensor com- putations in (4.7), (4.8), (4.9), and (4.10) are composed of standard operations over dense tensors and, therefore, can be implemented as first-order array pro- grams in an array programming language such as JAX [9]. In other words, we have essentially compiled the probabilistic model checking problem for DTMCs into array programs of the kind that are otherwise ubiquitous in machine learn- ing. These array programs can then be optimized by machine-learning compilers and exploit the massive parallelism offered by hardware accelerators such as GPUs. Finally, compiling to tensor computations makes it possible to search for model parameters satisfying reachability properties, via gradient-based optimiza- tion. Since the generated array programs are differentiable, we can compute gradients of distributional properties with respect to model parameters using automatic differentiation, which is readily supported by JAX. 5 Correctness of Tensor Probabilistic Model Checking In this section, we establish the correctness of tensor probabilistic model check- ing. In the following, letMbe a PML model, and letμ∈D(StateJMK). Theo- rem5.1establishes the correctness of the tensor-transformer interpretationLMM with respect to the standard interpretationStepJMK. Theorem 5.1.LMMΞ(μ) =Ξ( ∑ s μ(s)·StepJMKs).That is, the diagram below commutes. M ∑ s μ(s)·StepJ·Ks Ξ(·) L · M Ξ ( μ ) LetM= (StateJMK,InitJMK,StepJMK,GoalJMK)be the DTMC represented byMper Section 3.2. Letebe the goal expression specified inM. LetG : = GoalJMK= s∈StateJMK ∣ ∣ ∣ JeK(s)6= 0 be the set of goal states inM. Theo- rem5.2follows from Theorem5.1. Theorem 5.2.Pr M ( Ξ(μ)♢ ≤n e ) = Pr M ( μ ▶ ♢ ≤n G ) for anyn∈N. Proofs of Theorem5.1and Theorem5.2are available in the technical report [33]. Theorem 5.3establishes the ultimate correctness of tensor probabilistic model checking with respect to step-bounded reachability probabilities. Theorem 5.3.Pr M ( Ξ(InitJMK) ♢ ≤n e ) = Pr M ( ♢ ≤n G ) for anyn∈N. Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)15 Proof.Letι=InitJMKbe the initial-state distribution of the DTMCM. Start- ing from the left-hand side: Pr M ( Ξ(ι)♢ ≤n e ) = Pr M ( ι ▶ ♢ ≤n G ) (by Theorem5.2) = Pr M ( ♢ ≤n G ) (by Theorem2.1) The proof of Theorem5.3reveals thatPr M ( ι ▶ ♢ ≤n G ) bridges tensor-based reachability and the standard definition of reachability probability. 6 Accelerating Tensor Computations with JAX and XLA At this point, we have reduced the verification of DTMCs (with respect to step- bounded reachability properties) to dense tensor computations. But that alone does not guarantee efficiency. Efficiency hinges on how fast the tensor computations can run. In the last decade, training and inference in machine learning (ML) have driven significant advances in compiler optimizations for tensor computations, as well as hardware acceleration for them. It is thus natural to leverage off-the-shelf ML compilers and hardware accelerators to speed up tensor probabilistic model checking. Specifically, Tessa compiles DTMC models to array programs in JAX [ 9]. JAX, embedded in Python, is a popular array programming language for high- performance numerical computing and machine learning. A just-in-time com- piler traces the JAX program. The resulting intermediate representation is then handed off to XLA (Accelerated Linear Algebra) [ 2], a domain-specific compiler designed to optimize tensor computations. XLA performs whole-program optimizations that are critical for performance on modern hardware accelerators. Importantly, it applies kernel fusion, merging multiple element-wise operations (such as those in Sections 4.2and4.3) into a single GPU kernel. Fusion significantly reduces memory footprint and the pressure on memory bandwidth, as it avoids materializing intermediate results between operations to the GPU memory. Such optimizations enable our tensor- based verification algorithm to fully saturate the massive parallelism offered by GPUs. Moreover, by compiling to JAX, Tessa leverages ML compiler optimizations transparently. The user does not write any GPU kernel code or manage GPU memory explicitly. Future improvements in ML compiler technology and hard- ware accelerators will directly benefit Tessa without changes to its implementa- tion. 7 Evaluation We have implemented Tessa in Python. In this section, we evaluate Tessa on two fronts: model checking and parameter synthesis. 16Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang 7.1 Model Checking Evaluation methodology.We adopt benchmarks directly from Rubicon [19]. These benchmarks consist of DTMC models with dense transition dynamics. Future work could expand the scope of the evaluation to include a more com- prehensive set of benchmarks. Nevertheless, the current selection is already rep- resentative of challenging Markov chain verification tasks in the dense regime. All experiments ran on a machine with an Intel Core i7-7820X CPU,128 GB RAM, and an NVIDIA GeForce RTX 2080 Ti GPU (11 GBVRAM). At first glance, comparing CPU-bound methods against a GPU-accelerated tool might seem like comparing apples to oranges. However, this hardware distinction is pre- cisely the point. Prevailing methods rely on representations that induce irregular memory access. Consequently, they do not map naturally to GPUs. In contrast, Tessa generates dense tensor workloads at which GPUs excel. We are therefore comparingmethodologieson the hardware they naturally map to—rather than comparing the hardware per se. We note that the GPU used in our experiments is an older consumer-grade model in the hardware vendor’s lineup. While we already observe substantial speedups over state-of-the-art tools, access to frontline hardware is expected to yield even greater gains. Baseline methods.For Rubicon [19], the authors use Storm (Sparse) and Storm (MTBDD) as baselines in their evaluation. We include both, as well as Dice (via the Rubicon transpiler [19]). –Storm (Sparse). In this engine of the Storm model checker, the transition dy- namics of a DTMC model is represented as a sparse matrix in the standard compressed sparse row(CSR) format. The engine uses efficient sparse matrix– vector multiplication kernels from off-the-shelf libraries such as Eigen and Gmm++. –Storm (MTBDD). This engine of Storm represents the DTMC transition dy- namics using multi-terminal binary decision diagrams. Random access to en- tries of the transition matrix is not as efficient as in the sparse-matrix engine, but the MTBDD engine can be more memory-efficient for certain models. –Dice. In this approach, the DTMC model specification is first lowered to the Dice probabilistic programming language [18] using the Rubicon tran- spiler [19]. The Dice compiler then represents the set of paths from the initial state to the goal states effectively as a BDD. This approach is shown to excel on models with certain structures. All methods are configured to use the double-precision floating-point format. Results.We now report results for each of the benchmarks on which the Rubi- con/Dice method is evaluated [19]:Queues,Weather Factories, andHerman. Queues.TheQueuesmodel consists ofKqueues, each with capacity3. Tasks arrive probabilistically at every step. Three queues are designated type 1, while the remainder are type 2. The goal states are those in which all type 1 queues and at least one type 2 queue are full. We compute the probability of reaching Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)17 46810 0 50 100 150 K(H= 10) 46810 0 2 4 6 8 K(H= 10) 0200400 0 20 40 60 H(K= 9) 0200400 0 0.5 1 H(K= 9) Storm (MTBDD) Dice Storm (Sparse) Tessa (1st) Tessa (2nd) time (s) Figure 4. Scaling on theQueuesbenchmark. 1015 0 200 400 N(H= 10) 1015 0 2 4 N(H= 10) 0200400 0 200 400 H(N= 13) 0200400 0 1 2 H(N= 13) Storm (MTBDD) Dice Storm (Sparse) Geni Tessa (1st) Tessa (2nd) time (s) Figure 5. Scaling on theWeather Factoriesbenchmark. 591317 0 50 100 150 N(H= 100) 591317 0 1 2 N ( H = 100 ) 0200400 0 200 400 H(N= 17) 0200400 0 0.5 1 1.5 H(N= 17) Storm (MTBDD) Dice Storm (Sparse) Tessa (1st) Tessa (2nd) time (s) Figure 6. Scaling on theHermanbenchmark. 18Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang a goal state within horizonH. Figure4shows how each tool scales withK(left two plots) and withH(right two plots). The first plot compares all methods asKvaries, withH= 10fixed. All methods scale exponentially withK. On the hardest instance that any base- line method can solve (K= 10), Tessa shows over 100×speedup over the next fastest method. The second plot zooms in on the performance of Tessa. Since JAX and XLA perform just-in-time (JIT) compilation, we measure two runs to show the effect of JIT compilation: in the second plot, the difference between the two Tessa curves indicates the JIT compilation overhead. 3 The right two plots show how the methods scale asHvaries, withK= 9fixed. Storm (MTBDD) and Dice reach time limits at lowerKvalues. Storm (Sparse) and Tessa scale linearly withH, but Tessa is ~40×faster atH= 500. Weather Factories.This model consists ofNfactories, each in a binary state: striking or operational. Transition probabilities are local but conditioned on a Markov process, weather. We verify the reachability of the state where all factories are simultaneously striking within a given horizonH. Figure 5shows how each method scales on this model. We additionally in- clude the Geni probabilistic programming language (PPL) as a baseline; the benchmark is part of its evaluation suite [34]. The methodology of using Geni for DTMC model checking is similar to that of Dice: both repurpose a PPL for DTMC model checking. The difference is that Geni’s compiler targets generating functions, while Dice’s compiler targets BDDs. All methods scale exponentially withN. Unlike in theQueuesbenchmark, here Dice scales better than Storm (Sparse) asNincreases. On the hardest instance that any baseline method can solve (N= 16), Tessa shows over 100× speedup over the next fastest method, Geni. All methods (that run atN= 13) scale linearly withH. AtH= 500, Tessa shows over 100×speedup over the next fastest method, Storm (Sparse). Herman.Herman’s protocol [17] is a well-known example in the literature of probabilistic model checking [30]. It is a randomized self-stabilization algorithm for leader election in a distributed ring of processors. We verify the probability that a system ofNprocessors stabilizes within a horizon ofHsteps. Figure 6 shows how each method scales on this model. Due to state explosion, all methods scale exponentially withN, but the effect is not felt by Tessa until a largerN. On the hardest instance that any baseline method can solve (N= 17), Tessa shows over 100×speedup over the next fastest method, Storm (MTBDD). The right two plots in Figure 6show that Tessa scales effectively withHas well. AtH= 500, Tessa demonstrates over 300×speedup over the next fastest method, Storm (Sparse). While the speedup inherently includes the hardware 3 Technically, the Tessa (1st) curve also includes the time taken to compile the model spec- ification into JAX. Since this compilation happens entirely in Python, we measure it after warm-up runs and add it to the 1st-run time. Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)19 s 0 s 1 s 2 s 3 s 4 s 5 s 6 p1−p q p 1−q q p 1−q 1−p 1−p p 1−p 1−p p (a) 050100 0 0.2 0.4 KL 050100 0.4 0.6 0.8 optimization step p , q (b) 0.00.20.40.60.81.0 p 0.0 0.2 0.4 0.6 0.8 1.0 q 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 KL divergence (c) Figure 7. (a) A parametric DTMC encoding a Knuth–Yao die roller. In gray (resp. white) states, a coin of biasp(resp.q) is flipped. (b) KL divergence and parameter values as gradient descent progresses. (c) Optimization landscape over parameterspandqas a contour plot. As Figure 7(c)indicates, optimal values are p= 0.5andq= 0.5, which are found by gradient descent as shown in Figure7(b). advantage of a GPU, it underscores the value of mapping the verification problem to an accelerator-friendly representation. Discussion.We caveat that Tessa outpaces these baseline methods for mod- els thatfitwithin the VRAM limit. For sparse models, Storm (Sparse) and Storm (MTBDD) are in general more space-efficient than Tessa. Nevertheless, that Tessa achieves these speedups under the11 GBVRAM constraint indicates that the method is reasonably space-efficient for the class of models it targets (thanks to the Tessa implementation and XLA optimizations exploiting model structure to reduce memory footprint), whereas baseline methods may run out of memory or time out on the same dense models. 7.2 Parameter Search We further evaluate Tessa on parameter synthesis using the Knuth–Yao algo- rithm [ 24]. This model generates a target distribution using coin flips. Figure7(a) depicts the Markov chain [20,21]. The task is to find the coin biases (pandq) that produce this target distribution. We formulate this search as an optimiza- tion problem. The objective is to minimize the Kullback–Leibler (KL) divergence [26] between the model’s output and the target. Tessa compiles the model into a differentiable tensor program, whose result can be programmatically composed with distributional properties including but not limited to state reachability—in this case, the KL divergence. The resulting objective can then be directly composed with JAX’s automatic differentiation framework and gradient-based optimizers [ 11]. Figure7(b)shows that, starting from random initialization, gradient descent successfully converges to the opti- mal parameter values well within 100 steps in a few seconds. This experiment highlights theflexibilitywith which Tessa can be used to optimize for distribu- tional properties beyond state reachability. 20Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang 8 Related Work GPU-accelerated probabilistic model checking.Bošnački et al.[7,8] study CUDA-accelerated sparse matrix–vector multiplication for DTMCs. Our approach differs fundamentally fromBošnački et al.in data representation, mem- ory access patterns, engineering simplicity, and the class of models targeted. Bošnački et al.’s approach is best for DTMCs with sparse transition dynamics, representing the transition dynamics as a flattened sparse matrix. To mitigate the overhead of indirect memory access inherent of this storage format, they build on the modified CSR format [27] and develop custom CUDA kernels for sparse matrix–vector multiplication (SpMV). Češka et al. [ 35] adapt similar sparse- matrix techniques to parameter synthesis for continuous-time Markov chains through custom GPU kernels. Also using sparse matrices, Heemstra and Wijs [ 15] perform explicit state space exploration as well as model checking entirely on the GPU. In contrast, our approach is less suitable for sparse models but effective for dense ones. It compiles to JAX’s dense tensor operations and does not require GPU programming. Rather than materializing the model into a sparse matrix (i.e., data), we map the model to tensor transformations (i.e., code), which XLA can fuse and optimize. This allows our system to make good use of the high- throughput dense linear algebra units on modern accelerators, which are often underutilized in sparse regimes. Bak et al. [ 6] use GPUs to accelerate statistical model checking (SMC) of extended timed automata. SMC is fundamentally different from probabilistic model checking (PMC): while PMC computes exact probabilities by exhaustively exploring the state space, SMC estimates probabilities via sampling, thus trading accuracy for feasibility. Achieving low error margins in SMC is at the cost of high demands on computational resources for Monte Carlo simulations. State explosion.Techniques for mitigating state explosion have been heavily studied. They either compress, abstract, or prune the state space [ 5,4,28,23, 13,14,22,12]. State-of-the-art probabilistic model checkers (e.g., Storm [16] and PRISM [29]) integrate such sophisticated state-space reduction techniques. Compared to these techniques, Tessa takes an orthogonal approach. The translation to tensor computations is largely oblivious to the state-explosion issue. Rather, we rely on a tensor compiler for fusion and optimization. In a sense, we cast state-space reduction as compiler optimizations, offloading much of the complexity of extracting high performance to a mature compiler stack. Distribution transformers.Our tensor transformer semantics is akin to the distribution transformer semantics of Kozen [25], which has found many uses in the analysis of probabilistic models (e.g., [3,34]). We recast it and establish the formal ties between the tensor transformer semantics and DTMC model checking. Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)21 9 Conclusion We cast model checking of finite-horizon Markov chains as dense tensor compu- tations, for which compiler and hardware support is readily available. This new perspective delivers sizable performance gains for models with dense transition dynamics while maintaining mathematical soundness. We hope our approach makes a useful addition to the toolbox of probabilistic model checkers, extends their practical reach, and inspires future work on tensor-based verification meth- ods. Acknowledgments.We thank the anonymous reviewers for their valuable feedback. This work was supported in part by the Natural Sciences and Engineering Research Council of Canada. The views and opinions expressed are those of the authors and do not necessarily reflect the position of any funding agency. Disclosure of Interests.The authors have no competing interests to declare that are relevant to the content of this article. Data-Availability Statement.The artifact accompanying this paper is available athttps://doi.org/10.5281/zenodo.19802567[32]. The artifact includes the Tessa imple- mentation, the benchmarks used in the evaluation, and instructions to reproduce the experimental results. 22Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang A Proofs Theorem A.1(Restatement of Theorem2.1).Pr M ( ♢ ≤n G ) = Pr M ( ι ▶ ♢ ≤n G ) . Proof.We prove the stronger claim that for any distribution (or measure)μover S, ∑ s∈S μ(s)·Pr M ( s⊨ ♢ ≤n G ) = Pr M ( μ ▶ ♢ ≤n G ) .(A.1) The theorem follows immediately by settingμ=ι. We proceed by induction on n. Base Case ( n = 0 ). From (2.1), the left-hand side of (A.1) is: ∑ s∈S μ(s)·Pr M ( s⊨ ♢ ≤0 G ) = ∑ s∈G μ(s)·1 + ∑ s /∈G μ(s)·0 = ∑ s∈G μ(s). This matches the definition ofPr M ( μ ▶ ♢ ≤0 G ) in (2.4). Inductive Case.Assume (A.1) holds forn. We considern+ 1. We split the sum on the left-hand side based on membership inGand apply (2.2): ∑ s ∈S μ(s)·Pr M ( s⊨ ♢ ≤n+1 G ) = ∑ s ∈G μ(s)·1 + ∑ s /∈G μ(s)· ( ∑ s ′ ∈S η(s, s ′ )·Pr M ( s ′ ⊨ ♢ ≤n G ) ) = ∑ s∈G μ ( s ) + ∑ s ′ ∈S ( ∑ s /∈G μ(s)η(s, s ′ ) ) ·Pr M ( s ′ ⊨ ♢ ≤n G ) . Letμ ′ = ∑ s /∈G μ(s)η(s). The term in the parentheses is exactlyμ ′ (s ′ ). By the inductive hypothesis applied toμ ′ , the second term becomes: ∑ s ′ ∈S μ ′ (s ′ )·Pr M ( s ′ ⊨ ♢ ≤n G ) = Pr M ( μ ′ ▶ ♢ ≤n G ) . Substituting this back, we get: ∑ s∈G μ(s) + Pr M ( μ ′ ▶ ♢ ≤n G ) . This matches the definition ofPr M ( μ ▶ ♢ ≤n+1 G ) in (2.5). Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)23 Lemma A.2.For any modulem k , actiona, and local states∈State(Y k ), ∑ s ′ k ∈State(X k ) (Jm k Ka s)(s ′ k ) =EnJm k Ka s. Proof.We proceed by case analysis on whether the actionabelongs to the alphabet of modulem k . Case 1:a /∈Action(m k ). By definition, Jm k Ka s=δ s| X k , EnJm k Ka s= 1. We sum the probability mass over all possible next statess ′ k : ∑ s ′ k ∈State(X k ) (Jm k Ka s)(s ′ k ) = ∑ s ′ k ∈State(X k ) δ s| X k (s ′ k ). The Dirac delta functionδis a point mass distribution. Its sum over the entire domain is exactly1: ∑ s ′ k ∈State(X k ) δ s| X k ( s ′ k ) = 1 =EnJm k Ka s. Case 2:a∈Action(m k ). By definition, ∑ s ′ k ∈State(X k ) (Jm k Ka s)(s ′ k ) = ∑ s ′ k ∈State(X k ) ∑ [a]g→U∈Com(m k ,a) [JgKs6= 0]·(JUKs)(s ′ k ) = ∑ [a]g→U∈Com(m k ,a) [JgKs6= 0] ∑ s ′ k ∈State(X k ) (JUKs)(s ′ k ) . The termJUKsis a probability distribution overState(X k ). It is defined as a con- vex combination of point masses where the probabilitiesθ i sum to1. Therefore, its total mass over the state space is exactly1: ∑ s ′ k ∈State(X k ) (JUKs)(s ′ k ) = 1. We substitute this constant back into our equation: ∑ s ′ k ∈State(X k ) (Jm k Ka s)(s ′ k ) = ∑ [a]g→U∈Com(m k ,a) [JgKs6= 0]·1 =EnJm k Ka s. 24Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang Theorem A.3(Restatement of Theorem3.1).For anys∈State, StepJMKs is a probability distribution. Proof.To prove thatStepJMKsis a valid probability distribution, we must show that the sum over all possible next statess ′ ∈Stateequals1. We consider the two cases from the definition ofStepJMKsin (3.10): Case 1: ∑ a EnJMKa s= 0. By definition,StepJMKs=δ s . The sum of proba- bilities is trivially ∑ s ′ ∈State δ s (s ′ ) = 1. Case 2: ∑ a EnJMKa s >0. We sum the probabilities over alls ′ ∈State: ∑ s ′ ∈State (StepJMKs)(s ′ ) = ∑ s ′ ∈State ∑ a ( ⊗ K k=1 Jm k Ka s| Y k ) (s ′ ) ∑ a EnJMKa s = ∑ a ∑ s ′ ∈State ( ⊗ K k=1 Jm k Ka s| Y k ) (s ′ ) ∑ a EnJMKa s Because the state spaceStateis formed by disjoint sets of variablesX k , the sum of a tensor product of distributions over the joint state space equals the product of the sums of the individual marginal distributions: ∑ s ′ ∈State ( K ⊗ k=1 Jm k Ka s| Y k ) (s ′ ) = ∑ s ′ ∈State K ∏ k=1 ( Jm k Ka s| Y k )( s ′ | X k ) = K ∏ k=1 ∑ s ′ k ∈State(X k ) ( Jm k Ka s| Y k ) (s ′ k ) = K ∏ k=1 EnJm k Ka s| Y k ( by LemmaA.2) =EnJMKa s Replacing the inner sum in the original fraction gives ∑ a EnJMKa s ∑ a EnJMKa s = 1. Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)25 Lemma A.4.The mappingΞ:D <∞ (State)→Tensor R (State)is linear. That is, for any distributionsμ 1 , μ 2 ∈D <∞ (State)and scalarsc 1 , c 2 ∈R ≥0 , Ξ(c 1 μ 1 +c 2 μ 2 ) =c 1 Ξ(μ 1 ) +c 2 Ξ(μ 2 ). Proof.For any states∈State: Ξ(c 1 μ 1 +c 2 μ 2 )[s] = (c 1 μ 1 +c 2 μ 2 )(s)( def. ofΞ) =c 1 μ 1 (s) +c 2 μ 2 (s) =c 1 Ξ(μ 1 )[s] +c 2 Ξ(μ 2 )[s](def. ofΞ) = (c 1 Ξ(μ 1 ) +c 2 Ξ(μ 2 ))[s] Since this holds for all indicess, the tensors are equal. 26Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang Lemma A.5.For any expressione, set of variablesXcontaining the free variables ofe, and states∈State(X), LeM X [s] =JeKs. Proof.By structural induction on the expressione. Lemma A.6.For any updateu, local states∈State(Y k ), and variablex∈X k , LuM Y k x[s] = (JuKs)(x). Proof.Letube the updatex 1 : =e 1 · · ·x n : =e n . We proceed by case analysis on whether the queried variablexis actively updated. Case 1:x=x i for somei. LuM Y k x i [s] =Le i M Y k [s](by (4.1)) =Je i Ks( by LemmaA.5) = ( s| X k [x i 7→Je i Ks, . . .] ) (x i )( map lookup) = (JuKs)(x i )(by (3.5)) Case 2:x6=x i for alli. LuM Y k x[s] =I x∈Y k [s]( by (4.1)) =s(x)(by definition of index tensors) = ( s| X k [x 1 7→Je 1 Ks, . . .] ) (x)(byx6=x i for alli) = (JuKs)(x)( by (3.5)) Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)27 Lemma A.7.For any mixture of updatesU, states∈State(Y k ), and next states ′ ∈State(X k ), LUM[(s, s ′ )] = (JUKs)(s ′ ). Proof.LetUbeθ 1 :u 1 +· · ·+θ n :u n . We expand the tensor evaluation at the joint index(s, s ′ ): LUM[(s, s ′ )] = n ∑ i=1 θ i · ⊙ x∈X k [ Lu i M Y k x⊗1 X k =1 Y k ⊗I x∈X k ] [(s, s ′ )]( by (4.2)) = n ∑ i=1 θ i · ∏ x∈X k [( Lu i M Y k x⊗1 X k ) [(s, s ′ )] = ( 1 Y k ⊗I x∈X k ) [(s, s ′ )] ] (point-wise evaluation) = n ∑ i=1 θ i · ∏ x∈X k [ Lu i M Y k x[s]·1 = 1·I x∈X k [s ′ ] ] (point-wise equivalence) = n ∑ i=1 θ i · ∏ x∈X k [(Ju i Ks)(x) =s ′ (x)] ( by LemmaA.6and definition of index tensors) = n ∑ i=1 θ i ·[Ju i Ks=s ′ ]( point-wise equivalence implies state equivalence) = n ∑ i=1 θ i ·δ Ju i Ks (s ′ ) = (JUKs)(s ′ )( by (3.6)) 28Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang Theorem A.8.For any modulem k in modelM, actiona, global states∈State, and global next states ′ ∈State, Lm k Ma [( s| Y k , s ′ | X k )] = ( Jm k Ka s| Y k ) (s ′ | X k ). Proof.Lets in =s| Y k ∈State(Y k )ands out =s ′ | X k ∈State(X k )be the local state projections. We proceed by case analysis on whethera∈Action(m k ). Case 1:a /∈Action(m k ). Lm k Ma[(s in , s out )] = ⊙ x∈X k [ I x∈Y k ⊗1 X k =1 Y k ⊗I x∈X k ] [(s in , s out )](by (4.3)) = ∏ x∈X k [ I x∈Y k [s in ] =I x∈X k [s out ] ] (point-wise evaluation) = ∏ x∈X k [s in (x) =s out (x)](by definition of index tensors) = [ s in | X k =s out ] ( point-wise equivalence implies state equivalence) =δ s in | X k (s out ) = (Jm k Ka s in )(s out )( by (3.7)) Case 2:a∈Action(m k ). Lm k Ma[(s in , s out )] = ∑ [a]g→U∈Com(m k ,a) ([ LgM Y k 6= 0 ] ⊗1 X k ) LUM [(s in , s out )]( by (4.3)) = ∑ [a]g→U∈Com(m k ,a) [ LgM Y k [s in ]6= 0 ] ·1·LUM[(s in , s out )] (point-wise evaluation) = ∑ [a]g→U∈Com(m k ,a) [JgKs in 6= 0]·LUM[(s in , s out )](by LemmaA.5) = ∑ [a]g→U∈Com(m k ,a) [JgKs in 6= 0]·(JUKs in )(s out )(by LemmaA.7) = (Jm k Ka s in )(s out )(by (3.7)) Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)29 Theorem A.9(Restatement of Theorem5.1).LetMbe a PML model. Letμ∈D <∞ (StateJMK). Then, LMMΞ(μ) =Ξ ( ∑ s μ(s)·StepJMKs ) . Proof.LetT=Ξ(μ), which meansT[s] =μ(s)for alls∈State. We evaluate LMMTat an arbitrary indexs ′ ∈State: (LMMT)[s ′ ] = ( ∑ a P a ) [s ′ ] + ( T [ L=0 Vars(M) ]) [s ′ ](by (4.7)) = ( ∑ a ∑ s∈State μ(s) max(L[s],1) K ∏ k=1 Lm k Ma [( s| Y k , s ′ | X k )] ) +μ(s ′ )[L[s ′ ] = 0] (by (4.8) and def. of ) = ∑ s∈State μ(s) [ ∑ a 1 max(L[s],1) K ∏ k=1 Lm k Ma [( s| Y k , s ′ | X k )] + [L[s] = 0][s=s ′ ] ] ( by linearity of summation) = ∑ s∈State μ(s) [ ∑ a 1 max(L[s],1) K ∏ k=1 Lm k Ma [( s| Y k , s ′ | X k )] + [L[s] = 0]δ s (s ′ ) ] (since[s=s ′ ]equalsδ s (s ′ )) = ∑ s∈State μ(s) [ ∑ a 1 max( L [ s ] , 1) K ∏ k=1 ( Jm k Ka s| Y k )( s ′ | X k ) + [L[s] = 0]δ s (s ′ ) ] ( by TheoremA.8) = ∑ s∈State μ(s) [ ∑ a 1 max(L[s],1) ( K ⊗ k=1 Jm k Ka s| Y k ) (s ′ ) + [L[s] = 0]δ s (s ′ ) ] ( sinceX 1 , . . . , X K are disjoint sets of variables) = ∑ s∈State μ(s)· ∑ a ( ⊗ K k=1 Jm k Ka s| Y k ) (s ′ ) L[s] ifL[s]>0 δ s (s ′ )otherwise ( by analyzing cases forL[s]) = ∑ s∈State μ(s)· ∑ a ( ⊗ K k=1 Jm k Ka s| Y k ) (s ′ ) ∑ a EnJMKa s if ∑ a EnJMKa s >0 δ s (s ′ )otherwise ( sinceL[s] = ∑ a EnJMKa s) = ∑ s∈State μ(s)·(StepJMKs)(s ′ )( by (3.10)) 30Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang =Ξ ( ∑ s∈State μ(s)·StepJMKs ) [s ′ ](by LemmaA.4) Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)31 Theorem A.10(Restatement of Theorem5.2).For anyn∈N, Pr M ( Ξ(μ)♢ ≤n e ) = Pr M ( μ ▶ ♢ ≤n G ) . Proof.We proceed by induction onn. LetT=Ξ(μ). Base Case (n= 0).From (4.9): Pr M ( T♢ ≤0 e ) =〈∆ e ,T〉= ∑ s∈State [JeKs6= 0]·T[s]. Recall thatT[s] =μ(s)andG= s ∣ ∣ ∣ JeK(s)6= 0 . Thus[JeKs6= 0] = 1ifs∈ G and0otherwise. Substituting these into the sum: ∑ s∈State [JeKs6= 0]·μ(s) = ∑ s∈G μ(s). This matches the definition ofPr M ( μ ▶ ♢ ≤0 G ) in (2.4). Inductive Case.Assume the theorem holds forn. We expand the tensor reach- ability forn+ 1using (4.10): Pr M ( T♢ ≤n+1 e ) =〈∆ e ,T〉 ︸︷︸ Term A + Pr M ( LMM(∆ ¬e T)♢ ≤n e ) ︸︷︸ Term B . Similarly, we expand the distribution reachability using ( 2.5): Pr M ( μ ▶ ♢ ≤n+1 G ) = ( ∑ s∈G μ(s) ) ︸ ︷︸ Term C + Pr M ( ∑ s /∈G μ(s)η(s) ▶ ♢ ≤n G ) ︸︷︸ Term D . From the base case logic, we know Term A=Term C. It remains to show Term B=Term D. Letμ cont ∈D <∞ (State)be the distribution of probability mass continuing from non-goal states: μ cont : = ∑ s /∈G μ(s)δ s . Consider the argument toPr M in Term D: ∑ s /∈G μ(s)η(s) = ∑ s /∈G μ(s)·StepJMKs= ∑ s μ cont (s)·StepJMKs. Now consider the tensor argument in Term B. The masked tensor∆ ¬e T has entries: (∆ ¬e T)[s] = [J¬eKs6= 0]·μ(s) = μ(s)ifs /∈ G 0otherwise 32Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang This is exactly the tensor representation ofμ cont . That is,∆ ¬e T=Ξ(μ cont ). By Theorem5.1(commutativity), we have: LMM(∆ ¬e T) =LMMΞ(μ cont ) =Ξ ( ∑ s μ cont (s)·StepJMKs ) . Thus, the tensor argument in Term B isΞ( ∑ s μ cont (s)·StepJMKs). By the induction hypothesis: Pr M ( Ξ ( ∑ s μ cont (s)·StepJMKs ) ♢ ≤n e ) = Pr M ( ∑ s μ cont (s)·StepJMKs ▶ ♢ ≤n G ) . This proves Term B=Term D, completing the induction. Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)33 References 1.PRISM manual (2025), URLhttps://prismmodelchecker.org/manual, version 4.9 2.XLA: A machine learning compiler for GPUs, CPUs, and ML accelerators. https://github.com/openxla/xla(nd), accessed: 2026-01-10 3.S. Akshay, Krishnendu Chatterjee, Tobias Meggendorfer, Ðorđe Žikelić: MDPs as distribution transformers: Affine invariant synthesis for safety objectives. In: Int’l Conf. on Computer Aided Verification (CAV) (2023), https://doi.org/10.1007/978-3-031-37709-9_5 4.Luca de Alfaro, Marta Z. Kwiatkowska, Gethin Norman, David Parker, Roberto Segala: Symbolic model checking of probabilistic processes using MTBDDs and the Kronecker representation. In: Int’l Conf. on Tools and Algorithms for the Construction and Analysis of Systems (TACAS) (2000), https://doi.org/10.1007/3-540-46419-0_27 5.Christel Baier, Edmund M. Clarke, Vasiliki Hartonas-Garmhausen, Marta Kwiatkowska, Mark Ryan: Symbolic model checking for probabilistic pro- cesses. In: Int’l Colloquium on Automata, Languages and Programming (ICALP) (1997), https://doi.org/10.1007/3-540-63165-8_199 6.Oliver S. Bak, Mathias W. B. Christiansen, Oliver V. Eriksen, Sergio Feo- Arenis, Peter G. Jensen, Marcus D. Jensen, Simas Juozapaitis, Kim G. Larsen, Marius Mikučionis, Marco Muñiz, Andreas Windfeld: GPU accel- erating statistical model checking for extended timed automata. In: Prin- ciples of Verification: Cycling the Probabilistic Landscape (2025),https: //doi.org/10.1007/978-3-031-75775-4_12 7.Dragan Bošnački, Stefan Edelkamp, Damian Sulewski, Anton Wijs: GPU- PRISM: An extension of PRISM for general purpose graphics processing units. In: Proc. of the Ninth Int’l Workshop on Parallel and Distributed Methods in Verification, and Second Int’l Workshop on High Performance Computational Systems Biology (2010), https://doi.org/10.1109/PDMC-HiBi. 2010.11 8.Dragan Bošnački, Stefan Edelkamp, Damian Sulewski, Anton Wijs: Parallel probabilistic model checking on general purpose graphics processors. Int’l Journal on Software Tools for Technology Transfer (STTT)13(2011),https: //doi.org/10.1007/s10009-010-0176-4 9.James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake Vander- Plas, Skye Wanderman-Milne, Qiao Zhang: JAX: composable transforma- tions of Python+NumPy programs (2018), URLhttp://github.com/jax-ml/ jax 10.Carlos E. Budde, Christian Dehnert, Ernst Moritz Hahn, Arnd Hartmanns, Sebastian Junges, Andrea Turrini: JANI: Quantitative model and tool interaction. In: Int’l Conf. on Tools and Algorithms for the Construc- tion and Analysis of Systems (TACAS) (2017),https://doi.org/10.1007/ 34Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang 978-3-662-54580-5_9 11.DeepMind, Igor Babuschkin, Kate Baumli, Alison Bell, Surya Bhupati- raju, Jake Bruce, Peter Buchlovsky, David Budden, Trevor Cai, Aidan Clark, Ivo Danihelka, Antoine Dedieu, Claudio Fantacci, Jonathan Godwin, Chris Jones, Ross Hemsley, Tom Hennigan, Matteo Hessel, Shaobo Hou, Steven Kapturowski, Thomas Keck, Iurii Kemaev, Michael King, Markus Kunesch, Lena Martens, Hamza Merzic, Vladimir Mikulik, Tamara Nor- man, George Papamakarios, John Quan, Roman Ring, Francisco Ruiz, Al- varo Sanchez, Laurent Sartran, Rosalia Schneider, Eren Sezener, Stephen Spencer, Srivatsan Srinivasan, Miloš Stanojević, Wojciech Stokowiec, Luyu Wang, Guangyao Zhou, Fabio Viola: The DeepMind JAX Ecosystem (2020), URLhttp://github.com/google-deepmind 12.Tom Dijk, Jaco Pol: Multi-core symbolic bisimulation minimisation. In: Int’l Conf. on Tools and Algorithms for the Construction and Analysis of Systems (TACAS) (2016),https://doi.org/10.1007/978-3-662-49674-9_19 13.Ernst Moritz Hahn, Holger Hermanns, Björn Wachter, Lijun Zhang: PASS: Abstraction refinement for infinite probabilistic models. In: Int’l Conf. on Tools and Algorithms for the Construction and Analysis of Systems (TACAS) (2010),https://doi.org/10.1007/978-3-642-12002-2_30 14.Ernst Moritz Hahn, Holger Hermanns, Lijun Zhang: Probabilistic reach- ability for parametric Markov models. Int’l Journal on Software Tools for Technology Transfer (STTT)13(1) (2011),https://doi.org/10.1007/ s10009-010-0146-x 15.Jan Heemstra, Anton Wijs: GPUexplore prob : Markov chain state space con- struction and verification with GPUs. In: Int’l Conf. on Tools and Al- gorithms for the Construction and Analysis of Systems (TACAS) (2025), https://doi.org/10.1007/978-3-031-90660-2_8 16.Christian Hensel, Sebastian Junges, Joost-Pieter Katoen, Tim Quatmann, Matthias Volk: The probabilistic model checker Storm. Int’l Journal on Soft- ware Tools for Technology Transfer (STTT) (2022),https://doi.org/10.1007/ s10009-021-00633-z 17.Ted Herman: Probabilistic self-stabilization. Information Processing Letters 35(2) (1990),https://doi.org/10.1016/0020-0190(90)90107-9 18.Steven Holtzen, Guy Van den Broeck, Todd Millstein: Scaling exact infer- ence for discrete probabilistic programs. Proc. of the ACM on Programming Languages (PACMPL)4(OOPSLA) (Nov 2020), https://doi.org/10.1145/ 3428208 19.Steven Holtzen, Sebastian Junges, Marcell Vazquez-Chanlatte, Todd Mill- stein, Sanjit A. Seshia, Guy Van den Broeck: Model checking finite- horizon Markov chains with probabilistic inference. In: Int’l Conf. on Computer Aided Verification (CAV) (2021),https://doi.org/10.1007/ 978-3-030-81688-9_27 20.Nils Jansen, Sebastian Junges, Joost-Pieter Katoen: Parameter synthesis in Markov models: A gentle survey. In: Principles of Systems Design: Essays Tensor Probabilistic Model Checking of Finite-Horizon Markov Chains (Extended Version)35 Dedicated to Thomas A. Henzinger on the Occasion of His 60th Birthday (2022),https://doi.org/10.1007/978-3-031-22337-2_20 21.Sebastian Junges, Erika Ábrahám, Christian Hensel, Nils Jansen, Joost- Pieter Katoen, Tim Quatmann, Matthias Volk: Parameter synthesis for Markov models: covering the parameter space. Formal Methods in System Design62(Feb 2024),https://doi.org/10.1007/s10703-023-00442-x 22.Nishanthan Kamaleson, David Parker, Jonathan E. Rowe: Finite-horizon bisimulation minimisation for probabilistic systems. In: Model Checking Soft- ware (SPIN’16) (2016),https://doi.org/10.1007/978-3-319-32582-8_10 23.Joost-Pieter Katoen, Tim Kemna, Ivan Zapreev, David N. Jansen: Bisimu- lation minimisation mostly speeds up probabilistic model checking. In: Int’l Conf. on Tools and Algorithms for the Construction and Analysis of Systems (TACAS) (2007),https://doi.org/10.1007/978-3-540-71209-1_9 24.Donald E. Knuth, Andrew C. Yao: The complexity of nonuniform random number generation. In: Algorithms and Complexity: New Directions and Recent Results, Academic Press (1976) 25.Dexter Kozen: Semantics of probabilistic programs. Journal of Computer and System Sciences22(3) (1981),https://doi.org/10.1016/0022-0000(81) 90036-2 26.Solomon Kullback, Richard A. Leibler: On information and sufficiency. The Annals of Mathematical Statistics22(1) (1951),https://doi.org/10.1214/ aoms/1177729694 27.Marta Kwiatkowska, Rashid Mehmood: Out-of-core solution of large lin- ear systems of equations arising from stochastic modelling. In: Proc. of the Second Joint Int’l Workshop on Process Algebra and Probabilistic Methods, Performance Modeling and Verification (2002),https://doi.org/ 10.1007/3-540-45605-8_9 28.Marta Kwiatkowska, Gethin Norman, David Parker: Symmetry reduction for probabilistic model checking. In: Int’l Conf. on Computer Aided Verification (CAV) (2006), https://doi.org/10.1007/11817963_23 29.Marta Kwiatkowska, Gethin Norman, David Parker: PRISM 4.0: Verifica- tion of probabilistic real-time systems. In: Int’l Conf. on Computer Aided Verification (CAV) (2011), https://doi.org/10.1007/978-3-642-22110-1_47 30.Marta Kwiatkowska, Gethin Norman, David Parker: Probabilistic verifica- tion of Herman’s self-stabilisation algorithm. Formal Aspects of Computing 24(Jul 2012),https://doi.org/10.1007/s00165-012-0227-6 31.Jianlin Li, Nick Guo, Peter Ye, Yizhou Zhang: Tensor probabilistic model checking of finite-horizon Markov chains. In: Int’l Conf. on Computer Aided Verification (CAV) (2026),https://doi.org/10.1007/978-3-032-32537-2_24 32.Jianlin Li, Nick Guo, Peter Ye, Yizhou Zhang: Tensor probabilistic model checking of finite-horizon Markov chains (artifact) (2026),https: //doi.org/10.5281/zenodo.19802567, also available athttps://github.com/ tessa-cav26-ae/tessa-cav26-ae 36Jianlin Li, Nick Guo, Peter Ye, and Yizhou Zhang 33.Jianlin Li, Nick Guo, Peter Ye, Yizhou Zhang: Tensor probabilistic model checking of finite-horizon Markov chains (extended version). Technical Re- port CS-2026-03, School of Computer Science, University of Waterloo (2026) 34.Jianlin Li, Yizhou Zhang: Compiling with generating functions. Proc. of the ACM on Programming Languages (PACMPL)9(ICFP) (Aug 2025),https: //doi.org/10.1145/3747534 35.Milan Češka, Petr Pilař, Nicola Paoletti, Luboš Brim, Marta Kwiatkowska: PRISM-PSY: Precise GPU-accelerated parameter synthesis for stochas- tic systems. In: Int’l Conf. on Tools and Algorithms for the Construc- tion and Analysis of Systems (TACAS) (2016),https://doi.org/10.1007/ 978-3-662-49674-9_21