On-the-Fly Computation of Bisimilarity Distances
Giorgio Bacci Giovanni Bacci Kim G. Larsen Radu Mardare


On-the-Fly Computation of Bisimilarity Distances]On-the-Fly Computation of Bisimilarity Distances* 

G. Bacci]Giorgio Bacci	
Dept. of Computer Science, Aalborg University, Denmark	
{grbacci, giovbacci, kgl, mardare}@cs.aau.dk  
Work supported by 
Sapere Aude: DFF-Young
Researchers Grant 10-085054 of the Danish Council for Independent Research, 
by the VKR Center of Excellence MT-LAB and by the Sino-Danish Basic Research Center IDEA4CPS.
	

G. Bacci]Giovanni Bacci	

	



K. G. Larsen]Kim G. Larsen	

	



R. Mardare]Radu Mardare	

	










G.3,I.1.4,I.6.4 
*An earlier version of this paper appeared as <cit.>. The present version extends <cit.> by considering the case of CTMCs and improves the linear programming approach of <cit.>.






We propose a distance between continuous-time Markov chains (CTMCs) and study the problem of computing it by comparing three different algorithmic methodologies: iterative, linear program, and on-the-fly.

In a work presented at FoSSaCS'12, Chen et al. characterized the bisimilarity distance of Desharnais et al. between discrete-time Markov chains as an optimal solution of a linear program that can be solved by using the ellipsoid method.
Inspired by their result, we propose a novel linear program characterization to compute the distance in the continuous-time setting. Differently from previous proposals, ours has a number of constraints that is bounded by a polynomial in the size of the CTMC. This, in particular, proves that the distance we propose can be computed in polynomial time.

Despite its theoretical importance, the proposed linear program characterization turns out to be inefficient in practice. Nevertheless, driven by the encouraging results of our previous work presented at TACAS'13, we propose an efficient on-the-fly algorithm, which, unlike the other mentioned solutions, computes the distances between two given states avoiding an exhaustive exploration of the state space. 
This technique works by successively refining over-approximations of the target distances using a greedy strategy, which ensures that the state space is further explored only when the current approximations are improved.

Tests performed on a consistent set of (pseudo)randomly generated CTMCs show that our algorithm improves, on average, the efficiency of the corresponding iterative and linear program methods with orders of magnitude.


