Bisimulation Metrics are Optimal Transport Distances, and Can be Computed Efficiently
Authors: Sergio Calo
Anders Jonsson
Gergely Neu
Ludovic Schwartz
Javier Segovia-Aguas
Venue: NeurIPS 2024 poster (2024)
Keywords: Optimal transport
Markov chains
bisimulation metrics
Markov decision processes
Primary area: reinforcement_learning

Abstract: We propose a new framework for formulating optimal transport distances between Markov chains. Previously known formulations studied couplings between the entire joint distribution induced by the chains, and derived solutions via a reduction to dynamic programming (DP) in an appropriately defined Markov decision process. This formulation has, however, not led to particularly efficient algorithms so far, since computing the associated DP operators requires fully solving a static optimal transport problem, and these operators need to be applied numerous times during the overall optimization process. In this work, we develop an alternative perspective by considering couplings between a ``flattened'' version of the joint distributions that we call discounted occupancy couplings, and show that calculating optimal transport distances in the full space of joint distributions can be equivalently formulated as solving a linear program (LP) in this reduced space. This LP formulation formulation allows us to port several algorithmic ideas from other areas of optimal transport theory. In particular, our formulation makes it possible to introduce an appropriate notion of entropy regularization into the optimization problem, which in turn enables us to directly calculate optimal transport distances via a Sinkhorn-like method we call Sinkhorn Value Iteration (SVI). We show both theoretically and empirically that this method converges quickly to an optimal coupling, essentially at the same computational cost of running vanilla Sinkhorn in each pair of states. Along the way, we point out that our optimal transport distance exactly matches the common notion of bisimulation metrics between Markov chains, and thus our results also apply to computing such metrics, and in fact our algorithm turns out to be significantly more efficient than the best known methods developed so far for this purpose.

Bisimulation Metries are Optimal Transport Distances, and Can be Computed Efficiently
Sergio Calo Anders Jonsson Gergely Neu
Ludovic Schwartz Javier Segovia-Aguas
Universitat Pompeu Fabra, Barcelona, Spain
{sergio.cal calo, anders. jonsson gergely. neu, ludovic wwwwartz,aaaiaaiaae www..a..e..r.riaaaaaa.....aaa..a @upf edu
Abstract
We propose new framework for formulating optimal transport distances between Markov chains. Previously known formulations studied couplings between the en- tire joint distribution induced by the chains, and derived solutions via a reduction to dynamic programming (DP) in an appropriately defined Markov decision process. This formulation has, however, not led to particularly efficient algorithms so far, since computing the associated DP operators requires fully solving static optimal transport problem, and these operators need to be applied numerous times during the overall optimization process. In this work, we develop alternative perspective by considering couplings between a "flattened" version of the joint distributions that we call discounted occupancy couplings, and show that calculating optimal transport distances the full space of joint distributions can be equivalently formu- lated solving linear program (LP) in this reduced space. This LP formulation allows us to port several algorithmic ideas from other areas of optimal transport theory. In particular, our formulation makes possible introduce an appropriate notion of entropy regularization into the optimization problem, which in turn en- ables us to directly calculate optimal transport distances via a Sinkhorn-like method we call Sinkhorn Value Iteration (SVI). We show both theoretically and empirically that this method converges quickly to an optimal coupling, essentially the same computational cost of running vanilla Sinkhorn in each pair of states. Along the way, we point out that our optimal transport distance exactly matches the common notion of bisimulation metrics between Markov chains, and thus our results also ap- ply to computing such metrics, and in fact our algorithm turns out to be significantly more efficient than the best known methods developed so far for this purpose.
Introduction
Measuring distances between structured objects and sequences is an important problem in a variety of areas of science. The more structured the objects become, the harder it gets to define appropriate notions of distances, good notions of proximity need to take into account the possibly complex relationships between the constituent parts each object. The possibility that the objects in question may be random further complicates the picture, and in such cases becomes more natural to measure distances between the underlying joint probability distributions. Within the specific context of comparing stochastic processes, two natural notions of distance have emerged over the past decades: the notion of probabilistic bisimulation metrics that takes its root in modal logic and theoretical computer science (Sangiorgi 2009 Abate 2013 and the notion of cerimmal-transpon distances that originates from probability theory Villani 2009 Peyré and Cuturi] 2019]. In this paper, we show that bisimulation metrics are in fact optimal-Iranssoon distances, and make use of this observation to derive efficient algorithms for computing distances between stochastic processes.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

