We agree with the reviewer that these are important points and address them in the discussion after Theorem 3.3 in the current version of the manuscript (lines 235-244). There, we argue for the form of Theorem 3.3 rather than modifying Algorithm 1 because Theorem 3.3 also shows there exists a sequence of step sizes such that $(\lambda\_k,\nu\_k)$ are bounded. Additionally, none of our experiments enforce bounds (they implement PD-LMC exactly as described in Algorithm 1). Yet, we always observe the dual variables converging (see, e.g., Fig. 2 in the manuscript), suggesting that pathological cases are not common. That being said, the theory would not be complete without addressing these cases, which we discuss in the remarks of lines 235-244. We expand our treatment of these points, including the discussion below, in the revised version of the manuscript. If any point remains unclear, we are happy to provide further details.
To the reviewer's points:
1) **Bounds on $(\lambda^\star,\nu^\star)$**: Under Assumption 2.1, we provide such a bound on $\\|\lambda^\star\\|\_1$ in Lemma D.3. Explicitly, suppose there exists a $\mu^\dagger$ such that $\mathrm{KL}(\mu^\dagger \\| \pi) \leq C$, $\mathbb{E}\_{\mu^\dagger} [g\_i] \leq -\delta < 0$, and $\mathbb{E}\_{\mu^\dagger} [h\_j] = 0$ (Assumption 2.1). Then, by definition
$$
D^\star = d(\lambda^\star,\nu^\star) \leq \mathrm{KL}(\mu^\dagger \\| \pi) + \sum\_i \lambda^\star\_i \mathbb{E}\_{\mu^\dagger} [g\_i] + \sum\_j \nu^\star\_j \mathbb{E}\_{\mu^\dagger} [h\_j] \leq C - \\|\lambda^\star\\|\_1 \delta
$$
for the dual function in (7). By strong duality, $D^\star = P^\star \geq 0$, which yields
$$
\\|\lambda^\star\\|\_1 \leq \frac{C}{\delta}.
$$
A bound on $\nu^\star$, though more complicated, can be obtained under a similar constraint qualification (see, e.g., the general [Gauvin, "A necessary and sufficient regularity condition to have bounded multipliers in nonconvex programming," 1977]).
2) **Enforcing the bound**: Suppose $\\|\lambda^\star\\|\_1,\\|\nu^\star\\|\_1 \leq B$. Suffices it to replace steps 4 and 5 in Algorithm 1 by
$$
\lambda\_{k+1} = \Big[ \lambda\_k + \eta\_k g(x\_k) \Big]\_0^{B}
\quad \text{and} \quad
\nu\_{k+1} = \Big[ \nu\_k + \eta\_k h(x\_k) \Big]\_{-B}^{B},
$$
where $[z]\_{L}^U = \min(\max(z,L),U)$. This is a common solution used in the Euclidean case (see, e.g., [35]).
3) **How does this affect the guarantees?** It does not. Indeed, note that proof of Theorem 3.3 is based on the fact that (17) is a Lyapunov function. This is proved in Lemma C.1. When bounding the Euclidean norm of the dual variables (line 642), the first step is realizing that the projection $\lambda \mapsto [\lambda]\_0^B$ and $\nu \mapsto [\nu]\_{-B}^B$ are contractions, since $\lambda^\star\_i \in [0,B]$ and $\nu^\star\_j \in [-B,B]$. Hence, we can use
$$
\\| [\lambda]\_0^B - \lambda^\star \\|^2 \leq \\| \lambda - \lambda^\star \\|^2
\quad \text{and} \quad
\\| [\nu]\_{-B}^B - \nu^\star \\|^2 \leq \\| \nu - \nu^\star \\|^2
$$
to obtain (26). The proof then proceed as is and yields the exact same results as in Theorem 3.3, except that (13) [and similarly (14)] can now be rewritten as
$$
\frac{1}{K} \sum\_{k = 1}^K
\mathrm{KL}(\mu\_{k} \\| \mu^\star)
+ \frac{m}{2} W\_2^2(\mu\_k,\mu^\star)
\leq
\frac{R\_0^2}{\eta K}
+ \eta G^2 ( 3 + B^2 ).
$$