[
    [
    December 30, 2023
=====================






§ INTRODUCTION
 

Continuous-time Markov chains (CTMCs) are one of the most prominent models in performance and dependability analysis. They are exploited in a broad range of applications, and constitute the underlying semantics of many modeling formalisms for real-time probabilistic systems such as Markovian queuing networks, stochastic process algebras, and calculi for systems biology.
An example of CTMC is presented in Figure <ref>(left). Here, state s_1 goes to state s_3 and s_4 with probability 1/3 and 2/3, respectively. Each state has an associated exit-rate representing the rate of an exponentially distributed random variable that characterizes the residence-time in the state
[Note that the only residence time distributions that ensure the Markov property (a.k.a., memoryless transition probability) are exponential distributions.]. For example, the probability to move from s_1 to any other state within time t ≥ 0 is given by ∫_0^t 3e^-3x dx = 1 - e^-3t. A state with no outgoing transitions (as s_3 in Figure <ref>) is called absorbing, and represents a terminating state of the system.



A key concept for reasoning about the equivalence of probabilistic systems is Larsen and Skou's  probabilistic bisimulation for discrete-time Markov chains (MCs). This notion has been extended to several types of probabilistic systems, including CTMCs. In Figure <ref>(left) s_4 and s_5 are bisimilar. Moreover, although s_1 and s_2 move with different probabilities to states s_4 and s_5, their probabilities to reach any bisimilarity class is the same, so that, also s_1 and s_2 are bisimilar. 

However, when the numerical values of probabilities are based on statistical sampling or subject to error estimates, any behavioral analysis based on a notion of equivalence is too fragile, as it only relates processes with identical behaviors. 
This issue is illustrated in Figure <ref>(right), where the states t_1 and t_2 (i.e., the counterparts of s_1 and s_2, respectively, after a perturbation of the transition probabilities) are not bisimilar. A similar situation occurs considering perturbations on the exit-rates or on associated labels, if one assumes they are taken from a metric space.

This is a common issue in applications such as, systems biology <cit.>, planning <cit.>, games <cit.>, or security <cit.>, where one is interested in knowing whether two processes that may differ by a small amount in the real-valued parameters (probabilities, rates, etc.) have “sufficiently” similar behaviors. This motivated the development of the metric theory for probabilistic systems, initiated by Desharnais et al. <cit.> and greatly developed and explored by De Alfaro, van Breugel, Worrell, and others <cit.>. It consists in proposing a bisimilarity distance (pseudometric), which measures the behavioral similarity of two models. 
These pseudometrics, e.g., the one proposed by Desharnais et al., are parametric in a discount factor that controls the significance of the future in the measurements.

Since van Breugel et al. have presented a fixed point characterization of the aforementioned pseudometric in <cit.>, several iterative algorithms have been developed in order to compute its approximation up to any degree of accuracy <cit.>.

Recently, Chen et al. <cit.> proved that, for finite MCs with rational transition function, the bisimilarity pseudometrics can be computed exactly in polynomial time. The proof consists in describing the pseudometric as the solution of a linear program that can be solved using the ellipsoid method. Although the ellipsoid method is theoretically efficient, “computational experiments with the method are very discouraging and it is in practice by no means a competitor of the, theoretically inefficient, simplex method”, as stated in <cit.>. Unfortunately, in this case the simplex method cannot be used to speed up performances in practice, since the linear program to be solved may have an exponential number of constraints in the number of states of the MC.

In this paper, we introduce a bisimilarity pseudometric over CTMCs that extends that of Desharnais et al. over MCs, and we consider the problem of computing it both from a theoretical and practical point of view.
We show that the proposed distance can be computed in polynomial time in the size of the CTMC. This is obtained by reducing the problem of computing the distance to that of finding an optimal solution of a linear program that can be solved using the ellipsoid method. Notably, differently from the proposal in <cit.>, our linear program characterization has a number of constraints that is bounded by a polynomial in the size of the CTMC. This, in particular, allows one to avoid the use of the ellipsoid algorithm in favor of the simplex or the interior-point methods.

However, also in this case, the linear program characterization turns out to be inefficient in practice, even for small CTMCs. Nevertheless, supported by the encouraging results in our previous work <cit.>, we propose to follow an on-the-fly approach for computing the distance. 
This is inspired by an alternative characterization of the bisimilarity pseudometric based on the notion of coupling structure for a CTMC. Each coupling structure is associated with a discrepancy function that represents an over-approximation of the distance. The problem of computing the pseudometric is then reduced to that of searching for an optimal coupling structure whose associated discrepancy coincides with the distance.
The exploration of the coupling structures is based on a greedy strategy that, given a coupling structure, moves to a new one by ensuring an actual improvement of the current discrepancy function. This strategy will eventually find an optimal coupling structure. The method is sound independently from the initial starting coupling structure.
Notably, the moving strategy is based on a local update of the current coupling structure. Since the update is local, when the goal is to compute the distance only between certain pairs of states, the construction of the coupling structures can be done on-the-fly, delimiting the exploration only on those states that are demanded during the computation.

The efficiency of our algorithm has been evaluated on a significant set of randomly generated CTMCs. The results show that our algorithm performs orders of magnitude better than the corresponding iterative and linear program implementations. Moreover, we provide empirical evidence that our algorithm enjoys good execution running times.

One of the main practical advantages of our approach consists in that one can focus on computing only the distances between states that are of particular interest. This is useful in practice, for instance when large systems are considered and visiting the entire state space is computationally expensive. A similar issue has been considered by Comanici et al., in <cit.> in the case of Markov decision processes with rewards, who noticed that for computing the approximated pseudometric one does not need to update the current value for all the pairs at each iteration, but it is sufficient only to focus on the pairs where changes are happening rapidly. In our approach, the termination condition is checked locally, still ensuring that the local optimum corresponds to the global one. 
Our methods can also be used in combination with approximation techniques as, for instance, to provide a least over-approximation of the behavioral distance given over-estimates of some particular distances.



  
Synopsis: 
The paper is organized as follows. In Section <ref>, we recall the basic preliminaries on continuous-time Markov chains and define the bisimilarity pseudometric. Section <ref> is devoted to the analysis of the complexity of the problem of computing such a distance. Here, two approaches are considered: an approximate iterative method and a linear program characterization. In Section <ref>, we provide an alternative characterization of the distance based on the notion of coupling structure. This is the basis for the development of an on-the-fly algorithm (Section <ref>) for the computation of the pseudometric, whose correctness and termination is proven in Section <ref>. The efficiency of this algorithm is supported by experimental results, shown in Section <ref>. Final remarks and conclusions are in Section <ref>.



§ CONTINUOUS-TIME MARKOV CHAINS AND BISIMILARITY PSEUDOMETRICS
 

We recall the definitions of (finite) L-labelled continuous-time Markov chains (CTMCs) for a nonempty set of labels L, and stochastic bisimilarity over them. Then, we introduce a behavioral pseudometric over CTMCs to be considered as a quantitative generalization of stochastic bisimilarity.

Given a finite set X, a discrete probability distribution over it is a function μ X → [0,1] such that μ(X) =1, where μ(E) = ∑_x ∈ Eμ(x), for E ⊆ X. We denote the set of finitely supported discrete probability distributions over X by (X).

[Continuous-time Markov chain]

An L-labelled continuous-time Markov chain is a tuple = (S, A, τ, ρ, ℓ) consisting of a nonempty finite set S of states, a set A ⊆ S of absorbing states, 
a transition probability function τ S ∖ A →(S), an exit rate function ρ S ∖ A →[>0], and a labeling function ℓ S → L.


The labels in L represent properties of interest that hold in a particular state according to the labeling function ℓ S → L.  If s ∈ S is the current state of the system and E ⊆ S is a subset of states, τ(s)(E) ∈ [0,1] corresponds to the probability that a transition from s to arbitrary s' ∈ E is taken, and ρ(s) ∈[>0] represents the rate of an exponentially distributed random variable that characterizes the residence time in the state s before any transition is taken. Therefore, the probability to make a transition from state s to any s' ∈ E within time unit t ∈ is given by τ(s)(E) ·[ρ(s)]([0,t)), where [r](B) = ∫_B r e^-rx dx, for any Borel subset B ⊆ and r > 0.
Absorbing states in A ⊆ S are used to represent termination or deadlock states. An example of a CTMC is shown in Figure <ref>.

For discrete-time Markov chains, the standard notion of behavioral equivalence is probabilistic bisimulation of Larsen and Skou <cit.>. The following definition extends it to CTMCs. To ease the notation, for E ⊆ S, we introduce the relation ≡_E⊆ S × S defined by s ≡_E s' if either s,s' ∈ E or s,s' ∉ E.

[Stochastic Bisimulation] 
Let = (S, A, τ, ρ, ℓ) be a CTMC. An equivalence relation R ⊆ S × S is a stochastic bisimulation on  if whenever s R t, then

  
  * s ≡_A t, ℓ(s) = ℓ(t), and
  
  * if s,t ∉A, then ρ(s) = ρ(t) and, for all C ∈ S/_R, τ(s)(C) = τ(t)(C).

Two states s, t ∈ S are bisimilar with respect to , written s ∼_ t, if they are related by some probabilistic bisimulation on .

Intuitively, two states are bisimilar if they have the same labels, they agree on being absorbing or not, and, in the case they are non-absorbing, their residence-time distributions and probability of moving 
by a single transition to any given class of bisimilar states is always the same. 
As an example of two stochastic bisimilar states, consider s_1 and s_2 in the CTMC depicted on the left hand side of Figure <ref>. A bisimulation relation that relates them is the equivalence relation with equivalence classes given by s_1,s_2, s_3, and s_4, s_5.





 §.§ Bisimilarity Pseudometrics on CTMCs
 
In this section, we introduce a family of pseudometrics on CTMCs parametric in a discount factor λ∈ (0,1). Following the approach of <cit.>, given a CTMC = (S, A, τ, ρ, ℓ) we define a (1-bounded) pseudometric on S as the least fixed point of an operator on the set  of functions from S × S to [0,1]. This pseudometric is then shown to be adequate with respect to stochastic bisimilarity: we prove that two states are stochastic bisimilar if and only if they have distance zero.

Recall that d  X × X → is a pseudometric on a set X if d(x,x) = 0 (reflexivity), d(x,y) = d(y,x) (symmetry)  and d(x,y) + d(y,z) ≥ d(x,z) (triangular inequality), for arbitrary x,y,z ∈ X; it is a metric if, in addition, d(x,y) = 0 iff x = y. A pair (X,d) where d is a (pseudo)metric on X is called a (pseudo)metric space.

Hereafter we will assume L to be equipped with a 1-bounded metric
[Since the set S of states is assumed to be finite, one may assume the set of labels to be so as well. Thus, the metric d_L on labels can be bounded without loss of generality.] d_L  L × L → [0,1].

The operator we are going to introduce will use three key ingredients: the distance d_L between labels, a distance between residence-time distributions, and a distance between transition distributions. The first is meant to measure the static differences with respect to the labels associated with the states, the last two are meant to capture the differences in the dynamics, respectively, with respect to the continuous and discrete probabilistic choices. 

To this end, we consider two distances over probability distributions. The first one is the total variation metric, defined for arbitrary Borel probability measures μ, ν over  as  

    μν = sup_E |μ(E) - ν(E)|  ,

where the supremum is taken over the Borel measurable sets of .
The second one is the Kantorovich distance, which is based on the notion of coupling of probability measures, which we introduce next in the case of probability distributions over finite sets.
[Coupling]
Let S be a finite set, and let μ, ν∈(S). A probability distribution ω∈(S × S) is
said a coupling for (μ, ν) if, for arbitrary u,v ∈ S 

    ∑_v ∈ Sω(u, v) = μ(u) 
            and      ∑_u ∈ Sω(u, v) = ν(v)  .

In other words, ω is a joint probability distribution with left and right marginal, respectively, given by μ and ν. We denote the set of couplings for (μ, ν) by μν.


For a finite set S and a 1-bounded distance d  S × S → [0,1] over it, the Kantorovich distance is defined, for arbitrary distributions μ, ν∈(S) as follows
[The minimum can be used in place of an infimum thanks to the Fenchel-Rockafeller duality theorem (see <cit.>).]

    𝒦_d(μ, ν) = min∑_u,v ∈ S d(u,v) ·ω(u,v) ω∈μν .

Intuitively, 𝒦_d lifts a (1-bounded) distance over S to a (1-bounded) distance over its probability distributions. One can show that 𝒦_d is a (pseudo)metric if d is a (pseudo)metric.

Now, consider the following functional operator.
 
Let = (S, A, τ, ρ, ℓ) be CTMC and λ∈ (0,1) a discount factor.
The function Δ^_λ→ is defined as follows, for d  S × S → [0,1] and s, t ∈ S

    Δ^_λ(d)(s,t) =
      
        1    if s ≢_A t
    st   if s,t ∈ A
    maxst, λ·st   if s,t ∉ A

where []→ and ,  S × S → [0,1] are respectively defined by

    st = st  + (1 - st) ·𝒦_d(τ(s),τ(t))  ,  
    st = d_L(ℓ(s),ℓ(t))  ,
          and    st = [ρ(s)][ρ(t)] .


The functional Δ^_λ measures the difference of two states with respect to: their labels (by means of the pseudometric ), their residence-time distributions (by means of the pseudometric ), and their discrete probabilities to move to the next state (by means of the Kantorovich distance). If two states disagree on being absorbing (or not) they are considered incomparable, and their distance is set to 1. If both states are absorbing, they express no dynamic behavior, hence they are compared statically, and their distance corresponds to that occurring between their respective labels. Finally, if the states are non-absorbing, then they are compared with respect to both their static and dynamic features, namely, taking the maximum among their respective associated distances.

Specifically, the value st corresponds to the least probability that two transitions are taken independently from the states s and t at different moments in time. 
This value is used by the functional [] to measure the overall differences that might occur in the dynamics of the two states in combination with the Kantorovich distance between their transition probability distributions (see Remark <ref> for more details).

The set  is endowed with the partial order ⊑ defined as d ⊑ d' iff d(s,t) ≤ d'(s,t) for all s,t ∈ S and it forms a complete lattice. The bottom element 0 is the constant 0 function, while the top element is the constant 1 function. For any subset D ⊆, the least upper bound D, and greatest lower bound D are, respectively, given by ( D)(s,t) = sup_d ∈ D d(s,t) and ( D)(s,t) = inf_d ∈ D d(s,t), for all s,t ∈ S.



It is easy to check that, for any  and λ∈ (0,1), Δ^_λ is monotone (i.e., whenever d ⊑ d', then Δ^_λ(d) ⊑Δ^_λ(d')), thus, since (, ⊑) is a complete lattice, by Tarski's fixed point theorem Δ^_λ admits least and greatest fixed points.

[Bisimilarity distance] 
Let  be a CTMC and λ∈ (0,1).
The λ-disc­oun­ted bisimilarity pseudometric on ℳ,  
denoted by ^, is the least fixed point of Δ^ℳ_λ.


The rest of the section is devoted to show that the least fixed point ^ is indeed a pseudometric and, moreover, is adequate with respect to stochastic bisimilarity (Theorem <ref>). This justifies the definition above.
To this end we need some technical lemmas. In particular, we prove that Δ^_λ preserves pseudometrics (Lemma <ref>) and it is Lipschitz continuous (Lemma <ref>). 

Hereafter, unless mentioned otherwise, we fix a CTMC = (S, A, τ, ρ, ℓ) and a discount factor λ∈ (0,1). To ease the notation Δ^_λ, ^, and ∼_ will be denoted simply by Δ_λ, , and ∼, respectively, whenever  is clear from the context.

 
The operator Δ_λ preserves pseudometrics.


Let d  S × S → [0,1] be a pseudometric. We want to prove that Δ_λ(d) is a pseudometric. Recall that ,  S × S → [0,1] are pseudometrics. Thus, since the point wise maximum of pseudometrics is a pseudometric, it suffices to prove that [] preserves pseudometrics. Recall that, 𝒦_d (S) ×(S) → [0,1] is a pseudometric, since d is so. Thus, reflexivity and symmetry are immediate. The only nontrivial case is the triangular inequality. Let s,t,u ∈ S ∖ A, we want to prove st≤su + ut. First note that, for 0 ≤β≤ 1 and α' ≥α, the following holds:

    α + (1- α) β   
      
      
      = β + (1 - β) αdistributivity
       ≤β + (1 - β) α' α≤α' and 0 ≤β≤ 1
       = α' + (1 - α') β . distributivity

