Bisimulation Metrics are Optimal Transport Distances, and Can be Computed Efficiently

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.

Q1: It is not entirely clear why they need to duplicate the variables $(xy)$ and consider quadruples $(xy,x'y')$.

Indeed, 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.

Q2: Could you expand on what do you mean by "flatten"?

We 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.

Q3: ... Could you explain/expand?

This sentence was probably not as nicely phrased as it should have been. We will refine it for the final version.

Q4: Does gamma depend on t?

No, it doesn't (as is normally the case in the theory of discounted MDPs).

Q5: Is gamma^t \in (0,1)?

Yes; we will state this more prominently.

Q6: 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})$ ?

True; this line merely served to lay down notation for general stochastic processes. We will clarify.

Q7: Line 148 and equation (4): Which is the dependence of $d_\gamma$ on the labeling function $r$?

The labeling function impacts the choice of the function class F. Thanks for pointing this out!

Q8: Line 149: Which is the dependence of the set $\mathcal{F}_\gamma$ on $\gamma$?

Again, $\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.

Q9: It is clear from eq. (5)...

See our response to the previous two questions --- we will make this more clear in the final version.

Q10: Explain the pevious to last equality in (6). How do you make the variables $(x'y')$ to apprear?

We 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.

Q11: The term "discounted occupancy couplings" mentioned in the abtract is not excatly used in the main text, the authors mainly use just "occupancy couplings".

Thank you for pointing this out, we will smooth out the terminology!

Q12: $m$ is used for the size (line 87) and as an index.

Thanks 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).

Thank you for your positive evaluation of our work, as well as your insightful remarks! We respond to your questions below.

Re 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.

Q1) What is the best practice of bisimulation metric computation? Could your provide the comparison with your methods?

The 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.

Q2) What does the entropy means in the sense of the coupling between MDPs? Some explanations on this will be appreciated.

The 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.

Q3) How does this metric works with data contains noise as we know OT is sensitive to outliers?

This 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.

Thank you for your positive comments and your critical reading of our work! Regarding your questions:

1. 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.

2. 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.

Re complexity of SVI: This is a great point, we will update the paper accordingly for the final version!