The two distance notions have found strikingly different applications. Bisimulation emerged within the area of theoretical computer science as one the most important important concepts in concurrency theory and formal verification of computer systems [Park] 1981 Milner 1989]. and has been extended to probabilistic transition systems by Larsen and Skou 19891. Within machine learning, bisimulation metrics have become especially popular in the context of reinforcement learning (RL) due to the work Ferns et al. [2004), and have become one of the few standard tools of representation learning (Jiang] 2018 2024]. In particular, the work of et al. [2004) advocates for using bisimilarity for state aggregation, measuring similarities of states in terms of similarities of two chains Mx and My that only differ in their initial state. While this approach has inspired numerous follow-up works [Gelada et al., 2019 Castro 2020 Agarwal et al., 2021b Zhang et al.) 2021 Hansen- Estruch et al., 2022 Castro et al.] 2022], ,, ultimately this line of work has failed to discover efficient algorithms for computing bisimulation metrics and has largely resorted to heuristics for computing similarity metrics.
basis
On the other hand, optimal transport (OT) has found numerous applications in areas as diverse as economics [Galichon [016]. signal processing [Kolouri et al.] 2017]. or genomics [Schiebinger ct al. 2019 Within machine learning, it has been used for the similarly diverse areas of domain
adaptation Courty et al.] 2016), generative modeling [Arjovsky et al. 2017 Song et al. 2020 Shi
al. 2024], representation learning [Courty et al., 20188, and, perhaps most relevant to our work, as
a way of measuring distances between graphs Titouan et al.] 2019 Chen et al. 2022 Chuang and
Jegelka 2022]. The recent works Yi et al. [2021), I, Brugère et al. 2024] propose to define graph
distances via studying the behavior of random walks defined on the graph, thus reducing the problem of comparing graphs to comparing stochastic processes- exactly the subject of the present paper. Other applications of OT between stochastic processes include generative modeling for sequential data [Xu al. 2020), pricing and hedging in mathematical finance (Backhoff- Bacckhoff-Veraguss et al. 2017]. ,
and analyzing multistage stochastic optimization problems [Pflug 2010 Bartl and Wiesel 2022
It appears however that the literature on optimal transport for stochastic processes has apparently not yet discovered connections with bisimulation metrics and the rich intellectual history behind it. Also, applications of optimal transport for representation learning within the context of reinforcement learning appear to be nonexistent.
In this paper observe that, despite their apparent differences, bisimulation metrics and optimal transport distances are one and the same. Furthermore, we provide a new perspective on both OT distances and bisimulation metrics by formulating the distance metric as the solution of linear program (LP) in the space of "occupancy couplings", a dimensoammniional projection of the infinite- dimensional process laws. Building on tools from computational optimal transport (Peyré and Cuturi 2019] and reguopy-regularrized Markov decision processes (Neu et al., 2017 Geist et al. 2019]. we design algorithm that effectively combines Sinkhorn's algorithm [Sinkhorn and Knopp 1967 Cuturi 2013) with an mntropy-regularizee
version of the classic Value Iteration algorithm [ee man 1957 Neu et al. 2017]. Building on recent work on computational optimal transport [Altschller cet al. 2017 Ballu and Berthet 2023), we provide theoretical guarantees for the resulting algorithm (called Sinkhorn Value Iteration) and perform numerical studies that demonstrate its effectiveness for computing distances between Markov chains.
Notations. For a finite set S, we use As to denote the set of all probability distributions over S. We will denote infinite sequences by = (zg, .1....) and the corresponding subsequences as Tn = (1 ........ ). For two sets X' and V, we will often write X'Y' to abbreviate the direct- product notation x X), and for two indices x and y and a function f XY - Z, we will often write (zy) instead of [ar, y) to save space. Also, we will denote scalar products by (.-) and use Il-1l, to denote the nn-norm.
2 Preliminaries
We study the problem of measuring distances between pairs of finite Markov chains. Specifically, we consider two stationary Markov processes Mx = Px, x0,x) and My (),, Py,,y,,,.)), where
X and y are the finite state spaces with cardinalities m = | and n = =YY,
Px : X + Ax and Py: : y + Ay are the transition kernels that determine the evolution of the states as Px (2 [x)) = P P[X++1 =x'| X.= and Py(1//1) IP [Y+ = Y.=w) for all i, and
Vo,x and vo.y are the initial-state distributions with Xo vo.x and Yo ~ vo.y-

---

Without significant loss of generality, we will suppose that the initial states are fixed almost surely as Xo 2n and Yn = yo. and refer their corresponding joint distribution as v/- S These objects together define a sequence of joint distributions IP [Xo, XI...., (.....,,,) X.) = (ro, F1,.... , and P ((Y,, .,...., Vn,) (Vo. U1..... ...)..,,,))] for each n, which together define respective the laws of the infinite sequences X (X1, X2....) and Y (Y1, Y2....) via Kolmogorov's 's extension theorem. With slight abuse of notation, we will use M, and My to denote the corresponding measures that satisfy Mx(xx)) = IP X, En and My(U..) = P [Yx Jn. for any I E x y E you and n. The corresponding conditional distributions are denoted as Mx(In|Tn = P [ = In X = -n-1] and My(yalu. = P
[Yr, n = Yn Fn = Vn 1-
2.1
Optimal transport between Markov chains
Our main object of interest in this work is notion of optimal transport distance between infinite- horizon Markov chains. Several previous works have studied such distances (which are discussed in detail in Appendix A), and our precise definition we give below is closest to Moulos (202T), O D'Conoor et al. [2022] and Brugère et al. [2024). We consider Markov chains on state spaces where "ground metric" (or "ground cost")
X' x y + IR is available to measure distances between any two individual states x X and y E y. with the distance denoted as c(x, y). For any two sequences
F (I1, 02,....... and J (y1, 1/2 ...) we define the discounted total cost
(T,
Σ (e(w,,8));
where E (0. 1) is the discount factor that expresses the preference that two sequences be considered further apart if they exhibit differences at earlier times in terms of the ground cost c. Following the optimal-transport literature, we will consider distances between the stochastic processes Mx and My via the notion of couplings. To this end, we define coupling of Mx and My as a stochastic process evolving on the joint space X,,, with its law defined for all n as Mxy (I,Tn) = P[Xn In, Y Vn]. required to satisfy Σy, EYmh (xy(X) = (x( and Ex,ex (xy(Fn = (y(w,)- We will define the set of all such couplings as I1.
The notion of couplings defined above does not respect the temporal structure of the Markov chains Mx and My appropriately: while by definition the distribution of state X,, may only be causally influenced by past states X& with hs < n, the general notion of coupling above allows the state X,, to be influenced by future states Y& with works mentioned in the introduction), we will introduce the notion of bicausal couplings. A coupling Mxy called bicausal and only if it satisfies
n as well. To rule out this possibility (and following past
Σ Mx Zn 17. 17n-1) = Mx (Fa| and Vn
Σ Mx IF
= My(Uvl
for all sequences x, y E x x y and all n. Denoting the set of all bicausal couplings by Ih we define our optimal transport distance as
f f (XXX
(M My; c, zo. vo)
inf MxyElle
where we emphasize the dependence of the distance on zn, yo explicitly with our notation.
By noticing that the optimization problem outlined above can reformulated as a Markov decision process (MDP), Moulos (2021) has shown that the infimum in (2) is achieved within the family of Markovian couplings that satisfy
Σ Mxy(Tnyn Fn. - = Px(xxnEnn and Σ May(EENN 117-)) = Py(ynIy Vn
for all sequences of state pairs and all values of n. Furthermore, can seen that Markovian couplings are fully specified in terms of transition couplings of the form XY -> Axy, with T (2' v/yz standing for P [Xx++1, Y(+1) = (x',g') (Xz, Yt) = (Ie, yt)] under the law induced by the coupling. We say that transition coupling is valid if it satisfies the marginal constraints Ey ((x'''(ry) = (x'(() and Enr (x')'(xy) = Py(s/ lw).. Defining the set of such valid tran-
sition couplings in state pair zy by Hey {o E Axy Ey p(e'y') = Px (x), Ene Saappay' (z'v')
3

---

Py((/lv)), Moulos [2021] introduces an MDP M with an infinite action set corresponding to picking the joint next-state couplings Hyy. An optimal transition coupling can be found by solving the following Bellman optimality equations of the MDP M:
V" = c(zy) + inf Σ Eple'v'V
(3)
petlay x'y'
The infimum on the right-hand side is achieved by an optimal transition coupling π- (+(aay = arg min Exy (x'y') The solution V is unique and can be shown satisfy v*(xy) W,, (Mx, My; C, 2, y) for all zy E XY. For completeness, we include the precise definition of the MDP M and the proofs of these results in Appendix B
2.2 Bisimulation metrics
The notion of bisimulation metrics has been introduced by Desharnais et al. [1999] 2004 and van Breugel and Worrell [200T], with the purpose of defining distances between Markov chains, using methodology rooted in modal logic that at first may appear entirely different from the optimal-traspoott framework described above. We only give a very high-level overview of the classic logic-based characterization here (as the fine details are irrelevant to the final conclusion that this section is headed to), and refer the reader to the additional discussion in Appendix A for further reading. Desharnais et al. [9999 2004] considered labeled Markov chains where a labeling function r: X UU + I2 assigns labels to each state, and defined bisimulation metrics via
d. zo. yo) sup JEF,
If (zo) (y0)),
(4)
where F., is family of functional expressions generated by certain grammar, and fM, : X + IR
and fMy : y/- + R are the respective "interpretations" of each f E F on the Markov chains Mx and My. Building on this formulation, /van Breugel and Worrell (2001) have shown that the distance metric can be equivalently characterized by the solution of a fixed-point equation, whose expression was subsequently used by Ferns et al. (2004) to define bisimulation metrics for Markov decision processes where r takes the role of a reward function, and the evolution of states may be influenced by actions. The case of having no actions available corresponds to our setting, where their definition of a bisimulation metric simplifies to the solution of the fixed-point equation
U ('yy = (1 - >))(()) - - r(y)l + inf
Enldrs pEllry z'y
(5)
The solution to this system is unique and satisfies (ryyy) = dM My;r, Putting this result side-by-side with Equation (3), one can immediately realize that bisimulation metrics and our notion of optimal transport distances coincide when picking the ground cost function c(ay) = (1=)) r(y)). To our knowledge, this remarkable observation has not been publicly made anywhere in either the optimal-transport or the bisimulation-ned -etric literature, This connection has several important implications, which we have already discussed at some length in Section [ We relegate further discussion of these metrics in the light of this observation to Appendix
3 Optimal transport between Markov chains as a linear program
Optimal transport problems can typically be formulated as linear programs (LPs), since couplings can be characterized as joint distributions satisfying set of easily-expressed linear constraints (see, e.g., Chapter 3 in Peyré and Cuturi 2019). Our problem is no exception, and in fact the original problem statement of Equation 2 can be expressed in this form: the constraints defining 1le are all lincar in Mxy. However, Mxy is an anfinite-dimensiona object and thus this formulation is not instructive for developing computationally tractable algorithms. We address this problem in this section, where we define equivalent LP formulation that replaces the infinite-dimension optimization variable with an appropriate low-dimensional projection. Our framework builds on the classic LP formulation of optimal control in MDPs first proposed in the 1960's [Manne] 1960) de Ghellinck, 1960 d' ['penoox]] 1963 Denardo 1970], and covered thoroughly in several standard textbooks (e.g., Section 6.9 of Puterman 1994).
In order to set things up, we need to start with some important definitions. We say that transition coupling π x' Axy generates trajectory (Xo, Yo, X1, YI....) if (Xo. Yo) vo and the

---

subsequent state-pairy are drawn independently from the transition coupling as (X1+1, Yo+1) ~ (-I|,, Y,) for all Then, we define the occupancy coupling μº associated with this process as distribution over XY x XY, with each of its entries xy,r'y defined as
- 7) Σ (GN
where E, [-] emphasizes that the trajectory of state-pairs has been generated by я. In words, u² (xy, x' '²) is the discounted number of times that the quadruple (ry.x'y' is visited by the process. With this definition, it easy to notice that the objective optimized in Equation (2) can be rewritten as linear function of μ², Indeed, suppose that Mxy is the law of the process generated by as described above, so that we can write
/
(X,
= E
Σ 1=
C'dcX Y1)
= =EA
ET
Σ
E Tux. ((xxx)))caw))(lx))
Σ xV
E
80
((x d(ry)
1
Σ ay.r'y
a' 2'g (40.c)
1-7
1-7
We say that an occupancy coupling μ is valid if it is generated by valid transition coupling. It is easy to verify that every valid occupancy coupling E Axyxxy satisfies the following three constraints:
Σ
Culry, = Σ
El( y/ myy + (oyy)
(Vay (X),,
(6)
x'y
x'yy'
Σ
u(xy,, d'1)') = Σ
plry, x² Px(1(a)))) (Vx,s' EX xxx))
(7)
y
x²
Σ
u(zy,,'
Σ x²
play, x By) (Vu,3' x))
(8)
We refer to the first equality constraint as the flow constraint and the second and third ones as the transition coherence constraints for Mx and My, respectively We show in the following lemma that the above conditions uniquely characterize the set of valid occupancy couplings.
Lemma 1. A distribution μ E Axyxxy is a valid occupancy coupling associated with some transition coupling x : XY - Axy if and only if satisfies Equations
(8)
One important consequence of this result is that for each transition coupling, one can associate a transition coupling Ty with f= n., with entries satisfying Exty" "(ry,,"""") (x'y (my) plzy, x'y²" y/). A proof is provided in Appendix B.2 Having established in the previous section that stationary occupancy couplings are sufficient to achieve the supremum in the definition of the OT distance of Equation
(2). this result immediately implies the following.
Theorem 1. Let B, Sx and Sy respectively stand for the sets of all distributions μa E Axyxxy that
satisfy equations ( 7 and 8 Then,
inf W (Mx, My; c, no. 30) 1 - WEBnSxnSy (n, c).
Since the constraints on le in the above reformulation are all linear, the optimization problem stated above is obviously linear program. Notably, the resulting I.P is not the standard dual to the IPP associated with the MDP formulation introduced in Section [2. which would result in an infinite-dimensionall LP with one constraint associated with each of the continuously-vai
actions. Such an infinite-dimensions reformulation was previously considered by CChen et al. [2012] in the context of computing bisimulation metrics, who used design. Notably, our formulation results in imite-dimensional LP with [x p2 variables and |x| + | |>| + 2 constraints. Building on the celebrated result of Ye [2011] (and similarly to the work of Chen et al. 2012 mentioned above), can show that our LP can be solved in strongly polynomial time via an appropriate adaptation of the simplex method or the classic policy iteration
tool for analysis rather than algorithm
method of Howard 1960). We propose alternative methodology in the next section.
'henn 1. the first condition is known as detailed balance within statistical physics. Altogether, the constraints closely resemble what are often called the "Bellman flow constraints" in today's RL literature.
5

---

4 Sinkhorn Value Iteration
Solving ppiimll-transpon problems via standard LP solvers (like variations of the simplex method or interior-poimt methods) is known to be empirically hard, and thus seek alternatives to this approach towards optimizing our own LP defined the previous section. In computational optimal transport, a paradigm shift was initiated by Cuturi [0033] who successfully applied entropic regularization to the classic LP objective of optimal transport, and solved the resulting optimization method through an iterative algorithm called the Sinkhorn-Knoph method (sometimes simply called Sinkhorn's algorithm, due to Sinkhorn and Knopp] 1967). This method based on finding a feasible point in two-constraint problem by alternately satisfying one and the other, and has resulted in practical algorithms that were orders of magnitude faster than all previously studied methods. Entropy regular- ization and Sinkhorn's 's algorithm have thus became the most important comeertones of computational optimal transport. Drawing on the same principles as well as the theory of Markov decision processes (Neu et al., 2017 Geist et al., 2019], we develop computationally effective algorithm for computing optimal transport distances between Markov chains below.
4.1 Formal definition: Mirror Sinkhorn the space of occupancy couplings
Our method is an adaptation of a version of Sinkhorn's algorithm called Mirror Sinkhorn, first proposed and analyzed by Ballu and Berthet [2023). This method combines miiror-descent-sttllee updates [Nemirovski and Yudin [1883] Beck and Teboulle, (0033 with alternating projections to two convex sets whose intersection corresponds the feasible set we seek to optimize over. In our adaptation, we choose the two sets as
u(a"""",,, xy) + (1 ()))(()))))))))))))))))))))))) (x(2'(x) (Vag, 2²) }.
(7 EEEnea
Bx pu: Σ a(zy,z''))
xnyy"
that can be seen to be the set of distributions y that satisfy both Equations (6) and (7. and
By {Ix- u:: EMor,NV)
Σ ula""'' ry) + +(1 (1))))())))))))))))))))))))))- (v/ll) (Vay. y
which is the set of distributions z that satisfy both Equations (6) and (8) Naturally, the intersection of the two sets corresponds to valid occupancy couplings. It remains to define appropriate notion of entropy for the purpose of regularization. Following Neu et al. [211)) we will use the conditional relative entropy defined between two joint distributions p. μ² Axyxy
M(xy x/y/)/ Exiry (x",,x"y)
H (H(l(() =
Σ xy,z'''
play, y'y) log
1 ( x'v/// Emrgo y pl (xy, x² x"y"" Tu (1 (ry) T (x'y' ylay
Σ xy.x'y'
play, x' 3/) log
It is easy exercise to show that H is Bregman divergence that is convex in its first argument (see.
e.g., Appendix A.1 Neu et al. 2017). Note however that H. (ي]])²') can be zero even if μ * μí and thus it is not strongly convex in μ.
With these ingredients, we symbolically define our algorithm as calculating the sequence of updates
{ (1, c) + H (allmn))
p4++ = arg min N
for each 1: 1..... K 1, where μn the occupancy coupling associated to the trivial coupling Ti (1(12y) = (-(x) & Py(-lv) for each state pair zy E xY., 17 > 0 is a stepsize (or learning- -rate parameter and B, is chosen to Bx odd rounds and By in even rounds. By adapting tools from the theory of regulay-ee Markov decision processes, the updates can be computed in closed form by solving system of equations closely resembling the regularized Bellman equations. In particular, we define the Bellman-Sinkhorn operators for a given transition coupling as the operators

---

TF RXxx + RXV* and T5 R>Vx + RXXX acting on functions Vx E R and
Vy E RXXx respectively as
(-((clx)) c(xx)) + [Pdde"/VVale) ( (-1
(T (xy, (xY,')) = -
log x (<''y/1 (ry) (x(2'(r) 17 y¹
exp -11
and
T<'2// (ry) 2
exp -1 (cEe)) + EV x² Py(''l))
log (cy,)') (oy,')))----
V"
Then, for odd rounds, the updates can be calculated by solving the fixed-point equations Vv= - Tx Ve. defining the shorthand Quy z'y') c(ry) + y En. Px(a" (<'"a=""VVi (z'y/. ye'e''"z)" x²).. and subsequently updating the transition coupling Ty multiplicatively as
Th (2'''' (/yy exp (-TQu(ry, T (1'1/ (zy) = Ey Tx (1'y" mry) exp (~1 Q (x (x(x'(xx)
(10)
is easy to verify that this transition coupling satisfies) Ey Th+ (1'1)) 1/ |my) = P(x' (()) The updates for even rounds are computed analogously with the roles of X and y swapped. We respectively refer -iinhhnn-iikknnn equations for
to the fixed-point equations V Txx V2 and V = =75" Vi as the Mx and My, and the functions Vi and Q& as value functions. The following proposition (proved in Appendix E. formally establishes the equivalence between the two update rules.
Proposition 1. Let p and Th41 be specified for each k as in Equations Then, 1 - holds for all k.
and
10 respectively,
4.2 Practical implementation
The algorithm described above can be seen as performing online Mirror Sinkhorn updates in each state pair zy with a sequence of cost func- tions Qk. which are computed via solving the Bellman-Sinkhorn TV are easily seen be contractive with re- spect to the supremum norm with contraction factor y (as shown by a standard calculation in- cluded in Appendix [..3). 3. these equations can be solved by an adaptation of the classic Value Iteration method of Bellman| 1957]. Concretely, we repeatedly apply the Bellman-Sinkhon op- erators until the fixed point is reached up to sufficient precision (controlled by the number of update steps m.) We call the resulting method Sinkhorn Value Iteration (SVI,, and provide its pseudocode as Algorithm T
Algorithm 1: Sinkhorn Value Iteration
Input: Px, Py. C, 11, 7/ K, m
Initialise: T1 <- Px @ Py; for k = 1.... K
do
Bellman Sinkhorn equations. Since T and
if k: is odd then
Vx - (Tx ym (T Vx: (Bx projection.)
else
Vy (75² Th m Vyy: (By projection.)
end
T update (T+)::
(Equation TO
end
K Σk= (MTT):;
TT
Tout round (Ta VTan +- evaluate( Output: Touts V /oou
|Final coupling)
Notably, while SVI is defined in its abstract form as a sequence of updates in the space of occupancy couplings P&x its implementation only works with transition couplings тk. The final output of SVI is transition coupling Tout, obtained by computing the average T
1 Ex= 1.Tx of all occupancy couplings, computing Tyg and then rounding the result to a valid transition coupling. In particular we apply a simple rounding procedure due to Altschuler et al. [2017) individually on Tz (llay) for each state-pair zor-for the full details, see Appendix F.2 Besides Touts SVI also outputs estimate of V in the form of the value function V/Ku,, as denindd in Equation (TT) in Appendix B.I This function can be computed efficiently by solving the linear system of Bellman equations V The (xy) c(ry) + Exy Tour (2'() ' (xy) Tom (z'y').
(x''y')))
A number of small simplifying steps can be made to make the algorithm easier to implement. First, instead of obtaining मyx via the computationally expensive procedure described above, one can simply run the rounding procedure on the final transition coupling T and output the result. Second,
7

---

while theoretical analysis suggests setting m = oo in order to make sure that all projection steps are perfect, such exact computation may be unnecessary and inefficient in practice, and thus (much) smaller values can be used instead. Third, for small values of vi the softmax function used in the definition of the Bellman-Sinkhonn
operator can be accurately approximated by an average with respect to (x') |zy)/ (x [a), which suggests simple alternative to the projection steps. This approximates SVI similarly as to how the Mirror Descent Modified Policy Iteration method Geist et al.] [2019] approximates the mirror descent method of et al. [0177] (see also Azar et al. 2012). The resulting method (that we refer to as Sinkhorn Policy Iteration, or SPI) is presented in detail along with its theoretical analysis in Appendix D We study effects of these implementation choices
via a sequence of experiments in Section
4.3 Convergence guarantees
The following theorem establishes a guarantee on the number of iterations necessary for V V*=((zy) (zy) be an e-accurate approximation of the transport cost (Mx, My; c, x, y) for any 2, y.
Theorem 2. Suppose that Sinkhorn Value Iteration is run for K steps with regularization pa-
rameter 17 -|[l| (1-7)9 log x| and initialized with the uniform coupling defined for each vy as Ti (a'v 1/ any) = - Then, for any Tolo E Xy, the output satisfies V (2oyo) S (M,,,, My;c, 2o. vo) + if the number of iterations at least
||cl log | 324
K>
71-?
The proof is relegated to Section C and present similar performance guarantee for SPI in Appendix D Importantly, these guarantees technically only hold when setting m = 0o, which is a limitation we discuss in more detail in Section 6 The condition that m is chosen the uniform coupling is not necessary and simply made to make the statement easier to state. A more detailed statement of the bound is provided in Appendix C.4
5 Experiments
We have conducted a range of experiments on some simple environments with the purpose of illustrating the numerical properties of our algorithms and some aspects of the distance metrics studied. Due to space restrictions, we only report a very limited subsample of the results below, and refer the reader to Appendix F for the complete suite"]
One set of experiments we report here addresses the biggest open question left behind our theory: the effect of the choice of m on the quality of the updates. For this experiment, we use the classic "4-rooms" environment first studied by Sutton et al. 1999), and run both SVI (Algorithm [ and SPI (Algorithm 2 for range of different choices of m., and fixed 0.95. The results of this study are shown in Figure [] The plots indicate that the estimates produced by both algorithms converge towards the true distance at rate that is basically unaffected by m, and in particular even value of m = 1 remains competitive. This observation is consistent across all of our experiments. Also, the output of SPI appears to converge slightly more slowly towards the optimum in this experiment, but this observation is not entirely consistent and can be likely ascribed to the fact that the learning rate was not optimized favor either algorithm in this experiment. In most experiments, the two algorithms performed very similarly, up to some small occasional differences.
We have also conducted number of experiments to illustrate the potential of optimal-transport distances for comparing Markov chains of different sizes and transition functions. In the experiment we show here, we compare two Markov chains illustrated in Figure 2 The first Markov chain Mx is a simple, nine-state "gridworld" environment, which has its initial state in the upper left corner (denoted as so) and a reward of +1 in the lower left corner (shown in blue). The second, My. is an instance of the 4-rooms environment, where each room is a rotation of the aforementioned small grid. The transition kemnels in both environments are uniform distributions over the adjacent cells in the four principal directions. The plot shows the distances between the two chains as a function of the initial state of My, revealing an intuitive pattern of similarities that captures the symmetries of My.
The code is available at ht :/gg thub. www.coo///aaa/e Semm/eeriooalo/SII
8

---

10-2
10a²
Error
m:
Error
m: 1
10
m:
10-3
m: 2
m:
m: 5
m: 25
m:
10
10
10
10
10a
10a
Iterations
Iterations
(a) SVI (Algorithm
SPI (Algorithm)
Figure 1: Estimated transport cost as function k, for various choices of m and 1.
V (80.
X'
y
-1
+1
Figure 2: Visual representation the distances computed between the chains Mx and My.
Discussion
We discuss some further aspects of our framework and results below.
Representation learning for reinforcement learning, Among the numerous applications listed in Section and Appendix A the most interesting for using our metrics for representation learning in RL. As mentioned earlier, bisimulation metrics have been extensively used for this purpose in the past. In particular, almost all such work uses bisimulation metrics compare states within the same MDP and use the resulting similarity metrics for merging states that are at low distance (an approach called "state aggregation"). As our results highlight, this rather narrow view of what bisimulation metrics are capable of: they can define similarity metrics between processes that live on potentially different state spaces, which in particular can be used to select representations by minimizing the distance between high-dimensional process and a set of ow-dimensional representations. Curiously, our LP formulation may allow differentiating the distances with respect to the transition kernels, which we believe will be an important property for future developments in representation learning for RL.
Limitations of the theory. In their current form, theoretical guarantees in Theorems and
only apply to perfect projection and evaluation steps, corresponding to setting m = oo. We conjecture that this limitation
be addressed with more careful analysis, and results similar those of Theorems 2 and can be shown, potentially at the price of a worse dependence on the effective horizon 1
Scherrer et al. 2012 2015), by making of the techniques of Geist et al. 2009 proaaiiirogrrammngg algorithms.
and Moulin and Neu 2023 for analyzing regularized
From dynamic programming to learning from data. This paper focuses on computing distances between known Markov processes via dynamic- pwwarii-rrooraming-yyeeee applications however, the transition kernels are unknown, which requires the development of new tools. We are confident that our framework can serve as solid basis for such developments, and in particular that one can port many ideas from the field of reinforcement learning that is essentially all about turning jynamic-programminng methods into algorithms that can learn from interaction data. Additionally, we believe that our I.P formulation in Section 3 makes it much easier to import further ideas from computational optimal transport, and particular that stochastic optimization methods like those of Genevay et al. 2016 can be adapted to solving our linear programs.
methods. In the most interesting

---

Acknowledgments.
G. Neu would like to thank Marin Ballu and Quentin Berthet for clarifying number of details of their analysis of Mirror Sinkhorn, Tristan Brugère for pointing us to their publicly available code, and Anna Korba and Stefan Schrott for kindly helping out with some references about optimal transport between stochastic processes. G. Neu was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 950180). Anders Jonsson is partially supported by the EU 4CT-48 2020 project TAILOR (No. 952215), AGAUR SGR, and the Spanish grant PID2019-108141C -O0. Javier Segovia-Aguas is supported by MCIN/AEI /tt./11.1.1..// 13039/501 10001 1033 under the Maria de Maeztu Units of Excellence Programme (CEX2021-00) 195-M).
References
A. Abate. Approximation metrics based on probabilistic bisimulations for general state-space Markov processes: survey. Electronic Notes in Theoretical Computer Science, 297:3-25, 2013.
Y.
P. Bartlett, K. Bhatia, N. Lazic, C. Szepesvari, and G. Weisz. POLITEX: Regret bounds for policy iteration using expert prediction. In International Conference Machine Learning (ICML), pages 3692-3702, 2019.
P. Aczel. -ol-welllfounded sets. CSLI lecture notes, 14, 1988.
A. Agarwal, M. Kakade. J. Lee, and G. Mahajan. On the theory of policy gradient methods: Optimality, approximation, and distribution shift. Journal of Machine Learning Research, 22(98): 1-76, 2021a.
R. Agarwal, M. C. Machado, P. S. Castro, and M. G. Bellemare. Contrastive behavioral similarity embeddings for generalization in reinforcement learning. In International Conference on Learning Representations (ICLR), 2021b.
J. Altschuler, J. Niles-Weed, and P. Rigollet. Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration. Advances in Neural Information Processing Systems, 30, 2017.
M. Arjovsky, S. Chintala, and L. Bottou. Wasserstein generative adversarial networks. In International Conference on Machine Learning (ICML), pages 214-223, 2017.
M. G. Azar, V. Gómez, and H. J. Kappen. Dynamic policy programming. Journal of Machine Learning Rexearch, 13(NNv):3207-3222) 2012.
J. Backhoff- Veragua,, M. Beiglbock. Y. Lin, and A. Zalashko. Causal transport in discrete time and applications. SIAM Journal on Optimization, 27(4):2528 22-22222--26,, 2017.
J. Backhoff- Backhooff-Veraguas, D. Bartl, M. Beiglböck, and M. Eder. All adapted topologies are equal. Probability Theory and Related Fields, 178: 127- 2020.
M. Ballu and Q. Berthet. Mirror Sinkhorn: Fast online optimization on transport polytopes. In International Conference Machine Learning (ICML), 2023.
D. Bartl and J. Wiesel. Sensitivity of multiperiod optimization problems in adapted wasserstein distance. arXiv preprint arXiv:2208.0566 2022.
E. Bayraktar and B. Han. Fitted value iteration methods for bicausal optimal transport. arXiv preprint arXiv:2306. 22658, 2023.
A. Beck and M. Teboulle. Mirror descent and nonlinear projected subgradient methods for convex optimization. Operations Research Letters, 3773)1667-175 2003.
R. Bellman. Dynamic Programming. Princeton University Press, Princeton, New Jersey, 1957.
D. Bertsekas. Monotone mappings with application dynamic programming. SIAM Journal on 1503
Control and Optimization, 15(3):438-46) 464, 1977. doi:
10

---

G. Bian and A. Abate. On the relationship between bisimulation and trace equivalence in an approximate probabilistic context. In Foundations of Software Science and Computation Structures (FoSSaCS), pages 321-337, 2017.
T. Brugère, Z. Wan, and Y. Wang. Distances for Markov chains, and their differentiation. In International Conference on Algorithmic Learning Theory (ALT), pages 282-336, 2024.
X.-R. Cao. Single sample path-based optimization of Markov chains, Journal of optimization theory and upplications, 100:527-548, 1999.
P. S. Castro. Scalable methods for computing state similarity in deterministic Markov decision processes. In AAAI Conference on Artificial Intelligence (AAAI), pages 10069-10076, 2020.
Castro, Kastner, P. Panangaden, and M. Rowland. A kernel perspective on behavioural metrics for Markov decision processes. Transactions on Machine Learning Research, 2022.
D. Chen, F. van Breugel, and J. Worrell. On the complexity of computing probabilistic bisimilarity. In Foundations of Software Science and Computation Structure (FoSSaCS), pages 437-451, 2012.
S. Chen, S. Lim, F. Mémoli, Z. Wan, and Y. Wang. Weisfeiler-Lehmann meets Kttoov--Waseesste
In International Conference on Machine Learning (ICML), pages 3371-3416, 2022.
C.-Y. Chuang and S. Jegelka. Tree mover's 's distance: Bridging graph metrics and stability of graph
neural networks. Advances in Neural Information Processing Systems, 35:2944 35:2944-2957,, 2022.
N. Courty, R. Flamary, D. Tuia, and A. Rakotomamoniy. Optimal transport for domain adaptation. IEEE transactions pattern analysis and machine intelligence, 39(9):1853 1993-1853-1885, 2016.
N. Courty, Flamary, and M. Ducoffe. Learning Wasserstein embeddings. International Conference on Learning Representations (ICLR), pages 1-13, 2018.
M. Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in Neural Information Processing Systems, 26, 2013.
G. de Ghellinek. Les problèmes de décisions séquentielles. Cahiers du Centre d' Éttudes de Recherche Opérationnelle, 2:161-179, 1960.
E. V. Denardo. On linear programming in a Markov decision problem. Management Science, 16(5): 281-288, 1970.
d' d'Epenoux. probabilistic production and inventory problem. Management Science, 10(1): 98-108, 1963.
Desharnais, V. Gupta, Jagadeesan, and P. Panangaden. Metrics for labeled Markov systems, In International Conference Concurrency Theory (CONCUR), pages 258-273, 1999.
J. Desharnais, R. Jagadeesan, V. Gupta, and P. Panangaden. The metric analogue of weak bisimulation for probabilistic processes. In IEEE Symposium on Logic in Computer Science (LICS), pages 413 413-422, 2002.
Desharnais, V. Gupta, R. Jagadeesan, and P. Panangaden. Metrics for labelled Markov processes. Theoretical Computer Science, 318(3):323-33 2004.
S. Eckstein and G. Pammer. Computational methods for adapted optimal transport. The Annals of Applied Probability, 34(1A):675-713, 2024.
E. Even-Dar, S. M. Kakade, and Y. Mansour. Online Markov decision processes. Mathematics of Operations Research, 34(3):726-736, 2009.
N. Ferns, P. Panangaden, and D. Precup. Metrics for finite Markov decision processes. In Uncertainty in Artificial Intelligence (UAI), pages 162-169, 2004.
N. Ferns, P. S. Castro, Precup, and P. Panangaden. Methods for computing state similarity in Markov decision processes. In Uncertainty in Artificial Intelligence (UAI), pages 174-181, 2006.
11

---

M. Forti and F. Honsell. Set theory with free construction principles. Annali Scuola Normale Superiore, Pisa, Serie IV 5(3):493-522, 1983.
A. Galichon. Optimal Transport Methods in Economics. Princeton University Press, 2016.
M. Geist, B. Scherrer, and O. Pietquin. A theory of regularized Markov decision processes. In International Conference on Machine Learning (ICML), pages 2160-2169, 2019.
C. Gelada, S. Kumar, J. Buckman, O. Nachum, and M. G. Bellemare. Learning continuous latent space models for representation learning. In International Conference on Machine Learning (ICML), 2019.
A. Genevay, M. Cuturi, G. Peyré, and F. Bach. Stochastic optimization for large-scale optimal transport. Advances in neural information processing systems, 29, 2016.
A. Giacalone, C.-C. Jou, and S. A. Smolka. Algebraic reasoning for probabilistic concurrent systems.
In IFIP Conference on Programming Concepts and Methods, pages 443 443-458, 1990.
R. Givan, T. Dean, and M. Greig. Equivalence notions and model minimization in Markov decision processes. Artificial Intelligence, 147:163-223, 2003.
Hansen- Estrucc,, A. Zhang, A. Nair, P. Yin, and Levine. Bisimulation makes analogies in goal- conditioned reinforeement learning. In International Conference Machine Learning (ICML), pages 8407-8426, 2022.
R. A. Howard. Dynamic programming and Markov processes. John Wiley, 1960.
N. Jiang. Notes on state abstractions, 2018.
N. Jiang. note on loss functions and error compounding model-based reinforcement learning. arXiv preprint ar arXhv:2404.09946, 2024.
B. Jonsson and K. G. Larsen. Specification and refinement of probabilistic processes. In IEEE Symposium on Logic in Computer Science (LICS), pages 167-183, 1991.
S. Kakade. A natural policy gradient. Advances in Neural Information Processing Systems, 14: 1531-1538, 2001.
S. Kakade and J. Langford. Approximately optimal approximate reinforcement learning. In Interna- tional Conference on Machine Learning (ICML., pages 267-274, 2002.
L. V. Kantorovich. On the translocation of masses. In Dokl. Akad. Nauk. USSR (NS), volume 37, pages 199-201, 1942.
M. Kemertas and A. Jepson. Approximate policy iteration with bisimulation metrics. arXiv preprint 220.022202.0288., 2022a.
M. Kemertas and Jepson. Approximate policy iteration with bisimulation metrics. Transactions of Machine Learning Research, 2022b.
S. Kolouri, S. R. Park, M. Thorpe, Slepcev, and G. K. Rohde. Optimal mass transport: signal processing and machine-leaming applications. IEEE Signal Processing Magazine, 34(4):43-59, 2017.
B. Kruskal. Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis. Psychometrika, 29(1):1-27, 1964.
K. G. Larsen and Skou. Bisimulation through probabilistic testing. ACM Symposium on Principles of Programming Languages (POPL), pages 344-34440 144- 1989.
R. Lassalle. Causal Transport Plans and Their Monge-Kantonee Problems. Taylor & Francis, 2018.
T. Lattimore and C. Szepesvári. Bandit algorithms. Cambridge University Press, 2020.
12

---

A. S. Manne. Linear programming and sequential decisions. Management Science, 6(3):259-267, 1960.
R. Milner. Communication and Concurrency. Prentice Hall, 1989.
A. Moulin and G. Neu. Optimistic planning via regularized dynamic programming. International Conference Machine Learning (ICML), 2023.
V. Moulos. Bicausal optimal transport for Markov chains via dynamic programming. In IEEE International Symposium on Information Theory (ISIT), pages 1688-1693, 2021.
A. Nemirovski and D. Yudin. Problem Complexity and Method Efficiency in Optimization. Wiley Interscience, 1983.
G. Neu, Jonsson, and V. Gómez. A unified view of entropy-regularize Markov decision processes. arXiv preprint 1705:0770,,77788, 2017.
K. O'Connor, K. McGoff, and A. B. Nobel. Optimal transport for stationary Markov chains via policy iteration. Journal of Machine Learning Research, 23(1):2175-222 2022.
D. M. R. Park. Concurrency and automata on infinite sequences. In GI Symposium on Theoretical Computer Science, volume 104 of Lecture Notes in Computer Science, pages 167-183. Springer, 1981.
G. Peyré and M. Cuturi. Computational optimal transport. Foundations and Trends in Machine Learning, 1 11(--6)355-667, 2019.
G. C. Pflug. Version-independenne and nested distributions in multistage stochastic optimization. SIAM Journal on Optimization, (40-)1440-14200 2010.
G. C. Pflug and A. Pichler. A distance for multistage stochastic optimization models. SIAM Journal on Optimization, 22(1):1-23, 2012.
M. I.. Puterman. Markov Decision Processes: Discrete Stochastic Dynamic Programming. Wiley- Interscience, April 1994.
D. Sangiorgi. On the origins of bisimulation and coinduction. ACM Transactions on Programming Languages and Systems (TOPLAS), 31(4):1-41, 2009.
B. Scherrer. On the performance bounds of some policy search dynamic programming algorithms. arXiv preprini arXiv: 000.13006.0.00 0539, 2013.
B. Scherrer, V. Gabillon, M. Ghavamzadeh, and M. Geist. Approximate modified policy iteration. In International Conference on Machine Learning (ICML), pages 1207-1214, 2012.
B. Scherrer, M. Ghavamzadeh, V. Gabillon, B. Lesner, and M. Geist. Approximate modified policy iteration and its application to the game of Tetris. Journal of Machine Learning Research, 16: 1629-1676, 2015.
G. Schiebinger, J. Shu, M. Tabaka, B. Cleary, V. Subramanian, A. Solomon, J. Gould, S. Liu, S. Lin, P. Berube, L. Lee, J. Chen, Brumbaug, P. Rigollet, K. Hochedlinger, R. Jaenisch, A. Regev, and E. S. Lander. Optimal-transpor transport analysis of single-cell gene expression identifies developmental trajectories in reprogramming. Cell, 176(4):928-94) 900,-9433330 2019.
Y. Shi, V. De Bortoli, A. Campbell, and A. Doucet. Diffusion Schrödinger bridge matching. Advances in Neural Information Processing Systems, 36, 2024.
R. Sinkhorn and P. Knopp. Concerning nonnegative matrices and doubly stochastic matrices. Pacific Journal of Mathematics, 21(2):343-348, 1967.
Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole. bcore- generative modeling through stochastic differential equations. arXiv preprint arXiv:20/1. 133456 2020.
R. Sutton, D. Precup, and S. Singh. Between MDPs and semi-MDPs: framework for temporal abstraction in reinforcement learning. Artificial intelligence, 112(1-2):181-111 1999.
13

---

V. Titouan, N. Courty, R. Tavenard, and R. Flamary. Optimal transport for structured data with application on graphs. In International Conference on Machine Learning (ICML), pages 6275- 6284, 2019.
J. van Benthem. Modal Logic and Classical Logic. Bibliopolis, 1983.
F. van Breugel and J. Worrell. An algorithm for quantitative verification of probabilistic transition systems. In International Conference on Concurrency Theory (CONCUR), pages 336-350, 2001.
C. Villani. Optimal transport: old and new, volume 338. Springer, 2009.
T. Xu, L. K. Wenliang, M. Munn, and B. Acciaio. COT-GAN: optimal transport. Advances in neural information processing systems, 33:8798 13:8798-8809, 2020.
Generating sequential data via causal
Y. Ye. The simplex and policy-iteration methods are strongly polynomial for the Markov decision problem with fixed discount rate. Mathematics of Operations Research, 36(4):593-603, 2011.
B. Yi, K. O'Connor, K. McGoff, and A. B. Nobel. Alignment and comparison of directed networks via transition couplings of random walks. arXiv preprint arXiv:2106.071 106, 2021.
A. Zhang, R. McAllister, R. Calandra, Gal, and S. Levine, Invariant representations for reinforce- ment learning without reconstruction, In International Conference Learning Representations (ICLR), 2021.
14

---

A Extended discussion of related work
In this appendix we include a discussion of related work that could not be accommodated in the main text due to space limitations.
A.1 Bisimulation
The concept of bisimulation originated independently in modal logic van Benthem 1983
computer science [Park] [981 Milner 1989] and set theory [Forti and Honsell 1983 Aczel 1988 with firm roots in lixed-point theory. Bisimulation was originally devised as a tool for determining whether or not two processes are behaviorally equivalent in the sense that no test can distinguish between the labels they generate, Being much less demanding notion of equivalence than isomorphism (which is generally NP-hard verify), the notion of bisimulation has had significant impact in concurrency theory and formal verification of computer systems, and has become a standard tool for model checking. We refer the reader to the very enjoyable paper of [Sangiorgi 20099 for detailed history of bisimulation and related concepts.
Larsen and Skou [1899] developed theory of probabilistic bisimulation between stochastic processes. Their approach is based on comparing interpretations of logical formulas on stochastic processes, roughly saying that two processes are probabilistically bisimilar if all formulas acting on the sequence of labels encountered along the corresponding random trajectories follow the same probability distribution. Their definition can be most simply presented when the two processes live on the same state space and follow the same transition kernel, but are initialized at two different states. In this setup, bisimulation reduces to relation between individual states, which can be described formally as follows. Given a stationary Markov process Mx (X, Px, vo.x) as defined in Section 2 and label function f : X + IR a probabilistic bisimulation R is relation on X' x X' that satisfies the following property: two states z and x' are bisimilar (denoted x: ERr') if and only if L(x) = [(x')
and for each subset C in the partition X' \\ induced by R. holds that
Σ Px(x"<x)= = - Px(x"1x') EC
x"EC
Jonsson and Larsen 1991] define another similarity notion called "satisfaction relation" based on couplings, and prove that probabilistic bisimilarity and satisfaction are equivalent notions of similarity. These works show that bisimulation is indeed an equivalence relation, and that when two processes are initialized from two states within the same equivalence class (i.e., they are bisimilar), then they will not only transition to bisimilar states in the next step but will in fact continue to evolve in way that is indistinguishable based on the labels (in the sense that they will produce the same distribution over sequences of labes).
Giacalone al. 1990], Desharnais al. 1999 2004] and van Breugel and Worrell [2001] relax the restrictive notion of exact probabilistic bisimulation and introduce real-valued pscudometrics that measure the degree of bisimilarity between two states. These notions rely on real-valued labeling functions. Giacalone et al. [1900] gives a notion of e-bisimulation which relaxes the condition [(x) L(x') in the definition of the hard bisimulation relation given above, and only requires equality to hold up to some > O. Desharnais et al. [1999 2004] go further and define genuinely real-valued extension of bisimulation relations by defining bisimulation metrics as described in the main text (and particularly Equation (4). A fixed-point characterization of these bisimulation metrics was established by van Breugel and Worrell [2001] and Desharnais et al. [2002]. These works show that the resulting distance notion is in fact pseudometric, and two processes are bisimilar if and only if they are at distance zero. This justifies seeing bisimilarity metrics as "soft" extensions of the binary relation of bisimilarity.
Interestingly, the first bisimulation metrics all make use of concepts from optimal transport in some way or another: Desharnais et al. 1999 2002] already note that their definition is inspired by the Wasserstein distance (which they call the Hutchinson metric"), and the fixed- -oint characterization of van Breugel and Worrell [2001) also make use of this distance to compare transition kernels (cf. Equa- tion 5). Precisely, their definition that we recalled as Equation (4) is admittedly inspired by the Kan- torovich dual representation of the Wasserstein distance between probability distributions over metric spaces. To our knowledge, this connection with optimal transport has not been explored further the literature, and in particular no "primal" counterpart based on couplings has been discovered so far.
15

---

Givan et al. 2003 adapt probabilistic bisimulation to Markov decision processes (MDPs), requiring states (or stateactiion pairs) to have identical rewards and transition probabilities for two MDPs be bisimilar. |eenns et al. 2004 2006) introduce bisimulation metrics for MDPs, essentially using the fixed-point characterization of van Breugel and Worrell [200T] as starting point for their definition. Once again, their definition makes use of the Wasserstein distance between the transition kernels (called Kantorovich metric in the paper), but no deeper connection between bisimulation metrics and optimal transport is discussed.
Castro 2020 defines bisimulation metrics for MDPs with respect to a given policy π. which thus falls back to the standard definition of bisimulation metrics for Markov chains as studied in the works
of Desharnaiss et al., [1999] 2004) and van Breugel and Worrell (2001). Kemertas and Jepson [2022b) adjusts the definition of Ferns et al. [0044] [2066] by replacing the Wasserstein distance with the
regropy-regulariied Wasserstein distance proposed by buCuuuri [203)) which allows them to apply the
Sinkhorn algorithm to compute bisimulation metrics for MDPs. The resulting approach can be seen
to be nearly identical to the methods proposed by JyCCCnnnor et al.] [2022] and (Bugère et al. [2024) for approximately solving the fixed-point equations )3) in the context of optimal transport see Section A.2 for further discussion of these works.
Chen et al. 2012 study the computational complexity of computing bisimulation metrics. Similar to our work, the authors formulate a linear program that characterizes bisimulation metrics in an equivalent way to other common definitions, though the linear program has one constraint per next-state coupling, which makes their LP intractable as stated. They use their LP formulation as an analytic tool to show the existence of polynomial-time algorithm to solve the fixed-point equations (5) (which essentially amounts to Bellman's value iteration algorithm implemented in the MDP we describe in Appendix [B. Each step of the resulting algorithm solves one optimal-transpor problem per state pair via the network simplex algorithm, which is known to be impractical for this purpose in comparison with Sinkhorn-style methods [Cuturi] 2013]. In the context of optimal transport, a closely related linear program has been discovered by Backhoff- Backhoff-Veraguas et al. [2077), whose framework is more general in that it mostly focuses on general (potentially nan-Markovian) stochastic processes, but with the limitation that only finite horizons are considered. We discuss further developments on this topic in Section A.2
Finally, Bian and Abate 2017 show that e-bisimilar Markov processes generate distributions over finite-length trajectories that are close in total variation distance. Since this distance is special case of the Wasserstein distance (with the Hamming metric over sequences as ground metric), this result can be seen to establish some relation between OT distances and bisimulation, but the link is rather weak in the sense that no equivalence is shown between the two notions. Indeed, these results only imply that nearly-bisimilar processes generate nearly-identical trajectory distributions, but the reverse implication is not shown to hold.
A.2 Optimal Transport
Optimal transport VVillnni 2009] studies the problem of transporting mass between two density functions and g, given a cost function that measures the transport distance between any pair of points. In the classic formulation of Kantorovich [1942), the vehicle used to transport mass is a coupling, that is, a joint probability distribution whose marginals equal p and The problem of finding an optimal coupling that minimizes the total transport distance can be formulated as linear program. Historically, the resulting LPs have been solved via standard solvers like the network simplex method or interior-point methods, which lead to algorithms with polynomial runtime guarantees but rather poor empirical performance. Cuturri (2013) successfully advocated for adding entropy regularization to the standard LP objective, which enabled algorithms that are orders of magnitude faster than previously proposed methods.
In this work we consider a problem of optimal transport between stochastic processes with a temporal dimension. This topic has recently started to receive attention in the OT literature. mostly focusing on stochastic processes with finite horizon [Pilug and Pichler, 2012. [Bakhhoff- Veraguas et al.) 2017] Lassalle| 2018). In this setting (often called "adaptedd transport", "causal transport", or "bicausal transport"), ) the problem is to transport mass between joint distributions of sequences of elements, which can be formulated as an optimization problem over the set of causal couplings (i.e., the
set of couplings over joint distributions over sequences that respect the temporal order inherent in the process). Backhoff- Veraguas al. [017]] have observed that, due to the linearity of the
16

---

causality constraints, this optimization problem can be phrased as a linear program, which however is infiite-dimensional and thus intractable solve directly. They complement this view by providing principles for characterizing the structure of the optimal coupling, which, in the special case of Markov processes, boils down to the finite-horizon version of the fixed point equations (3). This development essentially mirrors the LP formulation and
principles put forth by Chen et al. [2012] in the context of computing bisimulation metrics (cf. the discussion in Section|A.1]
Still on the front of computing optimal transport distances, notable contribution is due to Eckstein and Pammer [2024], who propose and analyze a version of Sinkhorn's 's algorithm for optimal transport on the space of stochastic processes. When specialized Markov processes, their algorithm can be seen to be very closely related to ours, the technical explanation being that in finite-horize Markov processes the entropy of path distributions that they use as regularization can be seen to be equal to the conditional entropy that our method uses for the same purpose. The resulting algorithm performs iterative Bregman projections via backward recursion over the finite time horizon, with computational steps that are essentially identical to applying our Bellman-Sinkhorn operators. That said, their analysis relies very heavily on the finite-horizon structure of the problem and as such is not applicable in our considerably more challenging infinite-horizon problem setting.
The more recent works of Moulos 202T] O Connor et al.) [2022] Bayraktar and Han [2023 and
Brugère et al. [2024) have investigated optimal- transport distances between infinite-horiie Markoy chains. O'Connor al. [2022] considered the undiscounted version of our problem and proposed to compute optimal transition couplings via an adaptation of approximate policy iteration (cf. Scherrer
2013) to an appropriately adjusted version of the MDP we describe in Appendix B Their key algorithmic idea is approximating the greedy policy update steps by running Sinkhorn's is algorithm
for each pair of states. Essentially the same idea was used by Brugère et al. [2024) to solve the
discounted problem that is the subject of the present paper, with the difference that their method takes approximate value iteration as its starting point. Both of these approaches are closely related to the alternative fixed-point definition of bisimulation metrics using Sinkhorn divergences due to Kemertas and Jepson [022b]] as mentioned in Section While these approaches are nearly as effective as our Sinkhorn Value Iteration method in practice, their black-box use of Sinkhorn' 's algorithm make them difficult to analyze theoretically, and difficult to build further theory on.
To wrap up, let us mention some results that in sense have already foreshadowed our observation about the relation of OT distances and bisimulation metrics. First, we note that Yi et al. (202T), J Brugère et al. [024] proposed to study optimal transport distances of Markov chains defined over graphs as a means of studying the similarity of the underlying graphs. The purpose of these works was to define a notion of distance that is less demanding than isomorphism, but is still grounded in fundamental theory and can be computed effectively- which is precisely the reason that the notion of bisimulation [1981] and Milner 1989]. Finally, the work of Backhoff- Verffuas et al. (020) has established that all adapted topologies are equal" on the space of laws of stochastic processes, understood in the sense that a large number of topologies (including the one induced by optimal-transpon metrics defined in terms of bicausal couplings) in fact identical. While one may argue with their sweeping claim that all such topologies are equal, it may not be surprising in light of their results that the topology induced by bisimulation metrics is also identical to these well-studied topologies (as revealed to be true by our observations in this paper).
originally introduced in the 1980s in the context of formal verification by Park
17

---

B Optimal Transport as a Markov Decision Process
Consider the two Markov processes Mx and My with initial states 2g and Yo. respectively. To compute the optimal transport cost W, (Mx, My; c, 2o. yo) and the optimal transition coupling x², we can introduce a Markov decision process
M = (X A, 4, 7, C, vb) where:
XY is the state space, defined as the set joint states of Mx and My. A(zy) Hyy is the set of applicable actions state zy E XY,, corresponding the set of valid couplings of (x(-(a:) and Py(-lv).
d(-lzy, a) = a is the transition probability distribution, which is fully determined by the action Ila..
y is the discount factor,
c: xy - [0, o0) the cost function that maps joint states to positive real numbers,
1/u 820g0 is the initial state distribution.
The objective of the agent in this MDP is to select its sequence of actions Ao. A in way that minimizes the total discounted cost E 20 7'o(X, Y.)], where each state pair is drawn according to the action taken by the agent as (Xe, Y)) A.. The sequence of actions is generated by sequence of history dependent policies Te E HHD H4 - A((...Y.) where He = (Xo, , Yt). A first remark is that we can restrict ourselves to deterministic policies. Indeed, the action set is convex at every time step and both the reward and the transition probability distributions are linear in the action. A second remark is that bicausal couplings correspond exactly history- iistory-dependent deterministic policies. Indeed, by bicausality, the bicausal coupling Mxy is generated by the policy Mxy(22 (Zn, Yo,
We will denote this policy by TMay
Of special interest are stationary deterministic (or Markovian) policies of the form π XY A, as A, = π (X,, Yt). Such policies correspond exactly with transition couplings as defined the main text as mappings : XY - Axy the same type. Accordingly, will sometimes write т T-|2zy) to refer to the distribution (xy) A(ry) below.
mapping joint states (X,, Yt) to actions in A(X,, ,, Yt)
B.1 Value functions, optimal policies, and sufficiency of transition couplings
Each policy æ induces a value function V*, defined in each state xy E XY as
(T(ry) = Ex
Σ t=
NoYo = ry
(11)
where the expectation taken with respect to the stochastic process induced by the policy т. Here, X, and Y, are random variables representing the state of the two processes time l. In particular, have that for any bicausal coupling Mxy E The.
f (G 1 (Xx,,, (.)M2os88)
aa Σ Eree(Xe Yc)
=Etaxy
XoYo =x000
VTMxy (Io,
And as result can relate the optimal transport cost to the optimal value function of the MDP
/ Y)
W, WM My; 2o, Mn)
inf MxyElb.
inf Mxy
VaMxy (Iog0)
inf V V=(ooy)).
rellin
18

---

Building on classic results of MDP theory, it can be shown that there exists an optimal Markovian policy T" whose value function V = V* satisfies V (zy) < (z(zy) for all policies x and joint states zy, and said optimal value function V satisfies the Bellman optimality equations
(*(zy) = (VV))
where T the Bellman operator acting on a function V KR as
(V = c(ry) +
inf Σ E(('v))(rv')) (Vry). pElley x'y
This is precisely the set of equations in Equation (3). Since V is optimal, V*(yy) equals the optimal transport cost W, (Mx, My: c, z, y) for each state zy. An optimal policy л achieves the infimum in each state xy, with associated value function V" = V* These claims are summarized in the following theorem, stated in nearly identical form by Moulos [202]]
Theorem 3 (cf. Theorem in Moulos ) Under the above conditions, the following hold:
There exists an optimal Markovian transition coupling x such that for any policy и E HDD and any joint states xy, V (ry) S (xy))
The value function of π satisfies W. W.MMx, My;c, V* (xy).
There exists unique solution the Bellman optimality equation denoted V E R.Y...
We have that V V
The above theorem justifies considering Markovian transition couplings when computing the optimal transport cost (M My; c, Fo. yo).
Proof. One can straightforwardly check that our setting is an instance of optimal control prob- lems with additive cost functional studied in [ertsekas] (1977). In particular, the contraction assumption (Assumption C in the paper mentioned above) satisfied. We define pointwise V*(xy) = inf E (xy), and Proposition of Bertsekas [1977) shows that V the unique solution to the Bellman optimality equation.
Now, for the existence of an optimal stationary policy or transition coupling in our terminology, an additional technical condition on the action set must be carefully verified. For every xy E X), f E [0, oo) and k. we define the set
Un A) { a E A(zy) : c(x, u) + 7 (VV(
This set is compact the intersection of the compact set A(zy) with a closed set, the preimage of a closed set by a continuous application. Knowing that Us (xy, A) is compact, we can then apply Proposition 14 of Bertsekas [1977) and that gives us the existence of an optimal Markovian transition coupling x* that satisfies V = V*. In particular, for any policy Ihи and joint state zy, have that Va (xy) = V (ry) < (zy) by definition of V*
L
B.2 Occupancy measures and occupancy couplings
In finite Markov decision process, occupancy measures express the discounted number of times that given state and action are visited on expectation by the controlled stochastic process. This notion is transport problem, given that the
meaningfully applicable in the MDP formulation our action space is infinite. However, the closely related notion of occupancy coupling can be seen to play similar role that it allows expressing the waal--dicounee-ccoo objective as a linear function, and that the set of valid occupancy couplings can be fully characterized in terms of finite number of what follows, we prove this latter key property of occupancy couplings, stated as Lemma Tin in the main text.
linear constraints.
Proof of Lemma 7 We begin by showing that the occupancy coupling n² associated with any transi- tion coupling 1xxy satisfies the following system of equations (sometimes called the "Bellman
19

---

flow equations")
pe (oy..'yy) = π (2'y''lla
(r"""" xzy) + +1--7)u
(12)
8''a""
Indeed, this be shown to follow from the definition of occupancy couplings
4² (xy, d'y) = (1 (1-7)
yI,, = zy. x z'vl
1=
8
-11-7) Σ
7 7 ky [X,Y= ry
<<''y' y'
(1-7)u0(zy) (1 - 7) Σ
7'P [XeYe = ryl
t=1
(1-7)oo(zy) +7 Σ 7) Σ
1. [X 11 = x"y", X,Y, xyy
x"y²
<a'a yoy) (1 (1-7)uo(zy) +7 E 1. (x"""" vy)
x"y"
where in the first step we used the stationarity of the transition coupling т, then the definition of vo. followed by the law of total probability, and finally the stationarity of the Markov chain that allowed us to recognize μ (x"y", rry) in the last step. Now, summing both sides of Equation (12) for all r'y', we can confirm that μ indeed satisfies Equation (6). Furthermore, summing the two sides of Equation (12) over all x', we get
Σ (zy, (zy, 2'1)
Σ M(d' (oy)) (1 (1- + Enr xxy)
x"
Σ mlalvl ('(
E vn (xy,,""""")
x"v" ptt (xy,,"")""
Py(lll)) Σ
x"
where the first step in the second line follows from Equation and the 4th line comes from the fact
that (--lary) is coupling of Px(+-a:: and Py)())). This verifies that n²" satisfies Equation (8) and the same reasoning can be used to verify that it also satisfies Equation 7
Conversely, suppose that μ RRX satisfies Equations Ex'y plzy, x'y/) and let
and
Define Up (pyy))
play,rry' v2) Mylayy) Px (x
if v 00, otherwise
Ty ( (1a)
(x'
We will verify that Tu defines valid Markovian coupling and that is the state action occupancy measure of Ty. If v (wyy) = O, then Ex, nu (x' (2's s/lay) = Py(s/l)) En, Px Px = Py (v). y (zw) 4 0, then we have
Σ ne (x's''a)
Σ E(rry
Σ M(yy
(xy)
Vu (x()
x'x¹'
(,)) Py (ay) = Ps(s'ly)
where we have used that μ satisfies the constraint of Equation (8) In any case, we have that Ex ( ('y) Py(V'v). By symmetry, using that μ satisfies Equation 7 we also have
20

---

Ey u(x'y' my) = Px(x112), Hence we have demonstrated that Tμ is valid Markovian coupling of Px and Py. To proceed, observe that the occupancy coupling associated with Tu satisfies
Σ (xy, a''y'
Σ (x" (a"""""""" y", xr)) + (1 - (oy))
xl'ag"
We will verify that ne μ by showing that this system of equations has unique solution. In order to see this, let us recall the definition of Uu and reorder Equation (6 as
- (->))u(rr)) un (ry) - Σ /" ""mp"""""""" x"y"
Introducing the matrix Z E ( with entries 7 x' y²) Th x (zy) and representing the functions and Nn in matrix form, this system of equations can be written as
(1 - 7)0 (1-
Now, thanks to the Perron-Frobenius theorem, the stochastic matrix has spectral radius 1, and thus - (Z) is invertible, meaning that there is unique solution U (1 (1 - This in turn implies that μ ,,"",, thus verifying that is indeed an occupancy coupling induced by valid transition coupling n if satisfies Equations 66 This concludes the proof.
21

---

C The proof of Theorem 2
The proof is composed of two main parts: showing a bound on the regret of the iterates 11 H2; PKs and then accounting for the enrors incurred when rounding the average iterate Fx = === Pk a feasible occupancy coupling. We start by stating some general results that will be useful throughout the analysis, and then study the two sources of error mentioned above separately. For any j XXX''')'' define Ep((ay) Ex'.. w/ plx'y' V' sry), E" p(pr)) Ex' u(ay,, v²) and
Ty / = play.a'v if E m(rv) + 0, p(z))
(x(('(z))pp Py('l)) otherwise.
In particular, always have that p(xy,, x'y') E p(ay) (a'y'(zy). Furthermore, for any л
X' + Axy, we will denote by क = p(xT) the rounded transition coupling obtained by adapting the rounding procedure of Altschuler al. 2017], and will let p(n) denote the corresponding occupancy coupling μ²² (cf. Section E.2 for the description and analysis of this method). In particular, will always valid transition coupling associated with Px,Py. Finally we define Tout = p (TTT) and Pout ,==",, and we let дe" be an optimal occupancy coupling achieving the minimum in the problem formulation of Theorem
With these notations, we decompose the overall error of the output as
K (ux n², .c)
(Mout 1=,00 cc) (Mont c) + ("x ,,,c) (fomu (K. c) +
K
The first sum the right-hand side corresponds the rounding error, and the second one to the so- called regret of the sequence of iterates μıg- This latter sum can be controlled by adapting arguments from the classic analysis of mirror- 2003). with some ideas adopted from the Mirror Sinkhorn analysis of Ballu and Berthet] [2023] that will also come in handy for analyzing the rounding errors. We state these tools first below, and then analyze the two terms in the above decomposition separately.
methods [Nemirovski and Yudin 1983 [Beck and Teboulle
C.1 General tools
We begin with a version of the classic "three-poim identity" for mirror-descent methods (e.g.. Lemma 4.1 of Beck and Teboulle (003) adapted to our specific setting that involves alternating projections to the sets Bx and By. The result is similar to Lemma A.3 Ballu and Berthet [2023), which we reprove here with a more standard methodology (as used for proving, e.g., Theorem 28.4 of Lattimore and Szepesvári] (2200).
Lemma 2. Let 4 E Bx n By be arbitrary. Then,
(N++) - 1²,c) S H ((x) H (p=())) H (MA (|px)
Proof. We first recall that H. is the Bregman divergence associated with the conditional entropy p(zg.ayv) x'v')
function C(u) = Exy.s'g p(æy,22 a'y)') log Erng" ((zg.a"y") (cf. Appendix A.1 in Neu et al. 2017). For the actual proof, let us consider the case of & odd when uk+11 E Bx. Then, by defimiton, we have that 14++1 is the minimizer of (++i(u)= 77 (n, c) + H (μllиu) on this set. Noticing that the gradient of 4 at μ is written as ( = ne (C VC(ux), the first-order optimality condition over the convex set Bx implies that the following inequality holds for any μ 3x:
(nc + VC(mu+1) VC(uu), # px+1) 2 0.
In particular, using this result for n² (which is indeed in Bx). the claim follows from using the standard three- point identity of Bregman divergences that states
(V((u=-1) vC(u)), n² - pk+1) = H (μ² [[px) - H [[]]))))) - H (px+1llum)
Repeating the same argument for even rounds (and noticing that the comparator also satisfies j²* E By as needed for that case) completes the proof.
L
22

---

The following standard lemma will also be useful for studying various notions of distances between occupancy couplings: the total variation distance || μlll
Eng.rtv' (u(zy, μ (zy, x'y)),
alay.ayy') and the conditional relative the relative entropy (ulll²) = Enyary play, x' y log p plyy,rry entropy introduced earlier.
Lemma 3. For any two occupancy couplings μa and pl, we have
1 1l 2
H (allp')
-'1 s D (ullu() &
1-9
Proof. The first inequality is Pinsker's. For proving the second inequality, we let E and v/ = ETp' be the state occupancies associated with and μ', respectively. Then, we write the following:
(Hll))) D (vlle') + H (plle')
(by the chain rule of the relative entropy)
= D ((1 7)p0 + 71 - 7) + 7Bx') + H (uila')
(using that ys and p' are valid occupancy couplings)
$ (1 7)D (volloo) yD (Ej]) Ep') + H (nll)²)
(using the joint convexity of the relative entropy)
S SDD (wllp') + H (mll)²)
where the final step follows from using the data-processing inequality for the relative entropy. Reordering the terms concludes the proof.
C.2 Constraint violations
Let us begin with some definitions. First of all we introduce quantity that measures the extent to which an occupancy coupling μ violates the transition coherence constraints. Specifically, we will measure the violations of the X'-constraints by
(x(1)= = Σ
v (allry) Eu(ry,avu)
and the violations of the Y'-constraints by
(1) = Σ (.(l)))B//l))>>
aya'
and the overall constraint violations will be written as
8(μ) = (1) + (4)
Note that we have 6x (px) 0 for odd rounds and 6y = 0 for even rounds by definition of the updates. We also define the rounding error associated with an occupancy coupling μ as the average total variation distance between the transition coupling ग. and its rounded counterpart मे = p(Tx)):
(
Σ p(xy, 1 11 (+(+))))) 2300""
The first statement establishes link between the quality of an occupancy coupling and the occupancy coupling obtained after rounding the transition coupling satisfy the constraints.
Lemma 4. For any P. satisfying Equation
(6) we have
(p)
(o(H) - n. c) S |cllo (1 7)
23

---

Proof. Let मμ p(T)) f p(u), 24 = E" and vy = E" Furthermore, let us define the shorthand notation v πT to denote the composition of the state-pair distribution with the transition coupling as (v (T)
< (/(ry))(ayy) Then, we have
j - w|1 = 1107 o Ttp - U/ nulli = 1 TM + u/ th u/ mm Uu o Tuull
11 U/) u + Up (T rz)ll
S 11v7 -M,) # + 1|4 o ( rw)ll1
(using the triangle inequality)
112- - A(4) (())
(using the definition of A)
7 (using = (1 7)20 + F.j)
(14)
S p u|11 + A(p),
where in the last step we used that E is non-expansive with respect to the f-norm. Reordering gives
Aī 1ll1 - A(4)
and putting everything together proves the claim of the lemma.
The next result relates the rounding errors the constraint violations.
Lemma 5. For any Axyxxy, we have
A()) S 8(p)
The proof of this result builds on the error analysis of the rounding procedure of Altschuler al. [2017]. and can be found in Appendix F.2 Finally, the last technical lemma (inspired by Lemma A.4 Ballu and Berthet 2023) bounds the rounding errors in terms of the change rate of the occupancy couplings.
Lemma 6. For any k 2 1,
8(pk) & 2 min((1m - Pk+1ll1 11 - (46-111))
Proof. We study the case where the kth update is By projection. For this proof, it will be convenient to introduce the following notation. We define Ix : R - IR xyxx and (with some abuse of notation), Px xXXXXX
Xx'x the linear operators that respectively act on μ via the assignment (IxM) x²) = Egg plzy, "'ryy z'y') and (Px1) x²) Ex",y" p(ry, x"y"""""))))) " (x((()>>) This allows to write (x( = (Ix - -Px))1l so that we have the expression
8(1x) = (xx) + +)(()) = [|Iии Px
where we have also used (y())) that holds thanks to the fact that is even. Morcover, the (ke + 1)st update is Bx-preeeciion and thus have IxpR+1 Px Hence,
8(1x) I Pxuulli
I Ix + PxpR+1 - Px
S I (JL -M++1)ll + 1P pm)l
S 11 - MK+1ll1.
where have used the fact that both Ixx and Px non-expansions for the fy-norm in the last line, which follows from the Mata-processing inequality for the total variation distance. We then conclude the analysis for the even rounds by replacing p4+1 with P4-1 in the argument above, and repeating the same reasoning for odd rounds completes the overall proof.
Having established these elementary results, we now turn to addressing the main technical hurdle: bounding the cumulative rounding errors.
24

---

Theorem 4. The cumulative constraint violations of the iterates produced by Sinkhorn Value Iteration satisfy
(1 (n" (111)))))
16n) ||cl1 K
8(1k)
2m |cll
7)²
Proof. We start by applying Lemma
6 to show 8(ux) S 2 1 MM+11l2 which reduces our task to We this as follows, for any fixed ar >
bounding Ek= |1
18411ll1-
K
K
ak
M6+1ll
la - Mat1ll1 2
+
2
2a
k=1
(by the inequality arithmetic and geometric means)
K (uu+illAm)
aK
2
k=1
(Pinsker's inequality)
K
H (14++1 [(xx))
+
2
(1 7)
k=
(Lemma
aK
K
H. (1r² [ppp)) - H (1 (14411)
17
K
(c,p' pk+1))
2
(1 (1-7)ac
a 7)
A=1
k=1
(Lemma
aK
H H(" (()"|//1)> H (H* [1xx+1)
Σ
(c, p(ux++))) pm+11)
2
a(1 - 7) (since ()")c) S (p(#x+1).c)
(1 -
aK
H H1r (1111114)
" 1|ello
K
Σ A(Hx+1)
+
2
a(1 7) a(1 - 7)
x=1
(Lemma
< aK
+ H (141
2n ||clloc
Σ 8(1x+1)
K
2
a(1 -7)
a(1 - 7)
(Lemma 5
< aK H (144411144
4n ||c|loc
K
+
Σ km
llue -
2 a(1 -7) ax1 7)
where we have finally used that 6 S 2 mn - M4+1ll1 (Lemma 6). Now, we need to make sure Anjiella I turn this into meaningful result. In particular, setting Sallell (1-7)* guarantees
that a(1-7)²
that the constant in question equals 2. that we can reorder the terms to obtain
2H (n" ((m)
8n 8/lle K (1 7)H ((m)l))
|l,e Met1l11 S aK + a(1 -17) (1 (1- 4n ||clloc
This concludes the proof.
C.3 Regret analysis
In this section, we bound the regret of the iterates produced by Sinkhorn Value Iteration.
25

---

Theorem 5. The regret of the mirror Sinkhorn procedure satisfies
H (u" (11) 2 1|clloo
Pp 1¹, ,,,c) - "
Proof. We first apply Lemma 2 to obtain the bound
H (ullps) - H (ullma+1)) - H (uR+illAx)) + (c, ma)
(#z= 1 .) -
(13)
"
Adding both sides for all k = . K, we get
Σ Uux ,², H (pllm)) H (xllax+1) Ei= H (uR+illaR)
+ (C, -M1)
"
H (Mllpa) + 2|ccllo
This concludes the proof.
C.4 The proof of Theorem 2
The proof follows from applying the above results to bounding the rounding errors and the regret. The first of these is handled as follows:
A
8(Txx)
You AK.C) c) $ ||cll0o
22
1-7
1 1-7
Σ k== 8(1k))
H (n Kn
32n ||c| (1 7)3
s 2 ||cllo
K(1 7))
where the first and second inequalities respectively come from Lemmas and the third one from the convexity of 8, and the last inequality comes from the application of Theorem
The regret then bounded using Theorem
1
K
H (1m1)
2 K
(KK c) =
Σ
(pk 1=, , S
+
||el
K
Kn
Putting both bounds together, we obtain
(Nou -1²,,) = (Mom (xx,c) + (px (c,
27l (1t (11)))))
32n ||c|13 ||ccl
Kn (1 7)3 K
Now, to prove the actual claim of the theorem, we now let be the uniform transition coupling and Hi be the associated occupancy coupling. In this case, the conditional relative entropy can be upper bounded as 7 (f" ((m) < log [2 and we can further bound
(Mout 1²,c)) < 2og X1D
32m ||c ||
||ccllss
+
Kn
(1 7)
K
log |2'
2 ||c|l
< 16 ||clloc V
K(1 -7)3
K
log |2(1|
188 [|]|ccl
K(1 - 7)33
26

---

1 / (1-1)3 log |x|||> and the last line is where the second to last line is obtained by picking a|l=ll
K
obtained by noticing that : Finally, note that
S
log
R(1 7)3
holds whenever K 2 1, | 2, | 2 2 and Y E (0, 1).
log |
Viou (zoz0o (MM,,, My;c, vo) = (Pout - ,,,c) < 18 ||llo
K(1 7))
Now, it can be directly verified that the right-hand side indeed at most whenever K is greater than the expression given in the theorem, thus concluding the proof.
27

---

D Sinkhorn Policy Iteration
We describe here simple alternative Sinkhorn Value Iteration called Sinkhorn Policy Iteration (SPI). After introducing this method heuristically, we provide a formal performance analysis, and finally explain its relation to SVI.
The core concept underlying the definition of SPI is the notion of Q-functions, defined analogously to action-valur functions in an MDP. The Q-function associated with transition coupling π is function Q XxY x xy + IR, with each of its entries defined as
(ay x'yl
E Σ x'x Yi) (Xo, Yo) (z. (X,, (x)
Analogously to the results presented in Section B.I it is possible to show that the Q-function of a given transition coupling π satisfies the Bellman equations
Q (zy, drry) -7 Σ """"")""""""" x""",,
""yy"
and that the Q-value function of the optimal transition coupling Q" = Q" satisfies the Bellman optimality equations
(xy, = c(xy) +
inf Σ pell z"y"
pla" y Q ( v²)".
These are respectively the fixed points of the Bellman operator via
T RXXXY XVXXY defined
f) x' v/ c(ay) +
Σ <a v o' d' x", .")
a"ty"
and the Bellman optimality operator T R.VxXY xyxxY defined via
(T x'v/ c(zy) + +7 inf pell a""""
v(a"""))(
The system of equations TQ' = Q" is essentially as hard to solve (or even harder) than the Bellman optimality equations for V stated earlier. However, one can develop an algorithmic approach toward finding an optimal transition coupling by drawing inspiration from the literature on regularized dynamic programming.
In particular, we develop below analogue of entropy-regularrize policy iteration scheme that
is known under many names in the RI. literature: Natural Policy Gradients by Kakade [200T] and
Agarwal al. [2021a], MIP-Expper by Even-Dar al. [2009) J, Mirror-Deccent Policy Iteration by
Geist et al. [2019]. POLITEX by Abbasi-Yaddorri Yadkori et al. [2019] Policy Mirror Descent by [Agaraal
al. [2021a], and the list goes on. This method can be directly adapted our setting as follows.
Starting from an arbitrary initial transition coupling T1, SPI performs the following sequence of
updates for each is 1,2.....4. K:
Round the transition coupling Ta to म. p(Tx)),
update the Q-function by solving the -rentpoint equation Qs = TRxQU.
* if k is odd, then update the transition coupling
Th (xyyl (1y) exp( - Qu ('yl",
Tht (x' (wa) = Ey (e'y""" (oy) exp( - (xy, Px (z'lo))
else if is even, then update the transition coupling as (x'''' ( exp(= (*+1 lay) Enn Th. (x""" (yy) (xy x" Pv(l'l))
As in the case of SVI, it is easy to verify that these transition couplings satisfy the required marginal constraints. Furthermore, one can verify that the updates defined above exactly correspond to running an instance of Mirror Sinkhorn [Ballu and Berthei 2023) in each state-pair xy with the sequence of cost functions Q (ry, ), Q2(xy,,)) .)... Qx (xy, ). similarly how the entropy-regularized policy iteration methods run an instance of entropic mirror descent in each state. We discuss various aspects of this algorithm below.
28

---

D.1 Practical implementation
Like SVI, this method can be seen as perform- ing online Mirror Sinkhorn updates in each state pair zy with a sequence of cost functions QE. which are computed via solving the linear sys- tem of Bellman equations Q = Qr= The occupancy couplings produced by SPI de- noted by j4 = 1 and the final output of the method is produced by averaging these occupan- cies as Hout = Σズ= juk: and then extracting the transition coupling Tout you Notably, there no need to round this transition coupling since each jūg satisfies all constraints by con- struction, and so does their average pour due to the linearity of the constraints.
Algorithm 2: Sinkhorn Policy Iteration
Input: Px. Py. C, 11, 7 K, m
Initialise: T1 <- Px @ Py;
for k 1..., K 1 do
TE t P(TTa);
Q (T*c) m Q:
T update(T Q): (Equation 10
end
Mout +- K -k= Pki
1
Tout 5- Tpani
V +- evaluate(
Output:
Tout> VTou
Final coupling)
One advantage of SPI over SVI that finding an exact solution for the sxed-point equation Qx- T Q& is easier than computing the fixed points required by SVI, thanks to the fact that this is a linear system of equations. A downside of the method is that it requires to run the rounding procedure after each update, at least for the theoretical guarantees to remain valid. The impact of these steps may however be negligible in practical implementations. Similarly to SVI, the ideal updates of SPI can be approximated by applying the Bellman operator to the Q-functions only small number of times m, with m = do corresponding to the ideal implementation analyzed above. We present a pseudocode for SPI as Algorithm 2
D.2 Convergence guarantees
In what follows, we show the following performance guarantee for SPI.
Theorem 6. Suppose that Sinkhorn Policy Iteration is run for K steps with regularization pa-
1-2
8 log [X| and initialized with the uniform coupling defined for each
rameter 7 = अfo
K
ry, as Ti (1 my) (Mx, [[MMx,M;;ee zo. 30) + if the number of iterations is at least
- Then, for any Zulo E Xy, the output satisfies Va (ZuJo) S
5||c|13 log X1| K
(1 4)
As the analyses of all regularized policy iteration methods listed above, this one also starts with establishing the following claim that corresponds to the classic performance difference lemma (often attributed to Kakade and Langford 2002]] but proposed much earlier in works like Cao 1999 and even Howard 1960). To state the result, we define Ve(xy) = Exy T& (x'x y Quauuy(yyy. x x'y') and
the operators E_ Xyxxy - xy and E_ xVxxy R with each element given by (E_ V) x' z'y') = Vxxy and (E, (z x'y') = V(z'y'). Then, the following bound holds on the instantaneous regret of SPT in round k..
Lemma 7. (m ,=, c) (n", F._ Vv Qu)-
Proof. The proof follows from elementary properties of occupancy couplings. First, note that by the definition of the value function Ve, we have
(T c) (1 (1-7) (po, Vv).
Furthermore, by multiplying both sides of the Bellman equations Qx = c + Vv with μ³" and using the flow constraints _',,²= " μ¹ + (1 - 7)k,, obtain
( Q+) = (14", .. = c E_V)) - (1 (1-)) (1o, Vo)
The result follows after reordering the terms.
29

---

Given the above lemma, we can readily express the regret of SPI as follows:
K
,²,))
(p, E. - -Qu)
k=1
E ( (+(oy) 72 (+ay)) Qu(ry, ))
where we can recognize the regret of Mirror Sinkhorn in each state pair zy E XY. Thus, applying the bound of Theorem 3.1 of Ballu and Berthet [2033 each of these terms (while noting that |0x S ||clloo/ (1 - 7) holds for all k) gives
DEL (1r (+<x0)) || (-128))
9n |lel K 8(1 - 7)2
Σ (x' (-1yy)))) #x(-(zy) Qu .) =
and thus putting the bounds together we obtain
1 K
H. (4 (1(1)11)
9r ||el|
(Pout ,),
PR 1", S
8(1 7)
Now, setting Wi the uniform coupling, we can further upper bound the conditional relative entropy as H (f" mIm) s log |X|l) and after setting v Teell
| Blog
the bound becomes
K
6 |||c| oc
log [*]|
(Mou ,²,c)=
8K
Finally, note that
6 $||||
log |Xll
V (Fog6) (M,, My; c, Fo, vo, vo) =
1--- (Mou - c) I
(1 7)2
8K
Using a crude upper bound 36/8 < verifies the claim of Theorem 6
D.3 Relation to Sinkhorn Value Iteration
While on the surface, SPI may seem only loosely related to SVI, closer connection can be drawn by making the following observations. First, observe that the iunsition-coulling updates exactly match the updates of Sinkhorn Value Iteration, although there is apparent difference in how the Q-functions are defined. To expose the similarity between the two methods better, let us consider an even round in which Ty satisfies the Xmmmrginnal conditions Ey (x'yf |/ar) P(x" (z:). Thus,
introducing the notation Vx x') Ey Tu(''''lay) Qu (ry, ''y'), we can multiply the Bellman equations by Te (x'y) s/ (zg)// (aal (c:) and sum them up to obtain
ne (1'2/ (an) P(x1(z)
(x x²) =
c(xy) + Σ ('
y'
x"
The (1' (ry)
c(zy) +x "(Pa"la'WWaaa''ii
y
(zy) p(xl(x)
exp (-m (1x() + (
log
c(v))
where the approximation is accurate as 7 approaches zero. Thus, for small values 7, the Sinkhorn- Bellman operator used by Sinkhorn Value Iteration is an accurate approximation of the Bellman opera- tor used by Sinkhorn Policy Iteration, and thus one may reasonably expect their respective fixed points to be close as well (this intuition may, however easily fail and is not necessary for our analysis above).
30

---

E
Auxiliary proofs and technical results
In this appendix we prove several technical results from the main text.
E.1 Proof of Proposition
We prove the claim by showing that the solution of the constrained optimization problem of Equa- tion (2) is equivalent the transition-coupling update rule specified in Equation (TO). To this end, study the Lagrangian of the optimization problem (2) corresponding to the update for odd rounds, Bx. and note that the update rule for even rounds, By, be worked out analogously. By introducing Lagrange multipliers Vx (æy, x') for each constraint in Bx, we obtain the Lagrangian
1
C(u; Vx) (e, c) + (mlluss)
"
((15
+
Σ
u(a"'" v,y) + (1 - -7)0((o)) Px(x'(z) - Σ (llry,x'') x""""
Σ x(xy,x'))) ((xy)) +
Σ (xV Vx x0)
2y.2'y yn
Σ vo(zy) (x) (xxx))) - (H|l
+ (1-7)
"
quick calculation (cf. Appendix A.1 of Neu al. 2017) shows that the derivative of H (al))))) satisfies
0HL (pllA)) = log nu (x' //zy)) log nu (d'vll2)) log T (1's///// - log TT (x'v/ //rry), Oplzy, x'y))
where we have used that THe wa holds by definition of μu and mR. To proceed, for a fixed Vx, we set the gradient of the Lagrangian to zero and solve for the transition coupling T6+1, which gives
Tk+1 (x' 5/ly)) Ta(a l/yzy) exp (-v (c(es)) + PPlle"le)Vaaal
x!'
Then, the correct choice of Vx has to be such that the constraint Ey Th+1 (x'y (wy) = (x is satisfied. To see this, suppose that this condition indeed verified and that 1 is the occupancy coupling associated with Th+1- We need to show that M241 indeed verifies the condition defining 8x To this end, notice that μK+1 satisfies Equation (6). and thus the condition can be simply written as
Σ (ay' = Σ (utiloy,,""")) z"y"") (x(x)()),
x""""
plav.z'v)) which, after recalling the relation Th+1(a' (4'3/ (ry) Eomry (( can be indeed seen to hold Ey Thti 1 myy Px (x'lz) is true.
Enforeing this constraint gives the following expression for Ve(wy,
log Σ (2 (-1 (cden) + Pate""l"" Px(x'(a)
Vx(xy,x)) 7)
exp
+7
To conclude, we need to ensure that this system of equations has unique solution. In order to do this, we recall the definition of the Sallmnn-Sinkhor operator T xxxxx + xxxx acting on a function f as
Ta(xy' (ry) exp -17 (c(oo)+ + [e((())))(V))")))
(T ()(''
- -lgg 17
Px (x'lxx
31

---

With this notation, we can directly verify that Vx satisfies TxVVx Vx Vx. Furthermore, it can be shown that the operator xx" is a 7-comtraction in supremum norm (cf. Lemma { and thus it has a unique fixed point by the Banach fixed-point theorem. We finally note that defining Qu(xy, x'y/) = (c(ay) + +7 End Px(x" Px(a"l=")Vxxl Va'yy' ,'"), the transition-coupling update derived above can be rewritten as
Th (zy) exp ( (- Ey The (x' y/ wyy exp (ry, (x
Th+ (zy)
(14)
This concludes the proof.
E.2 Rounding procedure and proof of Lemma
We adapt the rounding procedure stated as Al- gorithm 2 of Aftcchule et al. [2017) and repro-
Algorithm 3: Rounding procedure for couplings
duced here as Algorithm g (where/ / denotes
Input:
approximate coupling F, margins p. !
X t- diag(min(p//C (F)))// 11,
element-wise division in the pseudocode). For- mally, for two probability distributions p E A(X) and E A()),, the set of valid cou- plings is Up = {P RXY, P = p; 1 +). For nonnegative matrix F e R. the rounding procedure outputs a valid coupling p(F, p.q) E Upgg which, by Lemma 7 of Altschhere al. [2017), satisfies
F' <- XF;
Y diag (miin(/i//(FF)) 1))1)):
F" t FYY;
e = F" 1, e 9 F"T 1;
Output: G t F" + emr,er,, ||erpp
|p(FF - Fll, s 2 (I 1-1 pll, + +11 1 -1- =111)
We will now define the rounding procedure for a (not necessarily valid) transition coupling æ XXVXXX by using the aforementioned rounding procedure at each state pair as
-(- = p( Px((1))
With some abuse of notation, we will write the resulting transition coupling as मे = p(T) and the associated occupancy coupling as ā n² p(μ).. Because of the correctness of the original rounding procedure, this transition coupling is valid, and so is the associated occupaney coupling. We can now proceed to the proof of Lemma
Proof of Lemma
Let μa E TR-X'xy>' and as before define vj (yy() = Ex'r' ualay, and
play.rry y') if (xyy) + 0, My (zg) Px (x) Py (w) otherwise.
T vf (zy)
For arbitrary state pairs xy, we use Lemma 7 of Altschuler et al. (017) to obtain that
Py(w'ly) - //
(+1r)) (-v(((11))11 2 Px((dl())= - Σ mnlad'ylaan
Σ
Σ y
+
Now, multiplying by (zy) and summing over zy, we get
Σ P/ (xv) 1| (+(vvl)
< 2
Σ un(zy) Px(a'(z) - Σ
+ Σ
(x(zy)
Py(o'ly)
Σ Ta (222////r))
aype!
11'
xyye'
Σ ayze
P (a(xy) Px(1(x))
Σ ((
Σ zuy
(3(
> Euloy,st')))
(p) + (y(p)= =5(m):
where we used the fact that uara,, x' "/²) = 11 (y(ly) T (z'y/ (zy) for any xy. a'y/ and the definitions 8x
by and 5. Finally, we notice that A1μ) = Exy up (zy) (+10)))) (11l1))) which concludes
the proof.
32

---

E.3 The contraction property of the Bellman-Sinkhowd operator
Lemma Let π : Xy Axy be arbitrary and consider the associated Gellman-Sinkhoon operator T acting on function f XY x y + as
exp (-1 (do))+ +1 Σ I(d""MMdd
(T"') = log
p))
22"
Then, T is a y-contraction for the supremum norm 11-1 that is, for any two functions f1. f2
XY x X, we have
| f1 - (f2)ll I > fell..
Proof. The claim easily follows from using the standard fact that the function 9p(2) = log En p(y")(( (y') 1-smooth with respect the supremum norm, so that for any vec- tors 2, we have (gp(z) - 9p(2]] < || - 2'1loo- To apply this result, we define q(3'(zy) = and (ry, x'y y/) c(ay) FrE..
(x'v) //llzy//((/'l c(zy) + Ez P(x" (2) x"), so that we can write
x") and 22(xy,
IT f T f2llo max (=1
Io(-1 (22(1ry
zy.r²
S max [ل - 22 x' )1 = =1|21 - 22llao S 7 llfi f2lloo
zy.
where the last step follows from the straightforward calculation
sup P(2M f2 llfi -flloo ay.e's y xt
This concludes the proof.
33

---

F
Additional experimental results
In this appendix we present the results of additional experiments not included in the main text.
E.1 Impact of the regularization parameter 7)
Besides the parameter m that we have already studied experimentally in Section 5 the only tuning parameter of SVI and SPI is the regularization parameter 1), which takes the role of learning rate. In this experiment, we study two random walks run on two separate 4-room environments [Sutton et al. 1999]. with two separate reward functions fx and ry that together define the ground cost function d/ y) |xx(a ry)l for each state pair. The results of this study for K = 2 10 iterations are shown Figure 3 The error computed as the difference between the distance estimate produced by the algorithm and near-optimal distance obtained by running Algorithm [T for very small value of 7 and a large number of iterations. We observe that higher values of rj lead to faster error reduction in the initial steps, but eventually prevent convergence to the true solution. In contrast, choosing smaller learning rates enables convergence to better solutions, at the cost of making slower progress initially, An intuitive explanation for this is that for larger values of 7). the iterates converge rapidly to the broad proximity of an optimal solution, but then reach a cycle where they continue to perform large updates to the transition coupling which results in large constraint violations, which necessitate large updates, and so on. This is formally verified by the fact that our bound on the constraint violation terms in Theorem 4 are increasing for large values of 1y. Notably, employing time-dependent learning-rate schedule (inspired by Ballu and Berthet 22033) with 7) V leads the best performance. This strategy leverages faster convergence initially and, for sufficiently large number of iterations, also achieves near-optimal solutions.
10
10
10
10-2
Error
Error
n: 10/KK
n): 10/VK
10
7: 0.1
10-3
0.1
7:
0:
n: 10
10
10
10
10a
10a
102
10a
10a
10"
10a
10a
10a
104
Iterations
Iterations
(a) Algorithm
(b) Algorithm
Figure 3: Error of estimated transport cost as a function k for various choices of
F22
Comparison with alternative methods
We now turn to studying the computational complexity of our algorithms, and compare them empiri- cally to some existing methods that have been proposed for computing optimal transport distances and bisimulation metrics between Markov chains. Specifically, we will focus here on two previous methods: the method of [onCono et al. 2022] that refer
"EntripicOTC" and Algorithm of [Brugère et al.] 2024 that we call "dWL.. We adapt both of these methods with some minor changes to our setting. First, we remove one scaling factor from the transport cost in the definition of the distance defined by |Brugère et al. [2024] so that it matches ours. Second, the algorithm of O'Connor al. 2022] is originally defined for the infinite-horizon average-cost case, and thus we made appropriate changes to adapt it to the discounted case by replacing their approximate policy evaluation step by T applications of the discounted Bellman evaluation operator. As pointed out in Appendix A.2, the resulting methods are closely related, and can be regarded as approximate dynamic programming methods for solving the MDP formulation of our optimal transport problem presented in Appendix B We also recall that the algorithm proposed by Kemertas and Jepson 2022al for the purpose of computing bisimulation metrics also falls into the same class approximate dynamic programming methods, and nearly matches the method of [Brugère et al. 2024 The comparison
34

---

V 0.75
y= 0.9
y = 0.99
1000.0
100.0
SPI
Time (s)
10.0
d0IC
1.0
0.01
10
15 Markov chain size
5
25
10 15
25
Figure 4: Comparison of the computational time of the different methods proposed to obtain a near-optimal solution for different values of 7. o. For each Markov chain size. the results obtained in 5 randomly generated instances are compared, showing the standard deviation in the plot. Data is displayed on log-log scale.
below is based on the original Python implementatior of Brugère et al. [2024) and our own Python adaptation of the MATLAB code of O'Connor al. [022].
One difficulty that we had to face in these experiments is having to tune various hyperparameters of each method (such number of iterations and regularization parameter), which can each influence the quality of the solution and the computation time. For a fair comparison between the methods, have adopted the following procedure obtain our results. First, we estimate ground truth obtained by running one of the algorithms for a very high number of iterations and very low level of regularization, and then use this ground truth the algorithms so that they are close to this value in as little wall-clock time as possible. While all the algorithms perform similar operations, their total runtimes turn out be rather different and heavily dependent on problem parameters such as the discount factor r. The comparison between the resulting runtimes of each method shown on Figure 4 as function of the size the Markov chains, and for three different choices of the discount factor, for a set of randomly generated MDPs (following the setup described Section 7.1 O'Connor et al. 2022).
comparator to adjust the hyperparameters of
First, we observe that EntropicOTC is a policy-iteration-lik method, and as such it needs fewer iterations than the rest of the algorithms we tested to converge to the true cost, but each iteration requires rumming an expensive policy evaluation subroutine until convergence. This computational cost eventually adds up in a way that this algorithm has always ended up being the slowest among all that we have tested, although its performance has proved notably robust to changes in the discount factor 7. Second, we note that the updates of dWL are much cheaper to compute, especially for large regularization parameters. However, the errors of this value-iteration-lie method compound much more rapidly than in the case of EntropicOTC, which makes especially hard to tune the hyperpa- rameters of this method. This problem is especially pronounced when the discount factor is large, which is the most interesting regime as leads to distances with much stronger discrimination power.
The plots shown on Figure indicate that our methods find optimal couplings consistently more efficiently than the other methods in the regime we studied, leading to up 10 times faster runtimes. For the case of sufficiently small values of y, the algorithm of Brugere et al. [2024] sometimes performs competitively, but the hyperparameters leading good performance are much harder find than in the case of our methods. In our experience, the massive speedup achieved by our methods can be largely ascribed to maintaining the transition couplings Tं between iterations, as opposed to computing these afresh by running Sinkhorn's 's algorithm from scratch for each update as done by all other competing methods. Adjusting these other algorithms by maintaining the couplings in memory and using them to warm-start the subsequent updates makes them competitive with our methods, and in fact doing so makes them quite similar to SVI and SPI. Our algorithms use such obscure implementation detail, which is
warm-starts as primary design choice opposed ultimately responsible for the computational efficiency of all these dynamii--rogrammiinn programming methods.
https://ge thub cot.oconn//c kevin/OTC
35

---

The computational time per iteration of each of these methods grows roughly at rate of n², with being the number of states of the Markov chains. In the case our methods, this is easily explained by noting that applying the Sallman-Sinkhone operators and updating the transition couplings requires | operations total. This matches the runtime necessary for running the regularized policy improvement subroutines employed by O'CConnor et al. [2022] and [Brugère et al. [2024], which consists of running an instance Sinkhorn's algorithm in each pair of states. The computational cost of all these methods can be improved by leveraging the sparsity of transition kernels: in particular, if at most S states are reachable with positive probability in both of the chains, the complexity of the updates be trivially improved to | We did not pursue this direction in experiments as our goal was to compare the basic versions of each studied method, and we believe that our conclusions would not be altered if we were to implement this improvement for all methods.
E3 Optimal transport distances similarity metrics
We finally provide a range of experiments that illustrate how the optimal-transpe distances we studied in this paper can be used to capture relationships and symmetries in groups of varying Markov chains. To this end, we have generated 35 different "4-room" instances and computed their pairwise distances. Each instance differs in its initial state and position of the obstacles (amounting to changes in the transition kernel), while maintaining a fixed reward function, with one reward located in each room except for the upper left room. Crossing door between each room results in negative reward. For each instance, we have studied the Markov chain induced by the corresponding optimal policy (which amounts to taking the shortest path toward the closest positive reward, modulo the additional randomness inherent to the transitions). Figure 5 presents a 3D visualization (where the z-dimension represented by color gradient) generated using Multidimensional Scaling (MDS) (Kruskal) 1964] based on the computed pairwise distances. One can observe a clear clustered structure, where instances with similar behaviors are grouped closely together. As the figure highlights, the resulting metrics capture the intuitive similarities and symmetries between each process, which indicates the potential usefulness of optimal transport distances and bisimulation metrics for comparing Markov chains under minimal structural assumptions made on the state spaces and the transition functions.
30
30
70
0-
-10
-10
-20
-30
-30
-20
10
Figure 5: Result of applying MDS the pairwise distances between the of 4-room instances studied. On the plot, the first two coordinates of the MDS embedding are used as the spatial coordinates, and the third coordinate is encoded via the color bar provided on the left-hand side of the axes. It can be observed how the elements in the same cluster present common features that differentiate them from those in another cluster. In the examples shown the figure we can see how the instances in which the closest reward involves crossing a door are concentrated in one cluster, while the instances in which the reward and the initial state are located the same room belong to different cluster. The remaining clusters correspond to having to cross two doors for reward (set of green points on the top), or having no reward that accessible from the initial state (set of blue points in the middle, with large negative z-coordinates).
36

---

NeurIPS Paper Checklist
Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
is fine include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
37

---

Justification:
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architccture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reprodueibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
38

---

Answer: [Yes]
Justification:
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co
public/guides/ wtt.:e/eeeeoeeSSueeeiioooooooooo
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https:
CodeSudmeemii www.iuuubiisseeee.eeeee ssi onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Question:
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [NA]
Justification: There is randomness in any of the experiments that could be meaningfully assessed via error bars.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the enror bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
39

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code Of Ethics
Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics http:pp://errrpss.c/////////cchc
Question:
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: There is no societal impact of the work performed.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no socictal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
40

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Lieenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the paper does existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. wtt/dat has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
41

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: There new assets introduced in the paper.
Guidelines:
The answer NA means that the paper does not release new assets,
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: the paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
42


Reviews:
[{"review_id": "OQMnLUPvYf", "reviewer": "Reviewer_aiTP", "summary": "This work integrates optimal transport with Markov chains by proposing an alternative joint distribution between Markov processes, namely \"discounted occupancy couplings\". They show that optimal transport distances can be computed as a a linear program (LP) in reduced space. This improves the computational efficiency of OT between markov decision process (MDP) over the previous methods. The previous methods, as they reviewed, often requires complex dynamic programming algorithms. They showed that the new formulation can be extended to the well-known entropy regularization problem, employing Sinkhorn-like iterations to solve, making the computation scalable with large problem. This paper provides both theoretical and experimental supports for the new formulation. In the end, they discuss the potential applications on RL models, limitations and challenges.", "questions": "- What is the best practice of bisimulation metric computation? Could your provide the comparison with your methods?\n- What does the entropy means in the sense of the coupling between MDPs? Some explanations on this will be appreciated \n- How does this metric works with data contains noise as we know OT is sensitive to outliers?", "limitations": "None", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The main contribution of this work is demonstration of optimal transport shares the same formulation as probabilistic bisimulation metrics, which is popular in the practice of reinforcement learning (RL). By establishing that solving for optimal transport is equivalent to computing bisimulation metrics, the authors creates a novel link between RL and optimal transport theory.\n- The formulation of OT on MDP is novel. This work provides an new prospective of Bisimulation metric between MDPs.\n- Theoretical analysis looks good.\n- Two computational solutions (Sinkhorn-like iterations) were proposed and tested in experiments sections.", "weaknesses": "The theoretical guarantees are based on the assumption of perfect projection steps $m=\\infty$, which is not practical. As stated in line 290, such exact computation is unnecessary in practice. This claim needs a further theoretical support.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive evaluation of our work, as well as your insightful remarks! We respond to your questions below.\n\nRe weakness: We agree that the analysis working only for the case $m=\\infty$ is the biggest limitation of our results. By analogy with modified policy iteration (that bridges the cases $m=1$ and $m=\\infty$ for standard dynamic programming), we expect that our results should generalize to finite values of $m$, but at the moment we do not have a proof. We nevertheless find the empirical results shown on Figure 1 to be encouraging. In any case, we hope that we were clear enough in stating this limitation in the paper, and hope that you agree that our results are still worthy of publication despite this limitation.\n\nQ1) What is the best practice of bisimulation metric computation? Could your provide the comparison with your methods?\n\nThe best method we are aware of for computing bisimulation metrics is due to Kemertas and Jepson (2022), which essentially aims to approximately solve Eq. (5) by using Sinkhorn's method to approximate the infimum on the right-hand side. Their method is essentially the same as the algorithm of Brugere et al. (2024) that we have discussed in some detail in the appendix. All similar methods are significantly slower in practice than our method. We will add further discussion about these matters to the final version of the paper.\n\nQ2) What does the entropy means in the sense of the coupling between MDPs? Some explanations on this will be appreciated.\n\nThe entropy we use is the Bregman divergence induced by the conditional entropy of $X',Y'$ given $X,Y$ if the tuple $XY,X'Y'$ is drawn from the occupancy measure. The conditional entropy (sometimes called the \"causal entropy\") plays a key role in the theory of entropy-regularized Markov decision processes, and can be shown to be the \"correct\" notion of entropy-regularization in that it induces the well-known \"soft Bellman equations\" as its dual (Neu, Jonsson, Gomez, 2017). We will add a few explanatory lines to the final version of the paper.\n\nQ3) How does this metric works with data contains noise as we know OT is sensitive to outliers?\n\nThis is a great question! So far we have not thought about this, as we have focused on computing distances between perfectly known Markov processes in this work. It remains to be seen how to compute distances between Markov chains based purely on data (that may include noise or outliers). In any case, we are sure that developing algorithms for this more realistic case will require the foundations laid down in the present paper.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "z0BiXWEX5M", "reviewer": "Reviewer_VaVz", "summary": "This submission \n\nNamely, the authors define a notion of optimal transport distance between Markov chains on state spaces with a ground metric. This notion of distance differs from standard optimal transport, as the set of couplings is restricted to the set of so-called bicausal couplings. Using the results of Moulos [2021], it is shown that this Markovian optimal transport can be characterized in terms of the solution of the Bellman optimality equations for a Markov decision process. \n\nUsing this equivalence, it is demonstrated that, in certain cases, bisimulation metrics can be though of, equivalently, as Markovian optimal transport with a specific cost function by exploiting their connection to solutions of certain fixed-point equations. \n\nNext, it is shown that the optimal value for the derived Markovian optimal transport problem can be computed by solving a finite-dimensional linear program, where the constraint set consists of the intersection of three sets. Rather than directly solving this linear program, the authors propose to regularize the problem using a conditional entropy (by analogy with entropic regularization of optimal transport distances which can then be solved efficiently using Sinkhorn iterations). To solve the entropy regularized problem, the Sinkhorn Value Iteration algorithm is proposed and the number of steps required to obtain a desired accuracy in estimating the Markovian optimal transport distance is provided. The authors also propose an analyze an alternative algorithm which they dub Sinkhorn Policy Iteration. \n\nThe paper concludes with some numerical experiments to illustrate the performance of the proposed algorithms.", "questions": "Apart from the points mentioned above, it would be helpful for the authors to state the complexity of the Sinkhorn Value Iteration algorithm in the main text (rather than only in the appendix).", "limitations": "The authors address the limitations of their work in the discussion.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "In my opinion, the submission is well-written and its contributions relative to the broader literature are clearly identified. \n\nThe main contributions of this work, identifying bisimulation metrics as a type of optimal transport problem and providing some new algorithms for estimating bisimulation metrics by using this connection is of interest and is, to my knowledge, novel.", "weaknesses": "1. It appears that the connection between the Markovian optimal transport problem and the finite-dimensional linear program provided in Theorem 1 only enables the computation of the optimal value for optimal transport problem, but does not allow the recovery of the optimal bicausal coupling. If this is the case, this should be further clarified in the text. \n\n2. While the analysis for the Sinkhorn Value Iteration is a nice addition, the guarantees are a bit confusing. Notably, if we wish to attain a precision of $\\epsilon$ in estimating the Markovian optimal transport one requires $K=O(1/\\epsilon^2)$ iterations and setting the regularization parameter to $\\eta =C/\\sqrt K=O(\\epsilon)$. However, if epsilon is sufficiently small, the objective in the update (10) will be dominated by the entropy which is minimized by $\\mu=\\mu_k$. Naively I would assume that $1/\\eta$ should be small so that minimization of the linear term dominates.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive comments and your critical reading of our work! Regarding your questions:\n\n1. Note that the argument $\\mu^*$ achieving the infimum exists, and an optimal transition coupling $\\pi^*$ can be decoded from it. Concretely, given the joint distriibution $\\mu^*$ over $\\mathcal{X}\\mathcal{Y}\\times\\mathcal{X}\\mathcal{Y}$, the corresponding transition coupling can be extracted as $\\pi^*(x'y'|xy) \\propto \\mu^*(xy,x'y')$ (see also line 776 in the appendix for the general formula relating occupancy couplings to transition couplings). We will clarify this in the final version. In any case, we note that our algorithms (SVI and SPI) do return transition couplings besides the transport distances.\n\n2. This observation is entirely correct: the individual updates do get smaller and smaller as $epsilon$ goes to zero. Note however that, at the same time, the number of iterations $K$ goes to infinity at a rate of $1/\\epsilon^2$, and as such the algorithm will make more and more updates that get smaller and smaller as the desired precision approaches zero. This is normal for algorithms that use a fixed learning rate. Following Ballu & Berthet (2023), it would be possible to extend our analysis to time-varying learning rates that start out large and then decay to zero over time, and we have observed in our experiments that such learning rates are indeed easier to tune and work with. We will add a more prominent comment about this in the final version.\n\nRe complexity of SVI: This is a great point, we will update the paper accordingly for the final version!", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "aHhA0f8hfN", "reviewer": "Reviewer_NKJ1", "summary": "The goal of this work is to study optimal transport (OT) distances between pairs of finite Markov chains, providing a novel relation between OT distances and probabilistic bisimulation metrics. The proposed linear program builds on ideas from optimal control in Markov decision processes, and the designed algorithm for solving the proposed linear program combines Sinkhorn's algorithm with an entropy-regularized version of the classic Value Iteration algorithm. Convergence guarantees and computational complexity analysis of the method are provided in the main text and in the appendix.", "questions": "- In the abstract, the authors say \"In this work, we develop an alternative perspective by considering couplings between a “flattened” version of the joint distributions that we call discounted occupancy couplings [...]\" Could you expand on what do you mean by \"flatten\"? \n\n- In the introduction, the authors say that \"The possibility that the objects in question may be random further complicates the picture, and in such cases it becomes more natural to measure distances between the underlying joint probability distributions.\" Could you explain/expand?\n\n- Section 2.1: \n\nDoes $\\gamma$ depend on $t$? \n\nIs $\\gamma^t\\in (0,1)$?\n\n- Line 99, if Markov chains are considered, isn't it that $M_{\\mathcal{X}}$ $(x_n | \\bar{x}_{n-1})$ $=$ $M_{\\mathcal{X}}$ $(x_n|x_{n-1})$ ?\n\n- Line 148 and equation (4): Which is the dependence of $d_\\gamma$ on the labeling function $r$?\n \nLine 149: Which is the dependence of the set $\\mathcal{F}_\\gamma$ on $\\gamma$?\n\nIt is clear from eq. (5) that $U^*$ does explicitly depend on $r$ and $\\gamma$. Could you expand on the equality $U^*=d_\\gamma$? \n\n- Explain the pevious to last equality in (6). How do you make the variables $x',y'$ to apprear?\n\n- Explain the last equality in (6). Why do you have $\\langle \\mu^\\pi,c\\rangle$ if c depends on $(x,y)$ and $\\mu^\\pi$ depends on $x,y,x',y'$?\n\nComments:\n\n- The term \"discounted occupancy couplings\" mentioned in the abtract is not excatly used in the main text, the authors mainly use just \"occupancy couplings\". \n\n- $n$ is used for the size $|\\mathcal{Y}|$ (line 87) and as an index. \n\nIsn't it $|\\mathcal{X}|=m=|\\mathcal{Y}|$? In the experiments, why the choice of $n$ is not analyzed too?\n\nTypos:\n\n- Abstract: line 17, \"method that we call\", instead of \"method we call\".\n\n- Notations: line 79, \"the corresponding finite subsequences\", instead of \"the corresponding subsequences\"\n\n- Line 39 (Introduction): \"important\" is repeated twice. \n\n- Notations: Add the definitions of $\\mathcal{X}^\\infty$, $\\mathcal{X}^n$, $\\Delta_{\\mathcal{X}}$.", "limitations": "N/A", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well-organized, the topic of study is very interesting, and the contribution is novel. Previous works on OT and bisimulation metrics are well-referenced and up-to-date.", "weaknesses": "The main theoretical contributions are presented in Section 3, where the authors introduce the definition of occupancy couplings $\\mu$. These are distributions over $\\mathcal{XY}\\times \\mathcal{XY}$. It is not entirely clear why they need to duplicate the variables $(x,y)$ and consider quadrupes $(xy,x'y')$.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive evaluation of our work and your very detailed reading! We will address your main points below, and will take the remaining minor comments into account when working on the final version of the paper.\n\nQ1: It is not entirely clear why they need to duplicate the variables $(xy)$ and consider quadruples $(xy,x'y')$.\n\nIndeed, defining occupancies as a function of $xy$ is perfectly sufficient for rewriting the optimization objective as a linear function. Duplication of the variables is necessary however for defining the constraints of Eq. (7-9): indeed, without the second set of variables is necessary for stating both the flow constraints and the transition coherence constraints. Consequently, without these additional variables, it would not be possible to characterize the set of occupancy measures via a set of linear constraints.\n\nQ2: Could you expand on what do you mean by \"flatten\"? \n\nWe mean that occupancy couplings are low-dimensional projections (\"flat\" representations) of the joint distribution of the infinite sequence of state pairs. Our main result is showing that this representation is sufficient as long as one is concerned with calculating OT distances. We will clarify in the final version.\n\nQ3: ... Could you explain/expand?\n\nThis sentence was probably not as nicely phrased as it should have been. We will refine it for the final version.\n\nQ4: Does gamma depend on t? \n\nNo, it doesn't (as is normally the case in the theory of discounted MDPs).\n\nQ5: Is gamma^t \\in (0,1)?\n\nYes; we will state this more prominently.\n\nQ6: Line 99, if Markov chains are considered, isn't it that $(x_n | \\bar{x}{n-1})=M{\\mathcal{X}}(x_n|x_{n-1})$ ?\n\nTrue; this line merely served to lay down notation for general stochastic processes. We will clarify.\n\nQ7: Line 148 and equation (4): Which is the dependence of $d_\\gamma$ on the labeling function $r$?\n\nThe labeling function impacts the choice of the function class F. Thanks for pointing this out!\n\nQ8: Line 149: Which is the dependence of the set  $\\mathcal{F}_\\gamma$  on $\\gamma$?\n\nAgain, $\\gamma$ appears in the concrete definition of the function class; we omitted this detail to maintain readability of the main text, but you are perfectly right that we should have at least hinted at this here. We will clarify.\n\nQ9: It is clear from eq. (5)...\n\nSee our response to the previous two questions --- we will make this more clear in the final version.\n\nQ10: Explain the pevious to last equality in (6). How do you make the variables $(x'y')$ to apprear?\n\nWe use the simple fact that $c(xy)$ does not depend on $x'y'$, and as such summing out $x'y'$ in the definition of $\\mu(xy,x'y')$ results in the same discounted sum of indicators as what appears in the last line of Eq.6. We will add a remark.\n\nQ11: The term \"discounted occupancy couplings\" mentioned in the abtract is not excatly used in the main text, the authors mainly use just \"occupancy couplings\".\n\nThank you for pointing this out, we will smooth out the terminology!\n\nQ12: $m$ is used for the size (line 87) and as an index. \n\nThanks for calling this to our attention! This is an unfortunate clash of notation that we didn't notice before and will fix it. In the experiments, we have used $m$ to denote the number of applications of the Bellman--Sinkhorn operators (see Alg 1).", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: Reviewers are unanimous that this is an interesting submission and that it should be accepted.
Decision: Accept (poster)