Thus, since  is a pseudometric, by triangular inequality and the above we have 

    st   = st  + (1 - st) ·𝒦_d(τ(s),τ(t))  def. []
       ≤su + ut + (1 - (su + ut) ) ·𝒦_d(τ(s),τ(t))
      *
    If we show that the last summand in (<ref>) is less than or equal to the sum of (1 - su) ·𝒦_d(τ(s),τ(u)) and (1 - ut) ·𝒦_d(τ(u),τ(t)),
    we get the following, and we are done
       ≤su  + (1 - su) ·𝒦_d(τ(s),τ(u)) +
       ut  + (1 - ut) ·𝒦_d(τ(u),τ(t)) 
       = su +   utdef. [] .

To this end, consider two cases. If su + ut > 1 then the inequality holds trivially, since 1 - (su + ut) < 0, so that the last summand in (<ref>) is negative and 1-su≥ 0. 
Instead, if su + ut≤ 1 then 1 - (su + ut) ≥ 0, so we have

    (1 - (su + ut) ) ·𝒦_d(τ(s),τ(t)) 
         ≤(1 - (su + ut) ) ·(𝒦_d(τ(s),τ(u)) + 𝒦_d(τ(u),τ(t)) ) 
      triang. 𝒦_d
         = 
      (1 - (su + ut) ) ·𝒦_d(τ(s),τ(u)) +
      (1 - (su + ut) ) ·𝒦_d(τ(u),τ(t)) 
         ≤(1 - su) ·𝒦_d(τ(s),τ(u)) +
      (1 - ut) ·𝒦_d(τ(u),τ(t))

and we are done.


The set  can be turned into a metric space by means of the supremum norm dd' = sup_s,t ∈ S |d(s,t) - d'(s,t)|. Next we show that the λ-discounted functional operator Δ_λ is λ-Lipschitz continuous, that is Δ_λ(d')Δ_λ(d)≤λ·d'd, for any d,d' ∈.
 
The operator Δ_λ is λ-Lipschitz continuous.

By <cit.>, to prove that Δ_λ is λ-Lipschitz continuous it suffices to show that, whenever d ⊑ d' then, for all s,t ∈ S, Δ_λ(d')(s, t) - Δ_λ(d)(s, t) ≤λ·d'd.
If s ≢_A t, then Δ_λ(d')(s, t) - Δ_λ(d)(s, t) = 1 - 1 = 0, so that, the inequality is satisfied. If s,t ∉ A, Δ_λ(d')(s, t) - Δ_λ(d)(s, t) = st - st = 0, and again the inequality holds trivially. The same happens when st≥λ·[d']st. Indeed, by monotonicity of [], st≥λ·st holds, so that, by definition of Δ_λ we have

    Δ_λ(d')(s,t) - Δ_λ(d)(s,t)
      = λ·(  st - st) = 0 ≤λ·d'd .

It remains the case when s,t ∉ A and λ·st≤st < λ·[d']st. 
Assume that 𝒦_d(τ(s),τ(t)) = ∑_u,v ∈ S d(u,v) ·ω(u,v), for some ω∈τ(s)τ(t), then we have

    Δ_λ(d')(s,t) - Δ_λ(d)(s,t)
         ≤λ·( [d']st - [d]st) 
       ≤λ·( 𝒦_d'(τ(s),τ(t)) - 𝒦_d(τ(s),τ(t)) ) 
       ≤λ·(
      ∑_u,v ∈ S d'(u,v) ·ω(u,v) - ∑_u,v ∈ S d(u,v) ·ω(u,v) ) 
       = λ·(
      ∑_u,v ∈ S (d'(u,v) - d(u,v)) ·ω(u,v) ) 
       ≤λ·(
      ∑_u,v ∈ Sd'd·ω(u,v) ) 
       = λ·(
       d'd·∑_u,v ∈ Sω(u,v) ) 
       = λ·d'd . to 221 pt


It is standard that  with the supremum norm forms a complete metric space (i.e, every Cauchy sequence converges). Therefore, since λ∈ (0,1), a direct consequence of Lemma <ref> and Banach's fixed point theorem is the following.
 
For any λ∈ (0,1),  is the unique fixed point of Δ_λ. Moreover, for any n ∈, and d  S × S → [0,1] we have Δ^n_λ(d)≤λ^n/1-λΔ_λ(d)d.


Now we are ready to state the main theorem of this section.
 
 is a pseudometric. Moreover, for any s,t ∈ S, s ∼ t if and only if (s,t) = 0


We first prove that  is a pseudometric. By Lemma <ref> and Banach's fixed point theorem, = _n ∈Δ^n_λ(0). Clearly, 0 is a pseudometric.  Thus by Lemma <ref>, a simple induction on n shows that, for all n ∈, Δ^n_λ(0) is a pseudometric. Since the least upper bound with respect to ⊑ preserves pseudometrics, we have that  is so.







Now we are left to prove that, for any s,t ∈ S, s ∼ t iff (s,t) = 0.
 0ex

  * (⇐) We prove that R =  (s,t) (s,t) = 0 is a stochastic bisimulation. Clearly, R is an equivalence. Assume (s,t) ∈ R, then, by definition of Δ_λ, one of the following holds:

  
  *  s,t ∈ A and st = 0; 
  
  *  s,t ∉ A, st = 0, and []st = 0.

If (<ref>) holds, by st = 0 we get that ℓ(s) = ℓ(t).
If (<ref>) holds, we have st = 0 and 𝒦_(τ(s),τ(t)) = 0. By st = 0 we get [ρ(s)] = [ρ(t)] and hence ρ(s) = ρ(t). By <cit.>, 𝒦_(τ(s),τ(t)) = 0 implies that, for all C ∈ S/_R, τ(s)(C) = τ(t)(C). Therefore R is a bisimulation.


  * (⇒) 
Let R ⊆ S × S be a stochastic bisimulation on , and define d_R  S × S → [0,1] by d_R(s,t) = 0 if (s,t) ∈ R and d_R(s,t) = 1 if (s,t) ∉ R. We show that Δ_λ(d_R) ⊑ d_R. If (s,t) ∉ R, then d_R(s,t) = 1 ≥Δ_λ(d_R)(s,t). If 
(s,t) ∈ R, then ℓ(s) = ℓ(t) and one of the following holds:

	
  *  s,t ∈ A;
	
  *  s,t ∉ A, ρ(s) = ρ(t) and,
	∀ C ∈ S/_R.  τ(s)(C) = τ(t)(C).

If (<ref>) holds, Δ_λ(d_R)(s,t) = st = 0 = d_R(s,t). 
If (<ref>) holds, by <cit.> and the fact that, for all C ∈ S/_R, τ(s)(C) = τ(t)(C), we have 𝒦_d_R(τ(s),τ(t)) = 0. Moreover st = 0. This gives  that Δ_λ(d_R)(s,t) = 0 = d_R(s,t).

Since ∼ is a stochastic bisimulation, Δ_λ(d_∼) ⊑ d_∼, so that, by Tarski's fixed point theorem, ⊑ d_∼. By definition of d_∼ and  ⊑ d_∼, s ∼ t implies (s,t) = 0.












































 
Observe that Theorem <ref> holds for alternative definitions of the functional Δ_λ. 
An example is given when the functional [] in the definition of Δ_λ is given as 

    [d]st = max1/M· |ρ(s) - ρ(t)|, 𝒦_d(τ(s),τ(t)) ,

where M = max_s ∈ Sρ(s) is used to rescale the symmetric difference |ρ(s) - ρ(t)| to a value within [0,1].
Another example is obtained by replacing the maximum above by a convex combination of the two values as below, 
for some α∈ (0,1),

    [d]st = α·1/M |ρ(s) - ρ(t)| + (1- α) ·𝒦_d(τ(s),τ(t))  ,


So, what does it make a proposal for Δ_λ preferable to another? Although Theorem <ref> is an important property for a behavioral pseudometric, it does not say much about the states that have distance different from zero. In this sense, a good behavioral metric should relate the distance with a concrete problem. Our definition of Δ_λ, for example, is motivated by a result in <cit.> that states that the total variation distance of CTMCs (more generally, on semi-Markov chains) is logically characterized as the maximal difference w.r.t. the likelihood for two states to satisfy the same Metric Temporal Logic (MTL) formula <cit.>. It turns out that when d_L is the discrete metric over L (i.e., d_L(a,b) = 0 if a = b, and 1 otherwise),  bounds from above the total variation distance. This relates  to the probabilistic model checking problem of MTL-formulas against CTMCs. The above alternative proposals for a distance do not enjoy this property.




§ COMPLEXITY AND LINEAR PROGRAMMING REPRESENTATION
 





In this section, we study the problem of computing the bisimilarity distance by considering two different approaches. The former is an iterative method that approximates  from below (resp. above) successively applying the operator Δ_λ starting from the least (resp. greatest) element in . The latter is based on a linear program characterization of  that is based on the Kantorovich duality <cit.>. In contrast to an analogous proposal in <cit.>, our linear program has a number of constraints that is polynomially bounded in the size of the CTMC. As a consequence, the bisimilarity distance  can be computed in polynomial time in the size of the CTMC.



 §.§ Iterative method
 
By Theorem <ref>, for any ϵ > 0, it follows that to get ϵ-close to , it suffices to iterate the application of the fixed point operator ⌈log_λϵ⌉ times.
 
For any ϵ > 0 and d  S × S → [0,1], Δ^⌈log_λϵ⌉_λ(d)≤ϵ.


By Theorem <ref> we have Δ^n_λ(d)≤λ^n/1-λΔ_λ(d)d and by
Δ_λ(d)d≤ 1, we have Δ^n_λ(d)≤λ^n/1-λ. For n = log_λ (ϵ - ϵλ), we have ϵ = λ^n/1-λ. Therefore, by Lemma <ref> and ⌈log_λϵ⌉≥log_λ (ϵ - ϵλ), we have Δ^⌈log_λϵ⌉_λ(d)≤ϵ.


By the above result, we obtain a simple method for approximating . If the starting point is 0 we obtain an under-approximation, whereas starting from 1 we get an over-approximation. Both the approximations can be taken arbitrary close to the exact value.



However, as shown in the following example, the exact distance value cannot be reached in general. This holds for any discount factor.
[<cit.>]
Consider the {red, blue}-labeled CTMC below. 


[
 label1/.style=circle, draw=red!50,fill=red!20, thin, circle split, inner sep=0.6mm,
 label2/.style=circle,draw=blue!50,fill=blue!20, thin, circle split, inner sep=0.6mm]
2.5cm
(0,0) node[label1] (s) s lower 1
  ((s) +(right:)) node[label1] (t) t lower 1
  ((t) +(right:)) node[label2] (u) u lower 1;

[-latex, font=]
  (s) edge[loop above] node[above] 1 (s)
  (t) edge[loop above] node[above] λ (t)
       edge node[above] 1-λ (u)
  (u) edge[loop above] node[above] 1 (u);



Let d_L  L × L → [0,1] be the discrete metric over L, defined as d_L(l,l') = 0 if l = l' and 1 otherwise.
One can check that (s,t) = λ - λ^2/1- λ^2 and, for all n ∈, Δ^n_λ(0)(s,t) ≤λ - λ^2n+1/1+λ. Since, for all n ∈, λ - λ^2n+1/1+λ < λ - λ^2/1- λ^2, we have that the fixed point cannot be reached in a finite number of iterations. 


In <cit.> is shown that the bisimilarity distance of Desharnais et al. <cit.> can be computed exactly by iterating the fixed point operator up to a precision that allows one to use the continued fraction algorithm to yield the exact value of the fixed point. This method can be applied provided that the pseudometric has rational values. In their case, this is ensured assuming that the transition probabilities are rational. Unfortunately, in our case this cannot be ensured under the same conditions. Indeed, the total variation distance between exponential distributions with rates r, r' > 0 is analytically 
solved as follows

    [r][r'] = 
    
    0    if r = r'
    | ( r'/r)^r/r - r' - 
    ( r'/r)^r'/r - r'|    otherwise

thus, even restricting to rational exit-rates and probabilities, the distance may assume irrational values. 
As a consequence, we cannot assume to compute, in general, the exact distance values.



 §.§ Linear Program Characterization
 
Our linear program characterization leverages on two key results. The first one is the uniqueness of the fixed point of Δ_λ (Theorem <ref>). The second one is a dual linear program characterization of the Kantorovich distance.

For S finite, d  S × S → [0,1], and μ, ν∈(S), the value 𝒦_d(μ,ν) coincides with the optimal value of the following linear program

    𝒦_d(μ,ν)    = 
    [t]
    min_ω    ∑_u,v ∈ S d(u,v) ·ω_u,v
       ∑_vω_u,v = μ(u)       ∀ u ∈ S 
       ∑_uω_u,v = ν(v)       ∀ v ∈ S  
       ω_u,v≥ 0       ∀ u,v ∈ S   .
    By a standard argument in linear optimization, the above can be alternatively represented by the following dual linear program𝒦_d(μ,ν)    = 
    [t]
    max_y    ∑_u ∈ S (μ(u) - ν(u) )· y_u 
        y_u - y_v≤ d(u,v)       ∀ u,v ∈ S   .

This alternative characterization is a special case of a more general result commonly known as the Kantorovich duality and extensively studied in linear optimization theory (see <cit.>).

Let n = |S| and h = n - |A|. Consider the linear program in Figure <ref>, hereafter denoted by D, with variables d ∈[]^n^2, y ∈[]^h^2 + n and  k,m ∈[]^h^2. The constraints of D are easily seen to be bounded and feasible. Moreover, the objective function of D attains its optimal value when the vectors k and m are maximized in each component[Indeed, for arbitrary x,y ∈^n and n ∈, if x_i ≤ y_i for all i = 1..n, then ∑_i =1^n x_i ≤∑_i =1^n y_i.]. Therefore, according to (<ref>), an optimal solution (d^*,y^*,k^*,m^*) ∈ D satisfies the following equalities

    ∀ s,t ∉A.       m^*_s,t   = min{st, λ( st + (1 - st) 𝒦_d^*(τ(s), τ(t)) ) } , 
       ∀ s,t ∉A.      k^*_s,t   = 𝒦_d^*(τ(s), τ(t))  .

By the above equalities and the constraints of D, it follows that d^* is a fixed point of Δ_λ. Since the distance is the unique fixed point of Δ_λ, d^* =.


 
Let (d^*,y^*,k^*,m^*) be a solution of D. Then, for all s,t ∈ S, d^*_s,t = (s,t).


Consider d^*∈[]^|S|^2 and the linear program above, hereafter denoted by D(d^*)

    *arg max_y,k,m    ∑_s,t ∉A k_s,t + m_s,t
        m_s,t≤st      ∀ s,t ∉A 
        m_s,t≤λ( st + (1 - st) k_s,t)        ∀ s,t ∉A 
        k_s,t = ∑_u ∈ S (τ(s)(u) - τ(t)(u))· y_u^s,t      ∀ s,t ∉A
        y_u^s,t - y_v^s,t≤ d^*_u,v      ∀ s,t ∉A, ∀ u,v ∈ S


Any feasible solution of D(d^*) can be improved by increasing any value of k_s,t or m_s,t for some s,t ∉ A. On the one hand, the value of k_s,t (for some s,t ∉ A) can be increased independently from that of (m_s,t)_s,t ∉ A since it only depends on the values of (y^s,t_u)_u ∈ S. For this reason, if we denote by k^*_s,t an optimal value for k_s,t we have the following equality

    k^*_s,t = max_y   ∑_u ∈ S (τ(s)(u) - τ(t)(u) )· y_u
        y_u - y_v≤ d^*_u,v      ∀ u,v ∈ S

Thus, by Equation (<ref>), we have that k^*_s,t = 𝒦_d^*(τ(s), τ(t)) for all s,t ∉A. 

On the other hand, the value of m_s,t (for some s,t ∉ A) is bounded from above by the constant st and the λ( st + (1 - st) k_s,t). The value of λ( st + (1 - st) k_s,t) increases with that of k_s,t, hence, if if we denote by m^*_s,t the optimal value for m_s,t, we have that

    ∀ s,t ∉A.       m^*_s,t   = min{st, λ( st + (1 - st) k^*_s,t) }
       ∀ s,t ∉A.      k^*_s,t   = 𝒦_d^*(τ(s), τ(t))   .

Let (d',y',k',m') be an optimal solution of D and (y”,k”,m”) be an optimal solution of D(d'). 
Notice that the optimal value ∑_s,t ∉A k”_s,t + m”_s,t of D(d') (i.e., ∑_s,t ∉A k”_s,t + m”_s,t) is greater than or equal to the optimal value of D (i.e., ∑_s,t ∉A k'_s,t + m'_s,t), since the constraints of D(d') are a subset of those of D.
Consider now two possible cases. 

(Case 1:  k' = k” and m' = m”) For s,t ∉ A we have that the following equalities hold

    d'_s,t   = st + λ( st + (1 - st)   k'_s,t) - m'_s,tby (d',y',k',m') feasible for D
       = st + λ( st + (1 - st)   k”_s,t) - m”_s,tby k' = k” and m' = m”
       = max{st ,  λ( st + (1 - st)   k”_s,t) }by (<ref>)
       = max{st, λ( st + (1 - st)  𝒦_d'(τ(s),τ(t)) ) }by (<ref>)
       = max{st, λ𝒯(d')(s,t) }by def. 𝒯
       = Δ(d')(s,t) by def. Δ_λ

From the above equality and feasibility of (d',y',k',m') for D, we have d'_s,t = Δ_λ(d')(s,t) for all s,t ∈ S. Thus, by Theorem <ref> we obtain d'_s,t = (s,t) for all s,t ∈ S. 

(Case 2: k' ≠ k” or m' ≠ m”) By the previous case we have that the optimal value of D is equal to ∑_s,t ∉A k”_s,t + m”_s,t. This is a sum of non negative values, indeed by (s,t) ≥ 0 for all s,t ∈ S and Equations (<ref>) and (<ref>), we have that k”_s,t≥ 0 and m”_s,t≥ 0 for all s,t ∉ S. Thus, if k' ≠ k” or m' ≠ m” we have ∑_s,t ∉A k”_s,t + m”_s,t > ∑_s,t ∉A k'_s,t + m'_s,t. This contradicts the assumption that (d',y',k',m') is optimal for D.

This proves that if (d',y',k',m') is an optimal solution of D then d' =.













































In <cit.> it has been shown that the bisimilarity distance of Desharnais et al. on discrete-time Markov chains can be computed in polynomial time as the solution of a linear program that can be solved by using the ellipsoid method.
However, in their proposal the number of constraints may be exponential in the size of the model. This is due to the fact that the Kantorovich distance is resolved listing all the couplings that correspond to vertices of the transportation polytopes involved in the definition of the distance <cit.>.

 In contrast, our proposal has a number of constraints and unknowns[Actually, the variables k are used only for ease the presentation but they can be removed by substitution.] bounded by 3 |S|^2 + |S|^4 and |S|^3 + 2 |S|^2, respectively.
This allows one to use general algorithms for solving LP problems (such as the simplex and the interior-point methods) that, in practice, are more efficient than the ellipsoid method.


Our LP characterization (Theorem <ref>) exploits the fact that, for λ∈ (0,1), the functional operator Δ_λ has unique fixed point (Theorem <ref>). In the so called undiscounted case, i.e., when λ = 1, the fixed point is not unique anymore and the LP characterization does not work directly. To obtain a similar LP characterization with a number of constraints polynomial in the size of the CTMC, by following <cit.> one may modify the operator Δ_1 to an operator Δ' that forces the distance to be 0 if the states are stochastic bisimilar, otherwise it behaves as Δ_1. The operator Δ' has unique fixed point when d_L is the discrete metric over L (i.e., d_L(a,b) = 0 if a = b, and 1 otherwise), but for other choices of d_L the uniqueness may not hold.


Moreover, this allows us to state the following complexity result.
 
 can be computed in polynomial-time in the size of .
 

By Theorem <ref>,  can be computed within the time it takes to construct and solve D. D has a number of constraints and unknowns that is bounded by a polynomial in the size of , therefore its construction can be performed in polynomial time. For the same reason, D admits a polynomial time separation algorithm: whenever a solution is given, its feasibility is checked by scanning each inequality in D; otherwise the first encountered inequality that is not satisfied is returned as a separation hyperplane. Therefore, the thesis follows by solving D using the ellipsoid method together with the naïve separation algorithm described above.



In Theorem <ref> the term “computed” should be replaced by “approximated”. We already discussed about the impossibility of obtaining an exact value of  (the total variations distance between exponential distributions may assume irrational values!). 

In fact, for the construction of the linear program D one has to use an approximated (rational) version of the coefficients st, for all s,t ∉ A, say e_s,t. To ensure that the optimal solution d^* of D is at most ε apart from  (i.e., |(s,t) - d^*_s,t | ≤ε, for all s,t ∈ S) one only needs that |st - e_s,t | ≤ε/2, for all s,t ∉ A. 













This can be done using the Newton-Raphson's iteration algorithm for approximating the n-square root of a number ζ∈ up to a precision ϵ > 0. This method is known to be polynomially computable in the size of the representation of ζ, n and ϵ <cit.>. Hence  can be approximated up to any precision ε > 0 in polynomial-time in the size of . 




§ ALTERNATIVE CHARACTERIZATION OF THE PSEUDOMETRIC
 

In the following, we propose an alternative characterization of the bisimilarity distance , based on the notion of coupling structure. Our result generalizes the one proposed in <cit.> for MCs, to the continuous-time settings.

[Coupling Structure] 
Let ℳ = (S, A, π, ℓ) be a CTMC. A coupling structure for  is a function (S∖ A) × (S∖ A) →(S × S) such that, for all s,t ∉ A, (s,t) ∈τ(s)τ(t).

Intuitively, a coupling structure for  can be thought of as an S × S-indexed collection of joint probability distributions, each having left/right marginals equal to τ.

The following definition adapts the definition of the operator Δ_λ (see Definition <ref>) with respect to the notion of coupling structure for a CTMC.

Let = (S, A, τ, ρ, ℓ) be CTMC,  a coupling structure for , and λ∈ (0,1) a discount factor. The function Γ^_λ→ is defined as follows, for d  S × S → [0,1] and s, t ∈ S

    Γ^_λ(d)(s,t) =
      
        1    if s ≢_A t
    st   if s,t ∈ A
    maxst, λ·st   if s,t ∉ A

where ,  S × S → [0,1] are as in Definition <ref> and []→ is given by

    st = st  + (1 - st) ·∑_u,v ∈ S d(u,v) ·(s,t)(u,v)   .


Recall that the Kantorovich distance between two distributions μ and ν is defined as 𝒦_d(μ,ν) = min_ω∑_u,v ∈ S d(u,v) ·ω(u,v), where the minimum is taken over all the possible couplings ω∈μν. Thus, the operator Γ^_λ can intuitively be thought of as a possible instance of Δ^_λ with respect to a fixed choice of the couplings given by . 

One can easily check that Γ^_λ is monotone, thus, by Tarski's fixed point theorem, it admits least and greatest fixed points. The least fixed point, in particular, will be denoted by C and referred to as the λ-discrepancy of 𝒞.























 
= minC𝒞 coupling structure for ℳ.


We first prove that δ_λ⊑C, for any coupling structure  for . By Tarski's fixed point theorem, it suffices to prove that, for any d  S × S → [0,1], Δ_λ(d) ⊑Γ^_λ(d). The only nontrivial case is when s,t ∉ A, which follows by definition of 𝒦_d, by noticing that [d]st≤[d]st and that the maximum is order preserving.
It remains to prove that the minimum is attained. To this end, define a coupling structure ^* as ^*(s,t) = ω_s,t, for s,t ∉ A, where ω_s,t∈τ(s)τ(t) is such that 𝒦_(τ(s),τ(t)) = ∑_u,v(u,v) ·ω_s,t(u,v). By construction, = Γ^^*_λ(), hence C^*⊑. Since ^* is a coupling structure for , by what we have shown above we also have ⊑C^*. Therefore, = C^*.




§ GREEDY COMPUTATION OF THE BISIMILARITY DISTANCE
 

Inspired by the characterization given in Theorem <ref>, we propose a procedure to compute the bisimilarity pseudometric that is alternative to those previously described. 

The set of coupling structures for  can be endowed with the preorder  defined as ' iff C⊑C'. 
Theorem <ref> suggests to look at all the coupling structures  for  in order to find an optimal one, i.e., minimal w.r.t. . However, it is clear that the enumeration of all the couplings is unfeasible, therefore it is crucial to provide an efficient search strategy which allows one to find an optimal coupling by exploring only a finite amount of them. Moreover we also need an efficient method for computing the λ-discrepancy associated with a coupling structure.



 §.§ Computing the Lambda-Discrepancy
 
In this section we consider the problem of computing the λ-discrepancy associated with a coupling structure. 

By Tarski's fixed point theorem, C corresponds to the least pre-fixed point of Γ^_λ, that is C = { d ∈ [0,1]^S × S|Γ^_λ(d) ⊑ d }. This allows us to compute the λ-discrepancy associated with  as the optimal solution of the following linear program, denoted by 𝐷𝑖𝑠𝑐𝑟_λ().

    *arg min_d    ∑_s,t ∈ S d_s,t
        d_s,t≥ 1       if s ≢_A t
        d_s,t≥st      if s ≡_A t
        d_s,t≥λ( st  + (1 - st) ·∑_u,v ∈ S d_u,v·(s,t)(u,v) )       if s,t ∉ A

𝐷𝑖𝑠𝑐𝑟_λ() has a number of inequalities that is bounded by 2 |S|^2 and |S|^2 unknowns, thus, it can be efficiently solved using the interior-point method. 

 
If one is interested in computing the λ-discrepancy for a particular pair of states (s,t), the method above can be applied on the least independent set of inequalities containing the variable d_s,t. Moreover, assuming that for some pairs the values associated to d are known, the set of constraints can be further decreased by substitution.







 §.§ Greedy Strategy for Optimal Coupling Structures
 


In this section, we propose a greedy strategy that moves toward an optimal coupling structure starting from any given one. Then, we provide sufficient and necessary conditions for a coupling structure, to ensure that its associated λ-discrepancy coincides with .

Hereafter we fix a CTMC =(S, A, τ, ρ, ℓ) and a coupling structure  for it. 
The greedy strategy takes a coupling structure and locally updates it at a given pair of states in such a way that it decreases it with respect to . 
For s,t ∉ A and ω∈τ(s)τ(t), we denote by [(s,t)/ω] the update of  at (s,t) with ω, defined as [(s,t)/ω](u,v) = (u,v), for all (u,v) ≠ (s,t), and [(s,t)/ω](s,t) = ω; it is worth noting that, by construction, [(s,t)/ω] is a coupling structure of .

The next lemma gives a sufficient condition for an update to be effective for the strategy.
 
Let s, t ∉ A and ω∈τ(s)τ(t). Then, for ℋ = [(s,t)/ω] and any λ∈ (0,1), if Γ^ℋ_λ(C)(s,t) < C(s,t) then H⊏C.

 
It suffices to show that Γ^ℋ_λ (C) ⊏C, i.e., that C is a strict post-fixed point of Γ^ℋ_λ. Then, the thesis follows by Tarski's fixed point theorem. 

Let u,v ∈ S. If u ≢_A v, then Γ^ℋ_λ(C)(u,v) = 1 = 
Γ^_λ(C)(u,v) =C(u,v). If u,v ∈ A, then Γ^ℋ_λ(C)(u,v) = uv = Γ^_λ(C)(u,v) =C(u,v). If u,v ∉ A and (u,v) ≠ (s,t), by definition of ℋ, we have that (u,v) = ℋ(u,v), hence
Γ^_λ(C)(u,v) = Γ^ℋ_λ(C)(u,v). The remaining case, i.e., (u,v) = (s,t), holds by hypothesis. This proves Γ^ℋ_λ (C) ⊏C.


Lemma <ref> states that 𝒞 can be improved w.r.t.  by updating it at (s,t), if s,t ∉ A and there exists a coupling ω∈τ(s)τ(t) such that the following holds

    ∑_u,v ∈ SC(u,v) ·ω(u,v) <
      ∑_u,v ∈ SC(u,v) ·(s,t)(u,v)  .

A coupling that enjoys the above condition is ω∈ TP(C, τ(s), τ(t)) where, for arbitrary μ,ν∈(S) and c  S × S → [0,1]

    TP(c, μ, ν) = *arg min_ω    ∑_s,t ∈ S c(u,v) ·ω_u,v
       ∑_vω_u,v = μ(u)       ∀ u ∈ S 
       ∑_uω_u,v = ν(v)       ∀ v ∈ S 
       ω_u,v≥ 0       ∀ u,v ∈ S   .

The above problem is usually referred to as the (homogeneous) transportation problem with μ and ν as left and right marginals, respectively, and transportation costs c. This problem has been extensively studied and comes with (several) efficient polynomial algorithmic solutions <cit.>.

This gives us an efficient solution to update any coupling structure, that, together with Lemma <ref> represents a strategy for moving toward  by successive improvements on the coupling structures.

Now we proceed giving a sufficient and necessary condition for termination.

 

If C≠, then there exist s,t ∉ A and a coupling structure ℋ = 𝒞[(s,t)/ω] for ℳ such that Γ^ℋ_λ(C)(s,t) < C(s,t).


We proceed by contraposition. If for all s,t ∉ A and ω∈τ(s)τ(t), Γ^ℋ_λ(C)(s,t) ≥C(s,t), then C = Δ_λ(C). Since, by Theorem <ref>, Δ_λ has a unique fixed point, C =.


The above result ensures that, unless 𝒞 is optimal w.r.t  
, the hypothesis of Lemma <ref> is 
satisfied, so that, we can further improve 𝒞 as aforesaid.



The next statement proves that this search strategy is correct.
 

= C iff there is no coupling ℋ for
ℳ such that Γ^ℋ_λ(C) ⊏C.


We prove: ≠C iff 
there exists ℋ such that Γ^ℋ_λ(C) ⊏C.
(⇒) Assume ≠C. 
By Lemma <ref>, there exist a pair of states s,t ∈ S and a coupling ω∈τ(s)τ(t) such that λ·∑_u,v ∈ SC(u,v) ·ω(u,v) < C(s,t). As in the proof of Lemma <ref>, we have that ℋ = 𝒞[(s,t)/ω] 
satisfies Γ^ℋ_λ (C) ⊏C.
(⇐) Let ℋ be such that 
Γ^ℋ_λ(C) ⊏C. 
By Tarski's fixed point theorem D⊏C. 
By Theorem <ref>, ⊑D⊏C.




 
Note that, in general there could be an infinite number of coupling structures for a given CTMC. However, for each fixed d ∈, the linear function mapping ω to ∑_u,v ∈ S d(u,v) ·ω(u,v) achieves its minimum at some vertex in the transportation polytope τ(s)τ(t). 
Since the number of such vertices is finite, the termination of the search strategy is ensured by updating the coupling structure using optimal vertices. This does not introduce further complications in the algorithm since there are methods for solving the transportation problem (e.g., Dantzig's primal simplex method <cit.>) which provide optimal transportation schedules that are vertices. 



§ THE ON-THE-FLY ALGORITHM
 


In this section we describe an on-the-fly technique for computing the bisimilarity distance  fully exploiting the greedy strategy of Section <ref>.





Let Q ⊆ S × S and consider the problem of computing (s,t) for all (s,t) ∈ Q. Recall that the strategy proposed in Section <ref> consists in a traversal _0 ▹_λ_1 ▹_λ⋯▹_λ_n of the set of coupling structures for  that starts from an arbitrary coupling structure _0 and leads to an optimal one _n. 
We observe that, for any i < n

  
  * the improvement of each coupling structure _i is obtained by a local update at some pair of states u,v ∉ A, namely _i+1 = _i[(u,v)/ω] for some ω∈ TP(_i,τ(u),τ(v));
  
  * the pair (u,v) is chosen according to an optimality check that is performed locally among the couplings in τ(u)τ(v), i.e., _i(u,v) ∉ TP(_i,τ(u),τ(v));
  
  * whenever a coupling structure _i is considered, its associated λ-discrepancy _i can be computed by solving the linear program 𝐷𝑖𝑠𝑐𝑟_λ(_i) described in Section <ref>.

Among the observations above, only the last one requires to look at the coupling structure _i.
However, as noticed in Remark <ref>, the value _i(s,t) can be computed without considering the entire set of constraints of 𝐷𝑖𝑠𝑐𝑟_λ(_i), but only the least independent set of inequalities that contains the variable d_s,t. 
Moreover, provided that for some pairs of states E ⊆ S × S the value of the distance is known, the linear program 𝐷𝑖𝑠𝑐𝑟_λ(_i) can be further reduced by substituting the occurrences of the unknown d_u,v by the constant (u,v), for each (u,v) ∈ E.

This suggests that we do not need to store the entire coupling structures, but they can be constructed on-the-fly during the calculation.

Specifically, the couplings that are demanded to compute _i(s,t) are only those _i(u,v) such that , where __i,E^* is the reflexive and transitive closure of __i,E⊆ S^2 × S^2, defined by

    (s',t') __i,E (u',v')    iff   _i(s',t')(u',v') > 0  and  (u',v') ∉ E  .




The computation of the bisimilarity pseudometric is implemented by 
Algorithm <ref>. It takes as input a finite CTMC 
= (S, A, τ, ρ, ℓ), a discount factor λ∈ (0,1), and a query set Q ⊆ S × S.
We assume the following global variables to store: 

  
  * 𝒞: the current (partial) coupling structure;
  
  * d: the λ-discrepancy associated with 𝒞;
  
  * ToCompute: the pairs of states for which the distance has to be computed; 
  
  * Exact: the set pairs of states (s,t) such that d(s,t) = (s,t), i.e., those pairs which do not need to be further improved[Actually, the set Exact contains those pairs such that |d(s,t) - (s,t)| ≤ϵ where ϵ corresponds to the precision of the machine. In our implementation ϵ = 10^-9.];
  
  * Visited: the set of pairs of states that have been visited so far.

Moreover, s,t will denote the set { (u,v) | (s,t) _,Exact^* (u,v) }.

At the beginning  (line <ref>–<ref>) both the coupling structure 𝒞 and the 
discrepancy d are empty, there are no visited states, no exact computed distances, and the pairs to be computed are those in the input query.

While there are still pairs left to be computed (line <ref>), we pick one (line <ref>), say (s,t). According to the definition of , if s ≢_A t then (s,t)=1; if s = t then (s,t) = 0 and if s,t ∈ A then (s,t) = st, so that, d(s,t) is set accordingly, and (s,t) is added to 𝐸𝑥𝑎𝑐𝑡 (lines <ref>–<ref>). Otherwise, if (s,t) was not previously visited, a coupling ω∈τ(s)τ(t) is guessed, and the routine 𝑆𝑒𝑡𝑃𝑎𝑖𝑟 updates the coupling structure  at (s,t) with ω (line <ref>), then the routine 𝐷𝑖𝑠𝑐𝑟𝑒𝑝𝑎𝑛𝑐𝑦 updates d with the λ-discrepancy associated with  (line <ref>). According to the greedy strategy,  is successively improved and d is consequently updated, until no further improvements are possible (lines <ref>–<ref>). Each improvement is obtained by replacing a sub-optimal coupling C(u,v), for some (u,v) ∈s,t, by one taken from TP(d, τ(u), τ(v)) (line <ref>).
Note that, each improvement actually affects the current value of d(s,t), since the update is performed on a pair in s,t. It is worth to note that  and 𝐸𝑥𝑎𝑐𝑡 are constantly updated, hence s,t may differ from one iteration to another.

When line <ref> is reached, for each (u,v) ∈s,t, we are guaranteed that d(u,v) = (s,t), therefore s,t is added to 𝐸𝑥𝑎𝑐𝑡 and, for these pairs, d will no longer be updated. At this point (line <ref>), the couplings associated with the pairs in 𝐸𝑥𝑎𝑐𝑡 can be removed from .
In line <ref>, the exact pairs computed so far are removed from ToCompute. 
Finally, if no more pairs need be considered, the exact distance on Q is returned (line <ref>).

Algorithm <ref> calls the subroutines 𝑆𝑒𝑡𝑃𝑎𝑖𝑟 and
𝐷𝑖𝑠𝑐𝑟𝑒𝑝𝑎𝑛𝑐𝑦. The former is used to construct and update the coupling structure 𝒞,
the latter to update the current over-approxi­ma­tion d during the computation.
Next, we explain how they work. 



𝑆𝑒𝑡𝑃𝑎𝑖𝑟 (Algorithm <ref>) takes as input a CTMC = (S, A, τ, ρ, ℓ), a pair of states s,t ∈ S, and a coupling ω∈τ(s)τ(t). In lines <ref>–<ref>, the coupling structure  is set to ω at (s,t), then (s,t) is added to 𝑉𝑖𝑠𝑖𝑡𝑒𝑑. The on-the-fly construction of the coupling structure is recursively propagated to the demanded successor pairs of (s,t) according to the information accumulated so far. 
During this construction, if some states with trivial distances are encountered, d and
Exact are updated accordingly (lines <ref>–<ref>).



𝐷𝑖𝑠𝑐𝑟𝑒𝑝𝑎𝑛𝑐𝑦 (Algorithm <ref>) takes as input
a discount factor λ∈ (0,1) and a pair of states s,t ∉ A. It constructs the least linear program obtained from 𝐷𝑖𝑠𝑐𝑟_λ(), that can compute (s,t) using the information accumulated so far (line <ref>). In lines <ref>–<ref> the current λ-discrepancy is updated accordingly; and those pairs (u,v) ∈s,t for which the current λ-discrepancy coincide with the distance are added to 𝐸𝑥𝑎𝑐𝑡.





Next, we present a simple example of Algorithm <ref>, 
showing the main features of our method:
(1) the on-the-fly construction of the (partial) coupling, and (2) the restriction only to those variables 
which are demanded for the solution of the system of linear equations.

[On-the-fly computation] 
Consider the CTMC in Figure <ref>, and assume we want to compute the λ-discounted bisimilarity distance between states s_1 and s_4, for λ = 1/2.

Algorithm <ref> starts by guessing an initial coupling structure 𝒞_0. This is done by considering only the pairs of states which are really needed in the computation. Starting from the pair (s_1,s_4) a coupling in ω_1,4∈τ(s_1)τ(s_4) is guessed as in Figure <ref> and assigned to _0(s_1,s_4). 
This demands for the exploration of the pairs  (s_2,s_3), (s_2,s_4), (s_1,s_2) and the guess of three new couplings ω_2,3∈τ(s_2)τ(s_3), ω_2,4∈τ(s_2)τ(s_4), and ω_1,2∈τ(s_1)τ(s_2), to be associated in _0 with their corresponding pairs. Since no other pairs are demanded, the construction of _0 terminates as shown in Figure <ref>. The λ-discrepancy associated with _0 for the pair (s_1, s_4) is obtained as the solution of the following reduced linear program  

    *arg min_d     ( d_1,4 + d_2,3 + d_2,4 + d_1,2 ) 
       d_1,4≥1/6
       d_1,4≥α/2 + (1-α)/2·( 4/9· d_2,3 + 17/63· d_2,4 + 11/63·d_4,4^=0) 
       d_2,3≥1/2
       d_2,3≥α/2 + (1-α)/2·( 1/4· d_1,2 + 1/4·d_2,2^=0 + 1/2· d_2,4) 
       d_2,4≥2/3
       d_2,4≥1/2·( 1/9· d_1,2 + 5/36· d_1,4 + 4/9· d_2,3 + 11/36· d_2,4) 
       d_1,2≥1/2
       d_1,2≥α/2 + (1-α)/2·( 5/7·d_2,2^=0 + 1/4· d_1,4 + 1/28· d_2,4)

where α = [15][9] = √(3/5)/25 (by Equation (<ref>)). Note that, the bisimilarity distance for the pairs (s_2,s_2) and (s_4,s_4) is always 0, thus d_2,2 and d_4,4 are substituted accordingly. The solution of the above linear program is 
d^_0(s_1,s_4) = α/2 + 5(1-α)/21, d^_0(s_2,s_3) = 1/2, d^_0(s_2,s_4) = 2/3, and d^_0(s_1,s_2) = 1/2.

Since, the λ-discrepancy for (s_2,s_3), (s_2,s_4), and (s_1,s_2) equals the distance  between their labels, it coincides with the bisimilarity distance, hence it cannot be further decreased. Consequently, the pairs of states are added to the set 𝐸𝑥𝑎𝑐𝑡 and their associated couplings are removed from _0. Note that, these pairs will no longer be considered in the construction of a coupling structure.

In order to decrease the λ-discrepancy of (s_1,s_4), Algorithm <ref> constructs a new coupling structure _1. According to our greedy strategy, _1 is obtained from _0 updating _0(s_1,s_4) (i.e., the only coupling left) by the coupling ω'_1,4∈τ(s_1)τ(s_4) (shown in Figure <ref>) that is obtained as the solution of a transportation problem with marginals τ(s_1) and τ(s_4), where the current λ-discrepancy is taken as cost function.
The resulting coupling does not demand for the exploration of new pairs in the CTMC, hence the construction of _1 terminates. The reduced linear program associated with _1 is given by 

    *arg min_d     d_1,4
       d_1,4≥1/6
       d_1,4≥α/2 + (1-α)/2·( 
    	1/9·d_2,2^=0  + 
    	4/9·d_2,3^=1/2  + 
    	10/63·d_2,4^=2/3  + 
    	2/7·d_4,4^=0)

whose solution is d^_1(s_1,s_4) = α/2 + 31(1-α)/189.

Solving again a new transportation problem with the improved current λ-discrepancy as cost function, we discover that the coupling structure _1 cannot be further improved, hence we stop the
computation, returning (s_1,s_4) = d^_1(s_1,s_4) = α/2 + 31(1-α)/189.


 
Algorithm <ref> can also be used for computing over-appro­xi­mated 
distances. Indeed, assuming over-estimates for some particular distances are already known, they
can be taken as inputs and used in our algorithm simply storing them in the variable d and treated
as “exact” values. In this way our method will return the least over-approximation
of the distance agreeing with the given over-estimates.
This modification of the algorithm can be used to further decrease the exploration of the CTMC.
Moreover, it can be employed in combination with approximated algorithms, having
the advantage of an on-the-fly state space exploration.



§ EXPERIMENTAL RESULTS
 

In this section, we evaluate the performance of the on-the-fly algorithm on a 
collection of randomly generated CTMCs[
The tests have been performed on a prototype implementation coded
in Wolfram Mathematica^  9 (available at 
<http://people.cs.aau.dk/ giovbacci/tools.html>)
running on an Intel Core-i7 3.4 GHz processor with 12GB of RAM.].



First, we compare the execution times of the on-the-fly algorithm with those of the iterative method
proposed in Section <ref>. 
Since the iterative method only allows for the computation of the distance for all state pairs at once,
the comparison is (in fairness) made with respect to runs of our on-the-fly algorithm with input 
query being the set of all state pairs. For each input instance, the comparison involves the following steps:
