Latent Guided Sampling for Combinatorial Optimization

## Abstract

Abstract Combinatorial Optimization problems are widespread in domains such as logistics, manufacturing, and drug discovery, yet their NP-hard nature makes them computationally challenging. Recent Neural Combinatorial Optimization methods leverage deep learning to learn solution strategies, trained via Supervised or Reinforcement Learning (RL). While promising, these approaches often rely on task-specific augmentations, perform poorly on out-of-distribution instances, and lack robust inference mechanisms. Moreover, existing latent space models either require labeled data or rely on pre-trained policies. In this work, we propose LGS-Net, a novel latent space model that conditions on problem instances, and introduce an efficient inference method, Latent Guided Sampling (LGS), based on Markov Chain Monte Carlo and Stochastic Approximation. We show that the iterations of our method form a time-inhomogeneous Markov Chain and provide rigorous theoretical convergence guarantees. Empirical results on benchmark routing tasks show that our method achieves state-of-the-art performance among RL-based approaches.

### 1 Introduction

Combinatorial Optimization (CO) consists of finding the best solution from a discrete set of possibilities by optimizing a given objective function subject to constraints. It has widespread applications across various domains, including vehicle routing (Veres and Moussa, 2019), production planning (Dolgui et al., 2019), and drug discovery (Liu et al., 2017). However, its NP-hard nature and the complexity of many problem variants make solving CO problems highly challenging.
Traditional heuristic methods (e.g., (Kirkpatrick et al., 1983; Glover, 1989; Mladenović and Hansen, 1997)) rely on hand-crafted rules to guide the search, providing near-optimal solutions with significantly lower computational costs.
Inspired by the success of deep learning in computer vision (Krizhevsky et al., 2012; He et al., 2016) and natural language processing (Vaswani et al., 2017; Devlin, 2018), recent years have seen a surge in learning-based Neural Combinatorial Optimization (NCO) approaches for solving CO problems, including the Travelling Salesman Problem (TSP) and the Capacitated Vehicle Routing Problem (CVRP).

These methods leverage neural networks to learn a policy that generates solutions, trained via either Supervised Learning (SL) (Vinyals et al., 2015; Joshi et al., 2019; Hottung et al., 2021; Fu et al., 2021; Joshi et al., 2022; Kool et al., 2022) or Reinforcement Learning (RL) (Bello et al., 2016; Nazari et al., 2018; Kool et al., 2019; Chen and Tian, 2019; Kwon et al., 2020; Hottung and Tierney, 2020; Grinsztajn et al., 2023; Chalumeau et al., 2023). SL-based methods often struggle to obtain sufficient high-quality labeled data, whereas RL-based approaches can surpass them by exploring solutions autonomously.
Despite their success on in-distribution problem instances, these methods often generalize poorly to out-of-distribution cases, limiting their applicability in real-world scenarios.
Moreover, once a policy is trained, inference typically relies on relatively simple strategies such as stochastic sampling (Kool et al., 2019; Kwon et al., 2020), beam search (Steinbiss et al., 1994), or Monte Carlo Tree Search (Browne et al., 2012). A more powerful alternative, representing the current state-of-the-art in search-based RL (Bello et al., 2016; Hottung et al., 2022), is to actively fine-tune the policy for each new problem instance. However, this approach introduces significant computational and practical challenges.

Recently, a few methods (Hottung et al., 2021; Chalumeau et al., 2023) have explored learning a continuous latent space for discrete routing problems, enabling latent space optimization during inference with any continuous optimization technique.
However, the method proposed by Hottung et al. (2021) is limited by its reliance on a labeled training set for SL, while the one of Chalumeau et al. (2023) enforces independence between the problem instance and the latent space structure.
In contrast, we introduce a latent space model that conditions the latent representation on problem instances, addressing the limitations of previous approaches and enabling more effective latent space optimization.
In addition, most prior inference methods, including those not based on latent spaces, rely on the augmentation trick (Kwon et al., 2020) which enhances performance by generating variations of the same problem, such as rotating the coordinates, a task-specific technique that is only applicable to certain routing problems in an Euclidean space. While augmentation can improve performance, achieving competitive results without it is equally crucial for certain problems.
To address this, we propose a novel guided inference method designed for latent-based models, based on Markov Chain Monte Carlo (MCMC) and Stochastic Approximation (SA). Our method provides theoretical convergence guarantees and outperforms most state-of-the-art NCO methods.

More precisely, our contributions are summarized as follows.

- •
We introduce LGS-Net, a novel latent space model for Neural Combinatorial Optimization that requires neither labeled data nor a pretrained policy, learns a structured, instance-conditioned latent representation, and is supported by a rigorous mathematical framework.
- •
We propose LGS, a new inference scheme based on interacting MCMC, which jointly samples from the learned distribution while updating parameters via Stochastic Approximation. Moreover, we establish that its iterates form a joint time-inhomogeneous Markov Chain over the latent and solution spaces, converging to the desired target distribution.
- •
We empirically demonstrate that our approach sets a new state-of-the-art among RL-based CO methods, consistently outperforming existing techniques across diverse problem types, both with and without the augmentation trick.

### 2 Related Work

Neural Combinatorial Optimization.
The application of neural networks to Combinatorial Optimization problems dates back to Hopfield and Tank (1985) who employed a Hopfield network to solve small instances of the TSP. In recent years, numerous approaches leveraging advancements in deep learning have been developed to address Combinatorial Optimization challenges. A notable early contribution was made by Vinyals et al. (2015), who introduced Pointer Networks inspired by sequence-to-sequence architectures (Sutskever et al., 2014). Pointer Networks are particularly well-suited for Combinatorial Optimization tasks where the output sequence length depends on the input. These models were trained in a supervised manner to solve TSP instances using ground truth optimal or heuristic solutions. However, the reliance on supervised training restricted the model’s ability to surpass the quality of the provided training data. To address this limitation, Bello et al. (2016) extended this approach by employing RL techniques, such as Reinforce and actor-critic algorithms, to train the model without requiring ground truth solutions.

The CVRP was first addressed by Nazari et al. (2018), where the recurrent neural network (RNN) encoder of the Pointer Network was replaced with element-wise projections. Subsequently, Kool et al. (2019) proposed a modified transformer architecture that employs self-attention mechanisms (Vaswani et al., 2017).
Building on the Attention Model (AM) proposed by Kool et al. (2019), various AM variants have since been explored (Kwon et al., 2020; Xin et al., 2020; Kim et al., 2021; Xin et al., 2021; Kwon et al., 2021). For instance, Kwon et al. (2020) introduced POMO, an attention-based model that incorporates a more robust learning and inference strategy grounded in multiple optimal policies. Similarly, Grinsztajn et al. (2023) proposed a framework using a population of agents with multiple decoders, focusing exclusively on learning from the best-performing agent at each iteration.
In addition to these methods, various architectures based on graph neural networks have been proposed, particularly for solving the TSP. Notable examples include a graph embedding network (Khalil et al., 2017), a graph attention network (Deudon et al., 2018), and a graph convolutional network (Joshi et al., 2019).

Several approaches combine heuristic algorithms, such as local search (Papadimitriou and Steiglitz, 1998; De Moura and Bjørner, 2008), with machine learning techniques to tackle routing problems. For example, Chen and Tian (2019); Lu et al. (2019) propose an RL-based improvement method that iteratively selects a region of the solution and applies a local heuristic determined by a trainable policy.
To improve the generalization ability of constructive methods during inference, various strategies have been introduced, such as Efficient Active Search (EAS) (Hottung et al., 2022) and Simulation-guided Beam Search (SGBS) (Choo et al., 2022). Notably, EAS builds on POMO by fine-tuning a subset of model parameters at inference time using Gradient Descent, in contrast to Active Search (Bello et al., 2016), which updates all model parameters.

Among the latent space models designed to map discrete routing problems to continuous spaces, Hottung et al. (2021) used a conditional variational autoencoder (CVAE) (Sohn et al., 2015) that maps solutions to a continuous latent space. However, their approach relies on supervised training, which is hindered by the significant cost of acquiring high-quality labeled data.
To overcome this limitation, Chalumeau et al. (2023) proposed COMPASS, an RL-based approach that adapts a pre-trained policy by learning the latent space. However, the latent space in COMPASS is independent of the problem instances, similar to a GAN (Goodfellow et al., 2014) without a discriminator.
In contrast, we introduce a latent space model that conditions the latent space on problem instances and removes the need for a pre-trained policy. This approach can be interpreted as a VAE with a modified encoder that conditions only on problem instances, rather than on both problem instances and solutions.
Furthermore, we propose an inference method based on MCMC and SA, rather than using Differential Evolution (DE) (Storn and Price, 1997) or Covariance Matrix Adaptation (CMA) (Hansen and Ostermeier, 2001), which were used in previous works on latent space models.

Monte Carlo Markov Chain.
The Metropolis-Hastings algorithm (Metropolis et al., 1953; Hastings, 1970) is one of the most widely used MCMC methods. Since its introduction, numerous variants have been developed, including the Gibbs sampler (Geman and Geman, 1984) and Hamiltonian Monte Carlo (Duane et al., 1987). These methods have been theoretically studied, particularly in terms of geometric ergodicity, under well-established drift and minorization conditions (Meyn and Tweedie, 1994; Baxendale, 2005).
A crucial factor influencing the performance of MCMC methods is the choice of the proposal distribution, which significantly affects the convergence rate (Gelman et al., 1997). Traditionally, tuning the proposal distribution relies on heuristics and manual adjustments. To address this limitation, adaptive MCMC methods have been developed, where the proposal distribution is adjusted dynamically based on previous samples (Haario et al., 2001). The ergodicity of adaptive MCMC methods incorporating Stochastic Approximation (Robbins and Monro, 1951) has been studied by Andrieu and Atchadé (2007); Andrieu and Moulines (2006).
Beyond adaptive MCMC, time-inhomogeneous MCMC methods, which extend adaptive MCMC, have also been explored (Andrieu et al., 2001; Douc et al., 2004).
In our setting, both the proposal and target distributions evolve dynamically, and existing results are therefore insufficient to establish convergence guarantees. Indeed, prior analyses (Andrieu et al., 2001; Douc et al., 2004) focus on continuous spaces and rely heavily on regularity assumptions, such as strict convexity and coercivity of the objective function, which do not hold in our case. To address this gap, we introduce new results for time-inhomogeneous MCMC methods, enabling us to derive convergence guarantees.

### 3 Notation and Background

#### 3.1 Notation

In the following, for all distribution $\mu$ (resp. probability density $p$) we write $\mathbb{E}_{\mu}$ (resp. $\mathbb{E}_{p}$) the expectation under $\mu$ (resp. under $p$). We may also write $\mathbb{E}_{x\sim\mu}$ for the expectation under $\mu$. Given a measurable space $(\mathsf{X},\mathcal{X})$, where $\mathcal{X}$ is a countably generated $\sigma$-algebra, let $\mathsf{F}(\mathsf{X})$ denote the set of all measurable functions defined on $(\mathsf{X},\mathcal{X})$.
Let $\mathsf{M}(\mathsf{X})$ be the set of $\sigma$-finite measures on $(\mathsf{X},\mathcal{X})$, and $\mathsf{M}_{1}(\mathsf{X})\subset\mathsf{M}(\mathsf{X})$ the probability measures. For all $f\in\mathsf{F}(\mathsf{X})$ and $\mu\in\mathsf{M}(\mathsf{X})$, we write $\mu(f)=\int f(x)\mu(dx)$.
For a Markov kernel $P$ on $(\mathsf{X},\mathcal{X})$ and $\mu\in\mathsf{M}_{1}(\mathsf{X})$, the composition $\mu P$ is defined as $\mu P:\mathcal{X}\ni A\mapsto\int\mu(dx)P(x,dy)\mathbf{1}_{A}(y)$.
For probability measures $\mu$ and $\nu$ defined on the same measurable space, the Total Variation (TV) is defined as
$\left\|\mu-\nu\right\|_{\mathrm{TV}}:=\sup_{A\in\mathcal{X}}|\mu(A)-\nu(A)|$. The $\mathrm{L}_{2}$-norm of a random variable $X$ is defined as $\|X\|_{\mathrm{L}_{2}}:=\left(\mathbb{E}[\|X\|^{2}]\right)^{1/2}$.
The Hadamard product of two vectors $u$ and $v$ is denoted by $u\odot v$.
For a sequence $\left(a_{m}\right)_{m\in\mathbb{N}}$, and all $u\leq v$, we write $a_{u:v}=\{a_{u},\dots,a_{v}\}$.
Table [4](#Ax1.T4) provides a summary of the notations used throughout the paper for ease of reference.

#### 3.2 Problem Setting

In a Combinatorial Optimization problem, the objective is to determine the best assignment of discrete variables that satisfies the constraints of the problem. Let $x$ represent a given problem instance and $y$ denote a solution. An instance $x=\{\mathsf{x}_{i}\}_{i=1}^{n}\in\mathsf{X}\subset\mathbb{R}^{n\times d_{x}}$ consists of a set of $n$ nodes, each represented by a feature vector $\mathsf{x}_{i}\in\mathbb{R}^{d_{x}}$, which encodes relevant information about the node. For a given instance $x$, we aim to find a solution $y^{*}$ that minimizes the associated cost function $C$:

$$ $y^{*}\in\underset{y\in\mathsf{Y}}{\mathrm{argmin}}\,C(y,x)\;,$ (1) $$

where $\mathsf{Y}$ denotes the discrete set of all feasible solutions for the given problem $x$. This setting covers problems such as the TSP, CVRP, Knapsack, and Job Scheduling. For instance, in TSP, $x$ represents the coordinates of all nodes and $\mathsf{Y}$ consists of all possible node permutations.
In CVRP, $x$ additionally includes demands, and $\mathsf{Y}$ comprises all feasible routes satisfying the capacity constraints. In both cases, the cost corresponds to the cumulative distance of the route. Further details on both problems are provided in Appendix [A.1](#A1.SS1).

#### 3.3 Constructive NCO Methods

Constructive NCO methods (Vinyals et al., 2015; Bello et al., 2016; Nazari et al., 2018; Kool et al., 2019; Kwon et al., 2020) generate solutions sequentially using a stochastic policy $p_{\theta}(y|x)$, which defines the probability of selecting a solution $y$ given a problem instance $x$. This policy is parameterized by $\theta\in\Theta$, where $\Theta$ is a parameter space, and factorized as:

$$ $p_{\theta}(y|x)=\prod_{t=1}^{T}p_{\theta}(\mathsf{y}_{t}|x,\mathsf{y}_{1:t-1})\;,$ $$

with the convention $p_{\theta}(\mathsf{y}_{1}|x,\mathsf{y}_{1:0})=p_{\theta}(\mathsf{y}_{1}|x)$, where $\mathsf{y}_{t}\in\{0,\ldots,n\}$ is the selected node at step $t$, $\mathsf{y}_{1:t-1}$ denotes the sequence of nodes selected up to step $t-1$, and $T$ is the total number of decoding steps.
Following Bello et al. (2016), writing $\mathbb{P}_{x}$ the distribution of the problem instances, the policy is trained via RL by minimizing an empirical estimate of the expected cost:

$$ $J(\theta)=\mathbb{E}_{x\sim\mathbb{P}_{x},y\sim p_{\theta}(.|x)}\left[C(y,x)\right]\;.$ $$

where $C(y,x)$ denotes the tour cost. This objective is optimized using RL techniques, such as REINFORCE (Williams, 1992) or Actor-Critic methods (Konda and Tsitsiklis, 1999).

### 4 Latent Guided Sampling

#### 4.1 Model

The proposed model introduces a continuous latent search space for routing problems similar to Hottung et al. (2021); Chalumeau et al. (2023), which can be efficiently explored by any continuous optimization method at inference time.
To achieve this, we model the target distribution that generates a solution $y$ given a problem instance $x$ as a latent-variable model:

$$ $p_{\theta,\phi}(y|x)=\int p_{\phi}(z|x)p_{\theta}(y|x,z)\,\mathrm{d}z\;.$ $$

The encoder $p_{\phi}(z|x)$ maps the problem instance $x$ to a continuous $d_{z}$-dimensional latent representation $z$. The decoder $p_{\theta}(y|x,z)$ then generates a solution $y$ conditioned on both $z$ and $x$. Both the encoder and decoder are parameterized by neural networks with learnable parameters $\phi\in\Phi$ and $\theta\in\Theta$, respectively.
The probability of the decoder generating a solution $y$ is then factorized as:

$$ $p_{\theta}(y|x,z)=\prod_{t=1}^{T}p_{\theta}(\mathsf{y}_{t}|\mathsf{y}_{1:t-1},z,x)\;,$ $$

where $\mathsf{y}_{t}\in\{0,\ldots,n\}$ is the selected node at step $t$, and $\mathsf{y}_{1:t-1}$ denotes the sequence of nodes selected up to step $t-1$. It is important to note that the encoder differs from the variational distribution $q_{\phi}(z|x,y)$ (Kingma and Welling, 2014); it corresponds to a CVAE-Opt (Hottung et al., 2021), which requires labeled data for training.

Figure: Figure 1: Our Latent Model Architecture with the Gaussian Reparameterization Trick
Refer to caption: /html/2506.03672/assets/x1.png

Our encoder architecture follows the general structure of Kool et al. (2019) but includes additional layers to compute the parameters of the encoder distribution.
To compute output probabilities, we use a single decoder layer with multi-head attention to enable efficient inference.
At step $0\leq t\leq T$, for all $i\in\{0,\ldots,n\}$, this layer computes the probability $p_{\theta}(\mathsf{y}_{t}=i|\mathsf{y}_{1:t-1},x,z)$ while masking nodes that lead to infeasible solutions. Details on the encoder and decoder are provided in Appendix [A.2](#A1.SS2).

#### 4.2 Training

During training, our objective is to minimize the cost $C$ while encouraging diversity in the generated solutions to improve inference efficiency. To achieve this, we introduce an entropic regularization term (Ziebart et al., 2008; Haarnoja et al., 2018) controlled by a parameter $\beta$. The training loss is given by:

$$ $\mathcal{L}(\theta,\phi;x)=\sum_{k=1}^{K}\mathbb{E}_{z^{k}\sim p_{\phi}(\cdot\mid x)}\left[\mathbb{E}_{y^{k}\sim p_{\theta}(\cdot\mid x,z^{k})}\left[w^{k}C(y^{k},x)\right]+\beta\mathcal{H}(p_{\theta}(\cdot\mid x,z^{k}))\right]\;.$ (2) $$

where $\mathcal{H}(p_{\theta}(\cdot|x,z))$ represents the entropy of the conditional decoder distribution $p_{\theta}(y|x,z)$. The loss is similar to that of Maximum Entropy RL (Ziebart et al., 2008), but with a weighted cost for each latent sample.
The weights $w^{k}$ are defined as $w^{k}=\mathrm{exp}(-C(y^{k},x)/\tau)$, where $\tau$ controls the importance assigned to each latent sample. When $\tau$ is large, the model treats all latent samples nearly equivalently, while smaller values of $\tau$ give greater priority to the best latent sample, as described in Grinsztajn et al. (2023). The parameter $\tau$ is gradually decreased, encouraging the model to initially explore a broader range of latent samples before focusing more on the most promising ones.
The training procedure is detailed in Appendix [A.3](#A1.SS3).

#### 4.3 Inference

Any search procedure strategy can be used at inference time to find the best solution while keeping the computational cost manageable.
Possible approaches include evolutionary algorithms such as DE and CMA-ES, as well as learnable methods like Active Search (Bello et al., 2016) and Efficient Active Search (Hottung et al., 2022).
We formulate inference as a sampling problem: given an instance $x$, and the learned encoder and decoder parameters $\theta$ and $\phi$, our goal is to sample from the distribution:

$$ $\pi_{\theta}(y|x)\propto\int\pi_{\theta}(z,y|x)\mathrm{d}z\;,\quad\text{where}\quad\pi_{\theta}(z,y|x)\propto p_{\phi}(z|x)p_{\theta}(y|z,x)e^{-\lambda C(y,x)}\;.$ (3) $$

We omit the explicit dependence on $\phi$ since $p_{\phi}$ serves as a prior over latent variables.
To favor lower-cost solutions, we introduce the reweighting factor $\mathrm{exp}(-\lambda C(y,x))$, where $\lambda$ controls the trade-off between likelihood and cost.
However, incorporating this reweighting renders the distribution in ([3](#S4.E3)) intractable to sample from directly. While methods such as MCMC (Hastings, 1970) can be used to approximate it, they are often inefficient in practice.
To address this challenge, we propose Latent Guided Sampling (LGS), a novel inference method designed for latent space models. LGS constructs sequences of latent samples and corresponding solutions by running multiple interacting Markov Chains to encourage better exploration, while simultaneously updating the model parameters $\theta$ via Stochastic Approximation to maximize the following test objective:

$$ $\mathcal{L}_{test}(\theta;x)=\mathbb{E}_{\pi_{\theta}(\cdot\mid x)}\left[C(y,x)\right]\;.$ $$

Since the solution quality depends on the trained parameters, it is natural to iteratively update $\theta$. In contrast, the encoder parameters $\phi$ are kept fixed to avoid the high computational cost associated with backpropagating through them.
The gradient is estimated using previously sampled latent variables:

$$ $H_{\theta}\left(x,\{\left(z^{k},y^{k}\right)\}_{k=1}^{K}\right)=\frac{1}{K}\sum_{k=1}^{K}\left(C(y^{k},x)-b(x)\right)\nabla_{\theta}\log p_{\theta}(y^{k}|x,z^{k})\;,$ (4) $$

where $b(x)$ is a baseline function that reduces variance. The proposed method is detailed in Algorithm [1](#alg1).
The value of $K$ should be selected to balance stability (reducing variance in gradient estimates), effective exploration of the latent space, and computational efficiency.
Although standard gradient updates are used in Algorithm [1](#alg1), any other optimizer such as Adam (Kingma and Ba, 2015) could be used.
In the next section, we establish that the sequence generated by our algorithm forms a Markov Chain and converges to the target distribution, depending on the optimality of $\theta$.

Figure: Algorithm 1 Latent Guided Sampling

### 5 Theoretical Results

Let $\mathsf{Z}\subset\mathbb{R}^{d_{z}}$ be the latent space and $\mathsf{Y}$ the solution space.
In this section, we present theoretical results on our inference method described in Algorithm [1](#alg1).

#### 5.1 Convergence Analysis for Fixed θ \theta

We first analyze convergence in the absence of the Stochastic Approximation step, that is, without updating the parameter $\theta$ (line 10 in Algorithm [1](#alg1)).
Specifically, we show that the sequences generated by our algorithm form a Markov Chain and exhibit geometric convergence to the joint target distribution defined in ([3](#S4.E3)).

###### Proposition 5.1

The sequence $\{(Z_{m},Y_{m}):m\in\mathbb{N}\}$ generated by Algorithm [1](#alg1) for a fixed parameter $\theta$ forms a Markov Chain with transition kernel $P_{\theta}$.

The explicit expression for $P_{\theta}$ is provided in Proposition [C.1](#A3.ThmProposition1).
Consider the following assumptions.

###### Assumption 1

The cost function $C$ is bounded. For all $\phi\in\Phi$, the encoder distribution $p_{\phi}$ is positive.
Furthermore, the decoder probability satisfies
for all $1\leq t\leq T$, and all $(\mathsf{y}_{1:t-1},x,z)$, $p_{\theta}(\mathsf{y}_{t}=i|\mathsf{y}_{1:t-1},x,z)>0$ for all feasible nodes $i\in\{0,\ldots,n\}$.

Since $\mathsf{Y}$ is discrete, the boundedness of $C$ is a natural assumption, analogous to bounded rewards in RL (Fallah et al., 2021). A Gaussian choice for $p_{\phi}$ is standard, enabling the reparameterization trick (Kingma and Welling, 2014) for efficient gradient backpropagation.
To ensure positivity of the decoder, a common approach is to use a softmax function in the final layer of the neural network, which is a standard practice in most architectures.

###### Assumption 2

The proposal density $q$ is positive and symmetric.

Assumption [2](#Thmassumption2) on the proposal density $q$ is commonly used in various sampling-based methods, such as Importance Sampling and MCMC (Douc et al., 2004). It holds for a wide range of distributions, including Gaussian, Laplace, and Uniform.

###### Theorem 5.1

Let Assumptions [1](#Thmassumption1) and [2](#Thmassumption2) hold. Then, the Markov Kernel $P_{\theta}$ admits a unique invariant probability measure $\pi_{\theta}$, defined in ([3](#S4.E3)). There exist constants $\rho_{1},\rho_{2}\in(0,1)$ and $\kappa_{1},\kappa_{2}\in\mathbb{R}_{+}$ such that for all $\mu\in\mathsf{M}_{1}(\mathsf{Z}\times\mathsf{Y})$, $\theta\in\Theta$, and $m\in\mathbb{N}$,

$$ $\|\mu P_{\theta}^{m}-\pi_{\theta}\|_{\text{TV}}\leq\kappa_{1}\rho_{1}^{m}\quad\text{and}\quad\|\mu P_{\theta}^{m}-\pi_{\theta}\|_{\mathrm{L}_{2}}\leq\kappa_{2}\rho_{2}^{m}\|\mu-\pi_{\theta}\|_{\mathrm{L}_{2}}\;.$ $$

Theorem [5.1](#S5.Thmtheorem1) shows that the Markov Chain generated by our algorithm with a fixed $\theta$ converges geometrically to the joint target distribution in both Total Variation and $\mathrm{L}_{2}$-distance. Specifically, when the initial distribution is $\mu=p_{\phi}p_{\theta}$, the theorem guarantees rapid mixing of the Markov Chain, provided that the model is well-trained.

###### Proposition 5.1

###### Assumption 1

###### Assumption 2

###### Theorem 5.1

#### 5.2 Convergence Analysis for Adaptive θ \theta

Incorporating SA steps introduce a time-inhomogeneous Markov Chain, where both the Markov kernel and the target distribution evolve dynamically. While only a few results are available in this setting (Douc et al., 2004), we establish new convergence results for time-inhomogeneous Markov Chains under assumptions adapted to our setting, without relying on strict convexity or coercivity of the cost function $C$.
For simplicity, we only present the results relating to our setting here; more general results are provided in Appendix [D.1](#A4.SS1).
To analyze the convergence, we introduce the following additional assumptions.

###### Assumption 3

There exists $L\in\mathsf{F}(\mathsf{X}\times\mathsf{Z}\times\mathsf{Y})$ such that for all $x\in\mathsf{X},y\in\mathsf{Y},z\in\mathsf{Z}$ and $\theta\in\Theta$,

$$ $\left\|\nabla_{\theta}\log p_{\theta}(y|z,x)\right\|\leq L(x,y,z)\;.$ $$

Assumption [3](#Thmassumption3) is commonly used in the analysis of convergence rates of policies (Papini et al., 2018; Surendran et al., 2025).
In Adaptive MCMC, the Lipschitz condition is often applied to the Markov kernel rather than the target distribution (Andrieu and Atchadé, 2007; Andrieu and Moulines, 2006). However, since here both the kernel and the target distribution evolve dynamically, this Lipschitz condition with respect to the target distribution is more appropriate.

###### Assumption 4

There exists $\theta_{\infty}\in\Theta$ and a positive sequence $\left(a_{m}\right)_{m\in\mathbb{N}}$, with $a_{m}\to 0$ as $m\to\infty$ such that

$$ $\left\|\theta_{m}-\theta_{\infty}\right\|_{\mathrm{L}_{2}}^{2}=O(a_{m})\;.$ $$

Notably, we do not require $\theta_{\infty}$ to be a unique minimizer; it can simply be a critical point, which is often the case when the objective function in inference is non-convex. With additional regularity assumptions on the objective, this condition can be verified (see Appendix [D.3](#A4.SS3)).

###### Theorem 5.2

Let Assumptions [1](#Thmassumption1) - [4](#Thmassumption4) hold. Then, there exist a constant $\rho\in(0,1)$ and a positive sequence $\left(b_{m}\right)_{m\in\mathbb{N}}$ such that for all $\mu\in\mathsf{M}_{1}(\mathsf{Z}\times\mathsf{Y})$, and $m\in\mathbb{N}$,

$$ $\displaystyle\mathbb{E}\big{[}\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{\infty}}\right\|_{\mathrm{TV}}\big{]}$ $\displaystyle=\mathcal{O}\bigg{(}\rho^{b_{m}}+\sum_{j=m-b_{m}}^{m-1}\gamma_{j+1}+a_{m}\bigg{)}\;.$ $$

Furthermore, if $\limsup_{m\to\infty}\left(b_{m}^{-1}+b_{m}/m+b_{m}\gamma_{m}\right)=0$, then:

$$ $\mathbb{E}\left[\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{\infty}}\right\|_{\mathrm{TV}}\right]\xrightarrow[m\to\infty]{}0\;.$ $$

Theorem [5.2](#S5.Thmtheorem2) establishes that the time-inhomogeneous Markov Chain generated by our algorithm converges to the joint target distribution $\pi_{\theta_{\infty}}$.
The bound in Theorem [5.2](#S5.Thmtheorem2) has three key components: $(i)$ the mixing error, which reflects how well the Markov Chain mixes from an arbitrary initial distribution; $(ii)$ the tracking error, which quantifies how much the stationary distribution shifts over time due to changes in the parameters; and $(iii)$ the optimization error, which measures the difference between the current parameters and their limiting value $\theta_{\infty}$.
These terms are interdependent: choosing a larger $b_{m}$ accelerates the convergence of the mixing error but may slow the convergence of the parameters, while the step size sequence $\gamma_{m}$ affects the convergence rate of the parameters $a_{m}$.
If $\limsup_{m\to\infty}\left(b_{m}^{-1}+b_{m}/m+b_{m}\gamma_{m}\right)=0$, the expected total variation distance between the Markov Chain and the target distribution tends to zero as $m\rightarrow\infty$, ensuring convergence.
If $\gamma_{m}=m^{-\gamma}$, then choosing

$$ $b_{m}=\left\lfloor\frac{-\gamma\log(m)}{\log(\rho)}\right\rfloor\;,$ $$

yields a convergence rate of $\mathcal{O}\left(m^{-\gamma}\log m+a_{m}\right)$.

###### Assumption 3

###### Assumption 4

###### Theorem 5.2

### 6 Experiments

In this section, we illustrate our method using two classic CO problems: the TSP and the CVRP. We evaluate performance using benchmark datasets from the literature Hottung et al. (2021), consisting of 1,000 instances drawn from a training distribution—100 nodes uniformly sampled within the unit square. To evaluate generalization, we also test on two out-of-distribution datasets with larger sizes of 125 and 150 nodes.
All experiments were conducted on a GPU cluster using a single NVIDIA RTX 6000 GPU. Our source code is publicly available at(^1^11https://github.com/SobihanSurendran/LGS), with training and inference details in Appendix [E](#A5). In our experiments, we perform SA steps at selected intervals, rather than at every iteration, to reduce the computational cost and enable more extensive exploration of the latent space (see Appendix [E.2](#A5.SS2)).

Baselines. We compare our model to a range of state-of-the-art Reinforcement Learning methods and industrial solvers. These include Concorde (Applegate et al., 2006), an exact solver specialized for the TSP, LKH3 (Helsgaun, 2017), a leading solver for CO problems, and Google OR-Tools (Perron and Furnon, 2019), a widely used suite of optimization tools. Among the machine learning-based methods, we evaluate our approach against POMO (Kwon et al., 2020), CVAE-Opt (Hottung et al., 2021), EAS (Hottung et al., 2022), and COMPASS (Chalumeau et al., 2023).

**Table 1: Experiment results on TSP without and with the augmentation trick.**
|  | Training distribution | Generalization |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | n = 100 | n = 125 | n = 150 |  |  |  |  |  |  |  |
|  | Method | Obj. | Gap | Time | Obj. | Gap | Time | Obj. | Gap | Time |
|  | Concorde | 7.752 | 0.00% | 8M | 8.583 | 0.00% | 12M | 9.346 | 0.00% | 17M |
|  | LKH3 | 7.752 | 0.00% | 47M | 8.583 | 0.00% | 73M | 9.346 | 0.00% | 99M |
| no aug. | POMO (greedy) | 7.785 | 0.429% | $<$1M | 8.640 | 0.664% | $<$1M | 9.442 | 1.022% | $<$1M |
| POMO (sampling) | 7.772 | 0.261% | 10M | 8.595 | 0.140% | 50M | 9.377 | 0.327% | 1H30 |  |
| CVAE-Opt | 7.779 | 0.348% | 15H | 8.646 | 0.736% | 21H | 9.482 | 1.454% | 30H |  |
| EAS | 7.767 | 0.197% | 20M | 8.607 | 0.280% | 30M | 9.387 | 0.434% | 40M |  |
| COMPASS | 7.753 | 0.014% | 20M | 8.586 | 0.035% | 30M | 9.358 | 0.128% | 40M |  |
| LGS-Net (ours) | 7.752 | 0.002% | 20M | 8.584 | 0.012% | 30M | 9.354 | 0.081% | 40M |  |
| aug. | POMO (greedy) | 7.762 | 0.132% | $<$1M | 8.607 | 0.280% | $<$1M | 9.397 | 0.541% | $<$1M |
| POMO (sampling) | 7.757 | 0.068% | 30M | 8.596 | 0.151% | 70M | 9.378 | 0.338% | 100M |  |
| EAS | 7.755 | 0.042% | 1H | 8.591 | 0.093% | 100M | 9.363 | 0.177% | 160M |  |
| COMPASS | 7.752 | 0.002% | 40M | 8.585 | 0.024% | 1H | 9.352 | 0.059% | 1H30 |  |
| LGS-Net (ours) | 7.752 | 0.000% | 40M | 8.583 | 0.001% | 1H | 9.349 | 0.027% | 1H30 |  |

**Table 2: Experiment results on CVRP without and with the augmentation trick.**
|  | Training distribution | Generalization |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | n = 100 | n = 125 | n = 150 |  |  |  |  |  |  |  |
| Method | Obj. | Gap | Time | Obj. | Gap | Time | Obj. | Gap | Time |  |
|  | LKH3 | 15.54 | 0.00% | 17H | 17.50 | 0.00% | 19H | 19.22 | 0.00% | 20H |
|  | OR Tools | 17.084 | 9.936% | 38M | 18.036 | 3.063% | 64M | 21.209 | 10.349% | 73M |
| no aug. | POMO (greedy) | 15.740 | 1.287% | $<$1M | 17.905 | 2.314% | $<$1M | 19.882 | 3.444% | $<$1M |
| POMO (sampling) | 15.633 | 0.598% | 10M | 17.687 | 1.069% | 12M | 19.597 | 1.961% | 17M |  |
| CVAE-Opt | 15.752 | 1.364% | 32H | 17.864 | 2.080% | 36H | 19.843 | 3.240% | 46H |  |
| EAS | 15.563 | 0.148% | 40M | 17.541 | 0.234% | 1H | 19.319 | 0.515% | 1H30 |  |
| COMPASS | 15.561 | 0.135% | 40M | 17.546 | 0.263% | 1H | 19.358 | 0.718% | 1H30 |  |
| LGS-Net (ours) | 15.524 | -0.102% | 40M | 17.496 | -0.022% | 1H | 19.286 | 0.343% | 1H30 |  |
| aug. | POMO (greedy) | 15.652 | 0.721% | 1M | 17.756 | 1.463% | 1M | 19.701 | 2.503% | 1M |
| POMO (sampling) | 15.567 | 0.174% | 40M | 17.595 | 0.543% | 1H15 | 19.476 | 1.332% | 2H |  |
| EAS | 15.508 | -0.205% | 80M | 17.466 | -0.194% | 2H10 | 19.212 | -0.041% | 3H20 |  |
| COMPASS | 15.531 | -0.057% | 80M | 17.512 | 0.068% | 2H10 | 19.318 | 0.509% | 3H20 |  |
| LGS-Net (ours) | 15.501 | -0.251% | 80M | 17.461 | -0.223% | 2H10 | 19.229 | 0.046% | 3H20 |  |

The average performance of each method is reported in Table [1](#S6.T1) (TSP) and Table [2](#S6.T2) (CVRP), both with and without the augmentation trick of Kwon et al. (2020). Overall, our approach achieves state-of-the-art performance across most settings.
For the TSP, our method produces near-optimal solutions even without augmentation, and consistently reaches optimality when the augmentation trick is applied. Moreover, it consistently outperforms other methods on out-of-distribution instances.
Latent-space-based models trained via RL (ours and COMPASS) outperform other baselines, underscoring the effectiveness of leveraging a learned latent representation to capture solution diversity without relying on problem-specific tricks. Importantly, our method surpasses COMPASS in all TSP settings.
Although EAS achieves reasonable performance, it is considerably more computationally expensive, as it requires gradient computation at every iteration.

For the CVRP without augmentation, our model again outperforms all baselines, including both COMPASS and EAS. It also surpasses the performance of LKH-3 on the instances with $n=100$ and $n=125$. When augmentation is applied, performance improves further. However, for $n=150$, EAS outperforms our method, likely due to the reduced number of latent samples imposed by the computational budget. In this case, exploring the latent space effectively may require a higher sample count. We note that, when solving one instance at a time (thus relaxing the budget constraint), our model can achieve even stronger performance under augmentation.
In summary, our method achieves the best results across all TSP settings and remains the top performer on CVRP in most cases. While COMPASS is the closest competitor for TSP, EAS performs more strongly than COMPASS on CVRP, but still falls short of our method in most settings.

Figure: Figure 2: Visualization of the 2-dimensional latent space ($d_{z}=2$) learned by our model on a problem instance. The plotted path illustrates the search trajectory leading to the best-found solution.
Refer to caption: /html/2506.03672/assets/x2.png

Table [3](#S6.T3) and Figure [3](#S6.F3) present the results of ablation studies, focusing on the comparison of different inference methods with our model on CVRP instances with $n=100$ for a fixed time. Among the methods evaluated, Parallel MCMC, Interacting MCMC, and LGS consistently outperform other techniques, particularly DE and CMA-ES, which are commonly used inference methods in continuous spaces. In contrast, Single MCMC struggles due to limited exploration, leading to poor performance. EAS also underperforms, as the initial particles are insufficiently effective, and adjusting the parameters does not significantly improve the solution. This highlights the critical importance of particle propagation and parameter learning in improving solution quality.

Figure: Table 3: Comparison of different inference methods with our model on CVRP with $n=100$
Refer to caption: /html/2506.03672/assets/x3.png

Additionally, the ”Sampling” method corresponds to direct sampling from the distribution defined in ([3](#S4.E3)), without the reweighting factor $\mathrm{exp}(-\lambda C(y,x))$. The results highlight the advantage of incorporating this reweighting factor, which improves the exploration of the latent space, as shown in Figure [2](#S6.F2).
While interacting MCMC benefits from faster mixing (Theorem [5.1](#S5.Thmtheorem1)), it is nevertheless outperformed by LGS, which achieves better cost convergence despite potentially slower mixing (Theorem [5.2](#S5.Thmtheorem2)). This contrast highlights the importance of updating $\theta$ during inference: without the SA step, interacting MCMC may converge to samples from a possibly inaccurate distribution, resulting in suboptimal solutions. Nonetheless, interacting MCMC remains competitive due to the mitigating effect of the reweighting factor.

### 7 Conclusion

This paper introduces LGS-Net, a novel latent space model for Neural Combinatorial Optimization that conditions directly on problem instances, thereby removing the need for labeled data and pretrained policies. We further propose a guided inference method that generates sequences of latent samples and corresponding solutions based on MCMC and SA. We establish that the iterates of our method form a time-inhomogeneous Markov Chain, with theoretical convergence guarantees.
We evaluate our approach on TSP and CVRP, setting a new benchmark for RL-based CO methods, both with and without domain-specific augmentations.
A promising direction for future work is to explore how frequently the parameters of the target distribution should be updated to balance between optimal convergence rate and computational efficiency.

### Acknowledgements

The PhD of Sobihan Surendran was funded by the Paris Region PhD Fellowship Program of Région Ile-de-France. We would like to thank SCAI (Sorbonne Center for Artificial Intelligence) for providing the computing clusters.

### References

- Andrieu and Atchadé (2007)
C. Andrieu and Y. F. Atchadé.
On the efficiency of adaptive mcmc algorithms.
*Electronic Communications in Probability*, 12:336–349, 2007.
- Andrieu and Moulines (2006)
C. Andrieu and É. Moulines.
On the ergodicity properties of some adaptive mcmc algorithms.
*Annals of Applied Probability*, 16(3):1462–1505, 2006.
- Andrieu et al. (2001)
C. Andrieu, L. A. Breyer, and A. Doucet.
Convergence of simulated annealing using foster-lyapunov criteria.
*Journal of Applied Probability*, 38(4):975–994, 2001.
- Applegate et al. (2006)
D. Applegate, R. Bixby, V. Chvatal, and W. Cook.
Concorde tsp solver, 2006.
- Baxendale (2005)
P. H. Baxendale.
Renewal theory and computable convergence rates for geometrically
ergodic markov chains.
*Annals of Applied Probability*, 15:700–738, 2005.
- Bello et al. (2016)
I. Bello, H. Pham, Q. V. Le, M. Norouzi, and S. Bengio.
Neural combinatorial optimization with reinforcement learning.
*arXiv preprint arXiv:1611.09940*, 2016.
- Browne et al. (2012)
C. B. Browne, E. Powley, D. Whitehouse, S. M. Lucas, P. I. Cowling,
P. Rohlfshagen, S. Tavener, D. Perez, S. Samothrakis, and S. Colton.
A survey of monte carlo tree search methods.
*IEEE Transactions on Computational Intelligence and AI in
games*, 4(1):1–43, 2012.
- Chalumeau et al. (2023)
F. Chalumeau, S. Surana, C. Bonnet, N. Grinsztajn, A. Pretorius, A. Laterre,
and T. Barrett.
Combinatorial optimization with policy adaptation using latent space
search.
In *Advances in Neural Information Processing Systems*,
volume 36, pages 7947–7959, 2023.
- Chen and Tian (2019)
X. Chen and Y. Tian.
Learning to perform local rewriting for combinatorial optimization.
In *Advances in Neural Information Processing Systems*,
volume 32, 2019.
- Choo et al. (2022)
J. Choo, Y.-D. Kwon, J. Kim, J. Jae, A. Hottung, K. Tierney, and Y. Gwon.
Simulation-guided beam search for neural combinatorial optimization.
In *Advances in Neural Information Processing Systems*,
volume 35, pages 8760–8772, 2022.
- De Moura and Bjørner (2008)
L. De Moura and N. Bjørner.
Z3: An efficient smt solver.
In *International conference on Tools and Algorithms for the
Construction and Analysis of Systems*, pages 337–340. Springer, 2008.
- Deudon et al. (2018)
M. Deudon, P. Cournut, A. Lacoste, Y. Adulyasak, and L.-M. Rousseau.
Learning heuristics for the tsp by policy gradient.
In *Integration of Constraint Programming, Artificial
Intelligence, and Operations Research: 15th International Conference, CPAIOR
2018, Delft, The Netherlands, June 26–29, 2018, Proceedings 15*, pages
170–181. Springer, 2018.
- Devlin (2018)
J. Devlin.
Bert: Pre-training of deep bidirectional transformers for language
understanding.
*arXiv preprint arXiv:1810.04805*, 2018.
- Dolgui et al. (2019)
A. Dolgui, D. Ivanov, S. P. Sethi, and B. Sokolov.
Scheduling in production, supply chain and industry 4.0 systems by
optimal control: fundamentals, state-of-the-art and applications.
*International journal of production research*, 57(2):411–432, 2019.
- Douc et al. (2004)
R. Douc, E. Moulines, and J. S. Rosenthal.
Quantitative bounds on convergence of time-inhomogeneous markov
chains.
*Annals of Applied Probability*, pages 1643–1665, 2004.
- Douc et al. (2018)
R. Douc, E. Moulines, P. Priouret, P. Soulier, R. Douc, E. Moulines,
P. Priouret, and P. Soulier.
*Markov chains: Basic definitions*.
Springer, 2018.
- Duane et al. (1987)
S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth.
Hybrid monte carlo.
*Physics letters B*, 195(2):216–222, 1987.
- Fallah et al. (2021)
A. Fallah, K. Georgiev, A. Mokhtari, and A. Ozdaglar.
On the convergence theory of debiased model-agnostic
meta-reinforcement learning.
In *Advances in Neural Information Processing Systems*,
volume 34, pages 3096–3107, 2021.
- Fu et al. (2021)
Z.-H. Fu, K.-B. Qiu, and H. Zha.
Generalize a small pre-trained model to arbitrarily large tsp
instances.
In *Proceedings of the AAAI conference on artificial
intelligence*, volume 35, pages 7474–7482, 2021.
- Gelman et al. (1997)
A. Gelman, W. R. Gilks, and G. O. Roberts.
Weak convergence and optimal scaling of random walk metropolis
algorithms.
*The annals of applied probability*, 7(1):110–120, 1997.
- Geman and Geman (1984)
S. Geman and D. Geman.
Stochastic relaxation, gibbs distributions, and the bayesian
restoration of images.
*IEEE Transactions on pattern analysis and machine
intelligence*, (6):721–741, 1984.
- Glover (1989)
F. Glover.
Tabu search—part i.
*ORSA Journal on computing*, 1(3):190–206,
1989.
- Goodfellow et al. (2014)
I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair,
A. Courville, and Y. Bengio.
Generative adversarial nets.
In *Advances in Neural Information Processing Systems*,
volume 27, 2014.
- Grinsztajn et al. (2023)
N. Grinsztajn, D. Furelos-Blanco, S. Surana, C. Bonnet, and T. Barrett.
Winner takes it all: Training performant rl populations for
combinatorial optimization.
In *Advances in Neural Information Processing Systems*,
volume 36, pages 48485–48509, 2023.
- Haario et al. (2001)
H. Haario, E. Saksman, and J. Tamminen.
An adaptive metropolis algorithm.
*Bernoulli*, 7(2):223–242, 2001.
- Haarnoja et al. (2018)
T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine.
Soft actor-critic: Off-policy maximum entropy deep reinforcement
learning with a stochastic actor.
In *International Conference on Machine Learning*, pages
1861–1870. Pmlr, 2018.
- Hansen and Ostermeier (2001)
N. Hansen and A. Ostermeier.
Completely derandomized self-adaptation in evolution strategies.
*Evolutionary computation*, 9(2):159–195,
2001.
- Hastings (1970)
W. K. Hastings.
Monte carlo sampling methods using markov chains and their
applications.
*Biometrika*, 1970.
- He et al. (2016)
K. He, X. Zhang, S. Ren, and J. Sun.
Deep residual learning for image recognition.
In *Proceedings of the IEEE conference on computer vision and
pattern recognition*, pages 770–778, 2016.
- Helsgaun (2017)
K. Helsgaun.
An extension of the lin-kernighan-helsgaun tsp solver for constrained
traveling salesman and vehicle routing problems.
*Roskilde: Roskilde University*, 12:966–980, 2017.
- Hopfield and Tank (1985)
J. J. Hopfield and D. W. Tank.
“neural” computation of decisions in optimization problems.
*Biological cybernetics*, 52(3):141–152,
1985.
- Hottung and Tierney (2020)
A. Hottung and K. Tierney.
Neural large neighborhood search for the capacitated vehicle routing
problem.
In *ECAI 2020*, pages 443–450. IOS Press, 2020.
- Hottung et al. (2021)
A. Hottung, B. Bhandari, and K. Tierney.
Learning a latent search space for routing problems using variational
autoencoders.
In *International Conference on Learning Representations*, 2021.
- Hottung et al. (2022)
A. Hottung, Y.-D. Kwon, and K. Tierney.
Efficient active search for combinatorial optimization problems.
In *International Conference on Learning Representations*, 2022.
- Huang and Belongie (2017)
X. Huang and S. Belongie.
Arbitrary style transfer in real-time with adaptive instance
normalization.
In *Proceedings of the IEEE international conference on computer
vision*, pages 1501–1510, 2017.
- Joshi et al. (2019)
C. K. Joshi, T. Laurent, and X. Bresson.
An efficient graph convolutional network technique for the travelling
salesman problem.
*arXiv preprint arXiv:1906.01227*, 2019.
- Joshi et al. (2022)
C. K. Joshi, Q. Cappart, L.-M. Rousseau, and T. Laurent.
Learning the travelling salesperson problem requires rethinking
generalization.
*Constraints*, 27(1):70–98, 2022.
- Karimi et al. (2019)
B. Karimi, B. Miasojedow, E. Moulines, and H.-T. Wai.
Non-asymptotic analysis of biased stochastic approximation scheme.
In *Conference on Learning Theory*, pages 1944–1974. PMLR,
2019.
- Khalil et al. (2017)
E. Khalil, H. Dai, Y. Zhang, B. Dilkina, and L. Song.
Learning combinatorial optimization algorithms over graphs.
In *Advances in Neural Information Processing Systems*,
volume 30, 2017.
- Kim et al. (2021)
M. Kim, J. Park, et al.
Learning collaborative policies to solve np-hard routing problems.
In *Advances in Neural Information Processing Systems*,
volume 34, pages 10418–10430, 2021.
- Kingma and Ba (2015)
D. P. Kingma and J. Ba.
Adam: A method for stochastic optimization.
In *International Conference on Learning Representations*, 2015.
- Kingma and Welling (2014)
D. P. Kingma and M. Welling.
Auto-encoding variational bayes.
In *International Conference on Learning Representations*, 2014.
- Kirkpatrick et al. (1983)
S. Kirkpatrick, C. D. Gelatt Jr, and M. P. Vecchi.
Optimization by simulated annealing.
*science*, 220(4598):671–680, 1983.
- Konda and Tsitsiklis (1999)
V. Konda and J. Tsitsiklis.
Actor-critic algorithms.
In *Advances in Neural Information Processing Systems*,
volume 12, 1999.
- Kool et al. (2019)
W. Kool, H. Van Hoof, and M. Welling.
Attention, learn to solve routing problems!
In *International Conference on Learning Representations*, 2019.
- Kool et al. (2022)
W. Kool, H. van Hoof, J. Gromicho, and M. Welling.
Deep policy dynamic programming for vehicle routing problems.
In *International conference on integration of constraint
programming, artificial intelligence, and operations research*, pages
190–213. Springer, 2022.
- Krizhevsky et al. (2012)
A. Krizhevsky, I. Sutskever, and G. E. Hinton.
Imagenet classification with deep convolutional neural networks.
In *Advances in Neural Information Processing Systems*,
volume 25, 2012.
- Kwon et al. (2020)
Y.-D. Kwon, J. Choo, B. Kim, I. Yoon, Y. Gwon, and S. Min.
Pomo: Policy optimization with multiple optima for reinforcement
learning.
In *Advances in Neural Information Processing Systems*,
volume 33, pages 21188–21198, 2020.
- Kwon et al. (2021)
Y.-D. Kwon, J. Choo, I. Yoon, M. Park, D. Park, and Y. Gwon.
Matrix encoding networks for neural combinatorial optimization.
In *Advances in Neural Information Processing Systems*,
volume 34, pages 5138–5149, 2021.
- Liu et al. (2017)
R. Liu, X. Li, and K. S. Lam.
Combinatorial chemistry in drug discovery.
*Current opinion in chemical biology*, 38:117–126,
2017.
- Lu et al. (2019)
H. Lu, X. Zhang, and S. Yang.
A learning-based iterative method for solving vehicle routing
problems.
In *International conference on learning representations*, 2019.
- Metropolis et al. (1953)
N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller.
Equation of state calculations by fast computing machines.
*The journal of chemical physics*, 21(6):1087–1092, 1953.
- Meyn and Tweedie (1994)
S. P. Meyn and R. L. Tweedie.
Computable bounds for geometric convergence rates of markov chains.
*The Annals of Applied Probability*, pages 981–1011, 1994.
- Meyn and Tweedie (2012)
S. P. Meyn and R. L. Tweedie.
*Markov chains and stochastic stability*.
Springer Science & Business Media, 2012.
- Mladenović and Hansen (1997)
N. Mladenović and P. Hansen.
Variable neighborhood search.
*Computers & operations research*, 24(11):1097–1100, 1997.
- Nazari et al. (2018)
M. Nazari, A. Oroojlooy, L. Snyder, and M. Takác.
Reinforcement learning for solving the vehicle routing problem.
In *Advances in Neural Information Processing Systems*,
volume 31, 2018.
- Papadimitriou and Steiglitz (1998)
C. H. Papadimitriou and K. Steiglitz.
*Combinatorial optimization: algorithms and complexity*.
Courier Corporation, 1998.
- Papini et al. (2018)
M. Papini, D. Binaghi, G. Canonaco, M. Pirotta, and M. Restelli.
Stochastic variance-reduced policy gradient.
In *International Conference on Machine Learning*, pages
4026–4035. PMLR, 2018.
- Perron and Furnon (2019)
L. Perron and V. Furnon.
OR-Tools.
[https://developers.google.com/optimization/](https://developers.google.com/optimization/), 2019.
- Robbins and Monro (1951)
H. Robbins and S. Monro.
A stochastic approximation method.
*The Annals of Mathematical Statistics*, pages 400–407,
1951.
- Roberts and Rosenthal (1997)
G. Roberts and J. Rosenthal.
Geometric ergodicity and hybrid markov chains.
*Electronic Communications in Probability*, 2(2):13–25, 1997.
- Sohn et al. (2015)
K. Sohn, H. Lee, and X. Yan.
Learning structured output representation using deep conditional
generative models.
In *Advances in Neural Information Processing Systems*,
volume 28, 2015.
- Steinbiss et al. (1994)
V. Steinbiss, B.-H. Tran, and H. Ney.
Improvements in beam search.
In *ICSLP*, volume 94, pages 2143–2146, 1994.
- Storn and Price (1997)
R. Storn and K. Price.
Differential evolution–a simple and efficient heuristic for global
optimization over continuous spaces.
*Journal of global optimization*, 11:341–359, 1997.
- Surendran et al. (2024)
S. Surendran, A. Fermanian, A. Godichon-Baggioni, and S. Le Corff.
Non-asymptotic analysis of biased adaptive stochastic approximation.
In *Advances in Neural Information Processing Systems*,
volume 37, pages 12897–12943, 2024.
- Surendran et al. (2025)
S. Surendran, A. Godichon-Baggioni, and S. Le Corff.
Theoretical convergence guarantees for variational autoencoders.
In *The 28th International Conference on Artificial Intelligence
and Statistics*, 2025.
- Sutskever et al. (2014)
I. Sutskever, O. Vinyals, and Q. V. Le.
Sequence to sequence learning with neural networks.
In *Advances in Neural Information Processing Systems*,
volume 27, 2014.
- Vaswani et al. (2017)
A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez,
Ł. Kaiser, and I. Polosukhin.
Attention is all you need.
In *Advances in Neural Information Processing Systems*,
volume 30, 2017.
- Veres and Moussa (2019)
M. Veres and M. Moussa.
Deep learning for intelligent transportation systems: A survey of
emerging trends.
*IEEE Transactions on Intelligent transportation systems*,
21(8):3152–3168, 2019.
- Vinyals et al. (2015)
O. Vinyals, M. Fortunato, and N. Jaitly.
Pointer networks.
In *Advances in Neural Information Processing Systems*,
volume 28, 2015.
- Williams (1992)
R. J. Williams.
Simple statistical gradient-following algorithms for connectionist
reinforcement learning.
*Machine learning*, 8:229–256, 1992.
- Xin et al. (2020)
L. Xin, W. Song, Z. Cao, and J. Zhang.
Step-wise deep learning models for solving routing problems.
*IEEE Transactions on Industrial Informatics*, 17(7):4861–4871, 2020.
- Xin et al. (2021)
L. Xin, W. Song, Z. Cao, and J. Zhang.
Multi-decoder attention model with embedding glimpse for solving
vehicle routing problems.
In *Proceedings of the AAAI Conference on Artificial
Intelligence*, volume 35, pages 12042–12049, 2021.
- Ziebart et al. (2008)
B. D. Ziebart, A. L. Maas, J. A. Bagnell, A. K. Dey, et al.
Maximum entropy inverse reinforcement learning.
In *Aaai*, volume 8, pages 1433–1438. Chicago, IL, USA, 2008.

## Part

## Part Supplementary Material

### Notation

**Table 4: Summary of notation used throughout the paper.**
| Object | Description |
| --- | --- |
| $x=\{\mathsf{x}_{i}\}_{i=1}^{n}\in\mathsf{X}\subset\mathbb{R}^{n\times d_{x}}$ | Problem instance |
| $y=\left(\mathsf{y}_{1},\cdots,\mathsf{y}_{T}\right)\in\mathsf{Y}\subset\left\{0,\cdots,n\right\}^{T}$ | Solution |
| $z\in\mathsf{Z}\subset\mathbb{R}^{d_{z}}$ | Latent variable |
| $(\mathsf{X},\mathcal{X})$ | Problem instance space $\mathsf{X}$ with Borel $\sigma$-algebra $\mathcal{X}=\mathcal{B}(\mathsf{X})$ |
| $(\mathsf{Y},\mathcal{Y})$ | Discrete solution space $\mathsf{Y}$ with the power set $\mathcal{Y}=\mathcal{P}(\mathsf{Y})$ |
| $(\mathsf{Z},\mathcal{Z})$ | Latent space $\mathsf{Z}$ with Borel $\sigma$-algebra $\mathcal{Z}=\mathcal{B}(\mathsf{Z})$ |
| $\mathbb{P}_{x}$ | Distribution over problem instances |
| $p_{\phi}(z|x)$ | Encoder distribution |
| $p_{\theta}(y|x,z)$ | Decoder distribution |
| $C$ | Cost function |
| $n$ | Number of nodes in the instance |
| $t$, $T$ | Decoding step index and horizon |
| $m$, $M$ | Inference step index and total iterations |
| $k$, $K$ | Particle index and total number of particles |
| $B$ | Batch size used during training |

For a given batch of problem instances, we denote by $x_{(i)}$ the $i$-th input in a training batch. The corresponding solution and latent variable samples are denoted by $y_{(i)}^{k}$ and $z_{(i)}^{k}$ respectively, where $k$ indexes multiple samples drawn for the same input $x_{(i)}$.
During inference, we denote by $y_{m}^{k}$ and $z_{m}^{k}$ the solution and latent variable of the $k$-th particle at the $m$-th inference iteration.

### Appendix A Problem and Model Description

#### A.1 Problem Setting

Traveling Salesman Problem (TSP).
A TSP instance $x=\{\mathsf{x}_{i}\}_{i=1}^{n}$ consists of a set of $n$ nodes, where the feature $\mathsf{x}_{i}$ corresponds to its coordinates $c_{i}\in\mathbb{R}^{2}$. The objective is to find a permutation $y=(\mathsf{y}_{1},\dots,\mathsf{y}_{n})$ of the nodes, where $\mathsf{y}_{t}\in\{1,\dots,n\}$ and $\mathsf{y}_{t}\neq\mathsf{y}_{t^{\prime}}$ for all $t\neq t^{\prime}$, that minimizes the total tour length:

$$ $C(y,x)=\sum_{i=1}^{n-1}\|x_{\mathsf{y}_{i+1}}-x_{\mathsf{y}_{i}}\|+\|x_{\mathsf{y}_{n}}-x_{\mathsf{y}_{1}}\|\;,$ (5) $$

where $\|\cdot\|$ denotes the Euclidean norm. Note that the number of decoder steps $T$ equals the number of nodes $n$, i.e., $T=n$.

Capacitated Vehicle Routing Problem (CVRP).
CVRP generalizes TSP by introducing a depot (indexed as 0) and multiple routes, each starting and ending at the depot. Each customer $i\in\{1,\dots,n\}$ has a demand $d_{i}>0$ and a location $c_{i}\in\mathbb{R}^{2}$, while the depot has $d_{0}=0$. A fleet of vehicles, each with a capacity $D>0$, serves the customers. The goal is to determine the minimum number of vehicles and the corresponding routes, ensuring that each customer is visited exactly once and that the total demand in each route does not exceed $D$: for any route $j$,

$$ $\sum_{i\in R_{j}}d_{i}\leq D\;,$ $$

where $R_{j}$ denotes the set of customers assigned to route $j$.

#### A.2 Model Architecture details

##### A.2.1 Encoder

Given $d_{x}$-dimensional input features $\mathsf{x}_{i}$, the encoder initially computes $d_{h}$-dimensional node embeddings $h^{(0)}_{i}$ through a learned linear projection using parameters $W_{0}$ and $b_{0}$:

$$ $h^{(0)}_{i}=W_{0}\mathsf{x}_{i}+b_{0}\;.$ $$

The embeddings are updated using $L$ attention layers, each consisting of two sublayers: a multi-head attention ($\mathsf{MHA}$) layer followed by a node-wise fully connected feed-forward ($\mathsf{FF}$) layer. Each sublayer adds a skip connection (He et al., 2016) and instance normalization ($\mathsf{InstanceNorm}$) (Huang and Belongie, 2017). Denoting $h^{(l)}_{i}$ as the node embeddings produced by layer $l\in\{1,\dots,L\}$, the updates are defined as follows:

$$ $\displaystyle\hat{h}_{i}^{(l+1)}$ $\displaystyle=\mathsf{InstanceNorm}\left(h^{(l)}_{i}+\mathsf{MHA}\left(h^{(l)}_{1},\ldots,h^{(l)}_{n}\right)\right),$ $\displaystyle h_{i}^{(l+1)}$ $\displaystyle=\mathsf{InstanceNorm}\left(\hat{h}_{i}^{(l+1)}+\mathsf{FF}(\hat{h}_{i}^{(l+1)})\right)\;.$ $$

Then, it computes an aggregated embedding $\bar{h}^{(L)}$ of the input graph as the mean of the final node embeddings $h^{(L)}_{i}$.
Finally, the encoder generates a latent space vector using a reparameterization trick:

$$ $z=\mu_{\phi}(x)+\sigma_{\phi}(x)\odot\varepsilon,\quad\varepsilon\sim\mathcal{N}(0,I_{d_{z}})\;,$ $$

where the mean $\mu_{\phi}(x)$ and log-variance $\log\sigma_{\phi}(x)^{2}$ of the conditional distribution $p_{\phi}(z|x)$ are given by:

$$ $\mu_{\phi}(x)=\mathsf{FF}\left(\bar{h}^{(L)}\right)\quad\quad\text{and}\quad\quad\log\sigma_{\phi}(x)^{2}=\mathsf{FF}\left(\bar{h}^{(L)}\right)\;.$ $$

##### A.2.2 TSP Decoder

The context $c_{t}$ for the TSP decoder at time $t$ is derived by combining the latent vector $z$ and the output up to time $t-1$. Specifically, the context is defined as:

$$ $c_{t}=\left[z,h^{(L)}_{\mathsf{y}_{t-1}},h^{(L)}_{\mathsf{y}_{0}}\right]\;,$ $$

where $h^{(L)}_{\mathsf{y}_{0}}$ and $h^{(L)}_{\mathsf{y}_{t-1}}$ represent the embeddings of the starting node and the previously selected node, respectively.

Computation of Log-Probabilities.
The output probabilities for the TSP decoder are computed using a single decoder layer with multi-head attention.
This layer computes probabilities while incorporating a masking mechanism:

$$ $p_{\theta}(\mathsf{y}_{t}=i|\mathsf{y}_{0:t-1},x,z)=\begin{cases}\text{softmax }\left(\omega\tanh\left(\frac{q^{T}_{(c_{t})}k_{i}}{\sqrt{d_{k}}}\right)\right)&\text{if }i\neq\mathsf{y}_{s}\quad\forall s<t\;,\\ 0&\text{otherwise}\;,\\ \end{cases}$ $$

where the query and key are given by $q_{(c_{t})}=\mathsf{MHA}\left(c_{t},\{h^{(L)}_{i}\}_{i=1}^{n},\{h^{(L)}_{i}\}_{i=1}^{n}\right)$ and $k_{i}=W^{K}h^{(L)}_{i}$ respectively.

##### A.2.3 CVRP Decoder

Similar to the TSP decoder, the context $c_{t}$ for the CVRP decoder at time $t$ is derived by combining the latent vector $z$ and the output up to time $t-1$. Specifically, the context is defined as:

$$ $c_{t}=\left[z,h^{(L)}_{\mathsf{y}_{t-1}},\widehat{D}_{t}\right]\;,$ $$

where we keep track of the remaining vehicle capacity $\widehat{D}_{t}$ at time $t$. At $t=1$, this is initialized as $\widehat{D}_{t}=D$, after which it is updated as follows:

$$ $\widehat{D}_{t+1}=\begin{cases}\max(\widehat{D}_{t}-d_{\mathsf{y}_{t},t},0)&\text{if }\mathsf{y}_{t}\neq 0\;,\\ D&\text{if }\mathsf{y}_{t}=0\;.\end{cases}$ $$

Computation of Log-Probabilities.
The output probabilities for the TSP decoder are computed using a single decoder layer with multi-head attention. This layer computes log-probabilities while incorporating a masking mechanism:

$$ $p_{\theta}(\mathsf{y}_{t}=i|\mathsf{y}_{0:t-1},x,z)=\begin{cases}\text{softmax }\left(\omega\tanh\left(\frac{q^{T}_{(c_{t})}k_{i}}{\sqrt{d_{k}}}\right)\right)&\text{if }i\neq\mathsf{y}_{s}\quad\forall s<t\quad\text{and}\quad d_{i,t}\leq\widehat{D}_{t}\;,\\ 0&\text{otherwise}\;,\\ \end{cases}$ $$

where the query and key are given by $q_{(c_{t})}=\mathsf{MHA}\left(c_{t},\{h^{(L)}_{i}\}_{i=1}^{n},\{h^{(L)}_{i}\}_{i=1}^{n}\right)$ and $k_{i}=W^{K}h^{(L)}_{i}$ respectively.

#### A.3 Training

The estimator of the gradient of the objective defined in ([2](#S4.E2)) is computed using the Monte Carlo method:

$$ $\displaystyle\widehat{\nabla}_{\theta}\mathcal{L}(\theta,\phi)$ $\displaystyle=\frac{1}{B}\sum_{i=1}^{B}\sum_{k=1}^{K}w_{(i)}^{k}\left(C(y_{(i)}^{k},x_{(i)})-b(x_{(i)})\right)\nabla_{\theta}\log p_{\theta}(y_{(i)}^{k}|x_{i},z_{(i)}^{k})$ (6) $\displaystyle\quad-\beta\frac{1}{B}\sum_{i=1}^{B}\sum_{k=1}^{K}\log p_{\theta}(y_{(i)}^{k}|x_{(i)},z_{(i)}^{k})\nabla_{\theta}\log p_{\theta}(y_{(i)}^{k}|x_{i},z_{(i)}^{k})\;,$ $$

$$ $\widehat{\nabla}_{\phi}\mathcal{L}(\theta,\phi)=\frac{1}{B}\sum_{i=1}^{B}\sum_{k=1}^{K}w_{(i)}^{k}\left(C(y_{(i)}^{k},x_{(i)})-b(x_{(i)})\right)\nabla_{\phi}\log p_{\phi}(z_{(i)}^{k}|x_{(i)})\;,$ (7) $$

where $b(x)$ denotes a baseline function that does not depend on $y$.
This update rule has an intuitive interpretation: it adjusts the parameters $\theta$ and $\phi$ in directions that favor solutions yielding the highest reward, while simultaneously constraining the latent space to remain bounded and encouraging diversity in the sampled trajectories. The training procedure is outlined in Algorithm [2](#alg2). The update step $\operatorname{ADAM}$ for the parameters $(\theta,\phi)$ corresponds to a single Adam update step (Kingma and Ba, 2015).

Figure: Algorithm 2 REINFORCE training

### Appendix B Preliminaries on Markov Chains

In this section, we use the following definitions.

- •
A sequence of random variables $\left\{X_{n},n\in\mathbb{N}\right\}$ is a Markov Chain with respect to the filtration $(\mathcal{F}_{n})_{n\geq 0}$ with Markov kernel $P:\mathsf{X}\times\mathcal{X}\rightarrow\mathbb{R}_{+}$ if for any bounded measurable function $f:\mathsf{X}\rightarrow\mathbb{R}$,
$\mathbb{E}\left[f\left(X_{n+1}\right)\mid\mathcal{F}_{n}\right]=Pf\left(X_{n}\right)=\int f(x)P\left(X_{n},\mathrm{d}x\right)\;.$
- •
Furthermore, the sequence $\left\{X_{n},n\in\mathbb{N}\right\}$ is a state-dependent Markov Chain if for any bounded measurable function $f:\mathsf{X}\rightarrow\mathbb{R}$,
$\mathbb{E}\left[f\left(X_{n+1}\right)\mid\mathcal{F}_{n}\right]=P_{\theta_{n}}f\left(X_{n}\right)=\int f(x)P_{\theta_{n}}\left(X_{n},\mathrm{d}x\right)\;,$
where $P_{\theta_{n}}:\mathsf{X}\times\mathcal{X}\rightarrow\mathbb{R}_{+}$ is a Markov kernel with controlled parameters $\theta_{n}\in\mathbb{R}^{d}$.

###### Definition B.1 (Invariant Probability Measure)

A probability measure $\pi$ on $(\mathsf{X},\mathcal{X})$ is called invariant for the Markov kernel $P$ if it satisfies $\pi P=\pi$.

If $\left\{X_{n}:n\in\mathbb{N}\right\}$ is a Markov Chain with Markov kernel $P$ and $X_{0}$ is distributed according to an invariant probability measure $\pi$, then for all $n\geq 1$, we have $X_{n}\sim\pi$.

###### Definition B.2 (Reversibility)

A Markov kernel $P$ on $(\mathsf{X},\mathcal{X})$ is said to be reversible with respect to a probability measure $\pi$ if and only if

$$ $\pi(\mathrm{d}x)P(x,\mathrm{d}y)=\pi(\mathrm{d}y)P(y,\mathrm{d}x).$ $$

###### Definition B.3 (Coupling of Probability Measures)

Let $(\mathsf{X},\mathcal{X})$ be a measurable space and let $\mu,\nu$ be two probability measures, i.e., $\mu,\nu\in\mathsf{M}_{1}(\mathsf{X})$. We define $\mathcal{C}(\mu,\nu)$, the coupling set associated with $(\mu,\nu)$, as follows:

$$ $\mathcal{C}(\mu,\nu)=\left\{\zeta\in\mathsf{M}_{1}(\mathsf{X}^{2}):\forall A\in\mathcal{X},\ \zeta(A\times\mathsf{X})=\mu(A),\ \zeta(\mathsf{X}\times A)=\nu(A)\right\}\;.$ $$

###### Definition B.4 (Total Variation Distance)

Let $(\mathsf{X},\mathcal{X})$ be a measurable space and let $\mu,\nu$ be two probability measures in $\mathsf{M}_{1}(\mathsf{X})$. The total variation norm between $\mu$ and $\nu$, denoted by $\|\mu-\nu\|_{\mathrm{TV}}$, is defined by

$$ $\displaystyle\|\mu-\nu\|_{\mathrm{TV}}$ $\displaystyle=2\sup\left\{|\mu(f)-\nu(f)|:f\in\mathsf{F}(\mathsf{X}),0\leq f\leq 1\right\}$ $\displaystyle=2\inf\left\{\zeta(\Delta):\zeta\in\mathcal{C}(\mu,\nu)\right\}\;,$ $$

where $\Delta(x,x^{\prime})=\mathbf{1}_{x\neq x^{\prime}}$ for all $(x,x^{\prime})\in\mathsf{X}^{2}$.

###### Assumption 5

Let $P$ be a Markov transition kernel on $(\mathsf{X},\mathcal{X})$. Suppose there exists a function $V:\mathsf{X}\to[0,\infty)$ satisfying $\sup_{x\in\mathsf{X}}V(x)<\infty$, and the following conditions hold.

- 1.
Minorization Condition. There exist $\mathcal{K}\in\mathcal{X}$, $\varepsilon>0$ and a probability measure $\nu$ such that $\nu(\mathcal{K})>0$ and, for all
$A\in\mathcal{X}$ and $x\in\mathcal{K}$,
$P(x,A)\geq\varepsilon\nu(A)\;.$
- 2.
Drift Condition (Foster-Lyapunov Condition). There exist constants $\lambda\in[0,1)$, $b\in(0,\infty)$ satisfying
$PV(x)\leq\begin{cases}\lambda V(x),&x\notin\mathcal{K}\;,\\
b,&x\in\mathcal{K}\;.\end{cases}$

###### Theorem B.1

(Meyn and Tweedie, 1994; Baxendale, 2005)
Let Assumption [5](#Thmassumption5) hold for a function $V:\mathsf{X}\to[0,\infty)$, where $\sup_{x\in\mathsf{X}}V(x)<\infty$. Then, the Markov kernel $P$ admits a unique invariant probability measure $\pi$. Moreover, $\pi(V)<\infty$ and there exist constants $(\rho,\kappa)\in(0,1)\times\mathbb{R}_{+}$ such that for all $\mu\in\mathsf{M}_{1}(\mathsf{X})$, and $m\in\mathbb{N}$,

$$ $\|\mu P^{m}-\pi\|_{\text{TV}}\leq\kappa\rho^{m}\mu(V)\;.$ $$

This result was originally stated and proven in (Meyn and Tweedie, 1994, Theorem 2.3) with explicit formulas for $\rho$ and $\kappa$, and was later improved in (Baxendale, 2005, Section 2.1).
Further details are available in (Meyn and Tweedie, 2012, Chapter 15) and (Douc et al., 2018, Chapter 15).

###### Definition B.1 (Invariant Probability Measure)

###### Definition B.2 (Reversibility)

###### Definition B.3 (Coupling of Probability Measures)

###### Definition B.4 (Total Variation Distance)

###### Assumption 5

###### Theorem B.1

### Appendix C Convergence Analysis for Fixed θ \theta

In this section, we prove that the iterates of Algorithm [1](#alg1) with $K=1$ and fixed $\theta$ (the exact algorithm is given in Algorithm [3](#alg3)) form a reversible Markov Chain (Proposition [5.1](#S5.ThmProposition1)) and exhibit geometric ergodicity toward the joint target distribution (Theorem [5.1](#S5.Thmtheorem1)). We then extend these results to the general case of arbitrary $K$ in Section [C.3](#A3.SS3).

Figure: Algorithm 3 Latent Guided Sampling ($K=1$)

#### C.1 Proof of Proposition 5.1

###### Proposition C.1

The sequence $\{(Z_{m},Y_{m}):m\in\mathbb{N}\}$ generated by Algorithm [3](#alg3) with $K=1$ and fixed $\theta$ forms a Markov Chain with transition kernel $P_{\theta}$. Moreover, $P_{\theta}$ is $\pi_{\theta}$-reversible and for all $z\in\mathsf{Z}$, $y\in\mathsf{Y}$, and $A\in\mathcal{Z}\times\mathcal{Y}$,

$$ $P_{\theta}\big{(}(z,y),A\big{)}=\int_{A}q(dz^{\prime}|z)p_{\theta}(dy^{\prime}|z^{\prime},x)\alpha(z,y,z^{\prime},y^{\prime})+\bar{\alpha}_{\theta}(z,y)\delta_{(z,y)}(A)\;,$ (8) $$

where

$$ $\displaystyle\alpha(z,y,z^{\prime},y^{\prime})$ $\displaystyle=\min\left(1,e^{-\lambda(C(y^{\prime},x)-C(y,x))}\frac{p_{\phi}(z^{\prime}|x)}{p_{\phi}(z|x)}\right)\;,$ $\displaystyle\bar{\alpha}_{\theta}(z,y)$ $\displaystyle=1-\int_{\mathsf{Z}\times\mathsf{Y}}q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)\alpha(z,y,z^{\prime},y^{\prime})\;.$ $$

###### Proof.

To compute the Markov kernel for the joint chain $\{(Z_{m},Y_{m}):m\in\mathbb{N}\}$, we introduce the filtration:

$$ $\mathcal{F}_{m}=\sigma(Z_{0},Y_{0},U_{1:m})\;,$ $$

where $U_{1:m}=\left(U_{1},\cdots,U_{m}\right)$ denotes the sequence of uniform random variables.
For all bounded or non-negative measurable function $h$ on $\mathsf{Z}\times\mathsf{Y}$ and all $m\in\mathbb{N}$,

$$ $\displaystyle\mathbb{E}\big{[}h(Z_{m+1},Y_{m+1})\mid\mathcal{F}_{m}\big{]}$ $\displaystyle\quad=\mathbb{E}\left[1_{\left\{U_{m+1}<\alpha(Z_{m},Y_{m},Z_{m+1}^{\prime},Y_{m+1}^{\prime})\right\}}h(Z_{m+1}^{\prime},Y_{m+1}^{\prime})\mid\mathcal{F}_{m}\right]$ $\displaystyle\qquad+\mathbb{E}\left[1_{\left\{U_{m+1}\geq\alpha(Z_{m},Y_{m},Z_{m+1}^{\prime},Y_{m+1}^{\prime})\right\}}h(Z_{m},Y_{m})\mid\mathcal{F}_{m}\right]$ $\displaystyle\quad=\int_{\mathsf{Z}\times\mathsf{Y}}q(\mathrm{d}z^{\prime}|Z_{m})p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)\alpha(Z_{m},Y_{m},z^{\prime},y^{\prime})h(z,y)+\bar{\alpha}_{\theta}(Z_{m},Y_{m})h(Z_{m},Y_{m})\;,$ $$

where

$$ $\bar{\alpha}_{\theta}(Z_{m},Y_{m})=1-\int_{\mathsf{Z}\times\mathsf{Y}}q(\mathrm{d}z^{\prime}|Z_{m})p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)\alpha(Z_{m},Y_{m},z^{\prime},y^{\prime})\;.$ $$

Thus, $\{(Z_{m},Y_{m}):m\in\mathbb{N}\}$ forms a Markov Chain with the transition kernel given, for all $z\in\mathsf{Z}$, $y\in\mathsf{Y}$, and $A\in\mathcal{Z}\times\mathcal{Y}$, by

$$ $P_{\theta}\big{(}(z,y),A\big{)}=\int_{A}q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)\alpha(z,y,z^{\prime},y^{\prime})+\bar{\alpha}_{\theta}(z,y)\delta_{(z,y)}(A)\;.$ $$

The reversibility of the Markov transition kernel $P_{\theta}$ follows directly from its construction as a Metropolis–Hastings algorithm (Douc et al., 2018). Nonetheless, we include a proof specific to our setting by verifying the detailed balance condition:

$$ $\pi_{\theta}(z,y)P_{\theta}\big{(}(z,y),(z^{\prime},y^{\prime})\big{)}=\pi_{\theta}(z^{\prime},y^{\prime})P_{\theta}\big{(}(z^{\prime},y^{\prime}),(z,y)\big{)}\;.$ (9) $$

Define the ratio in the acceptance probability $\alpha$ as $r$:

$$ $r(z,y,z^{\prime},y^{\prime})=e^{-\lambda\big{(}C(y^{\prime},x)-C(y,x)\big{)}}\frac{p_{\phi}(z^{\prime}|x)}{p_{\phi}(z|x)}\;.$ $$

We separate the analysis in two cases depending on the value of $r(z,y,z^{\prime},y^{\prime})$.

Case 1.
If $r(z,y,z^{\prime},y^{\prime})\leq 1$, then $\alpha(z,y,z^{\prime},y^{\prime})=r(z,y,z^{\prime},y^{\prime})$ and $\alpha(z^{\prime},y^{\prime},z,y)=1$. Thus,

$$ $\displaystyle\pi_{\theta}(z,y)P_{\theta}\big{(}(z,y),(z^{\prime},y^{\prime})\big{)}$ $\displaystyle=p_{\phi}(z|x)p_{\theta}(y|z,x)e^{-\lambda C(y,x)}P_{\theta}\big{(}(z,y),(z^{\prime},y^{\prime})\big{)}$ $\displaystyle=p_{\phi}(z|x)p_{\theta}(y|z,x)e^{-\lambda C(y,x)}$ $\displaystyle\quad\times q(z^{\prime}|z)p_{\theta}(y^{\prime}|z^{\prime},x)e^{-\lambda(C(y^{\prime},x)-C(y,x))}\frac{p_{\phi}(z^{\prime}|x)}{p_{\phi}(z|x)}$ $\displaystyle=p_{\phi}(z^{\prime}|x)p_{\theta}(y^{\prime}|z^{\prime},x)e^{-\lambda C(y^{\prime},x)}q(z^{\prime}|z)p_{\theta}(y|z,x)$ $\displaystyle=\pi_{\theta}(z^{\prime},y^{\prime})P_{\theta}\big{(}(z^{\prime},y^{\prime}),(z,y)\big{)}\;.$ $$

Case 2.

If $r(z,y,z^{\prime},y^{\prime})>1$, then $\alpha(z,y,z^{\prime},y^{\prime})=1$ and $\alpha(z^{\prime},y^{\prime},z,y)=r(z^{\prime},y^{\prime},z,y)$. Similarly,

$$ $\displaystyle\pi_{\theta}(z^{\prime},y^{\prime})P_{\theta}\big{(}(z^{\prime},y^{\prime}),(z,y)\big{)}$ $\displaystyle=p_{\phi}(z^{\prime}|x)p_{\theta}(y^{\prime}|z^{\prime},x)e^{-\lambda C(y^{\prime},x)}P_{\theta}\big{(}(z^{\prime},y^{\prime}),(z,y)\big{)}$ $\displaystyle=p_{\phi}(z^{\prime}|x)p_{\theta}(y^{\prime}|z^{\prime},x)e^{-\lambda C(y^{\prime},x)}$ $\displaystyle\quad\times q(z|z^{\prime})p_{\theta}(y|z,x)e^{-\lambda(C(y,x)-C(y^{\prime},x))}\frac{p_{\phi}(z|x)}{p_{\phi}(z^{\prime}|x)}$ $\displaystyle=p_{\phi}(z|x)p_{\theta}(y|z,x)e^{-\lambda C(y,x)}q(z|z^{\prime})p_{\theta}(y^{\prime}|z^{\prime},x)$ $\displaystyle=\pi_{\theta}(z,y)P_{\theta}\big{(}(z,y),(z^{\prime},y^{\prime})\big{)}\;.$ $$

Since the detailed balance condition holds in both cases, we conclude that $P_{\theta}$ is $\pi_{\theta}$-reversible.
∎

###### Proposition C.1

###### Proof.

#### C.2 Proof of Theorem 5.1 for K = 1 K=1

###### Proof.

This follows from Theorem [B.1](#A2.Thmtheorem1), provided that we verify the minorization and drift conditions of Assumption [5](#Thmassumption5). We consider the set

$$ $\mathcal{K}=\left\{(z,y)\in\mathsf{Z}\times\mathsf{Y}\mid\|z\|^{2}\leq R^{2},C(y)\neq\max_{y^{\prime}\in\mathsf{Y}}C(y^{\prime})\right\}$ $$

which is compact since $\mathsf{Y}$ is finite. We denote by $B(c,R)=\{z\in\mathsf{Z}\mid\|z-c\|^{2}\leq R^{2}\}$ the ball centered at $c$ with radius $R$.

Minorization Condition.

The Markov transition kernel is given by: for all $A\in\mathcal{Z}\times\mathcal{Y}$,

$$ $\displaystyle P_{\theta}((z,y),A)$ $\displaystyle\quad=\int_{A}q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)\min\left(1,e^{-\lambda(C(y^{\prime},x)-C(y,x))}\frac{p_{\phi}(z^{\prime}|x)}{p_{\phi}(z|x)}\right)+\bar{\alpha}_{\theta}(z,y)\delta_{(z,y)}(A)$ $\displaystyle\quad\geq\int_{A\cap\mathcal{K}}q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)\min\left(1,e^{-\lambda(C(y^{\prime},x)-C(y,x))}\frac{p_{\phi}(z^{\prime}|x)}{p_{\phi}(z|x)}\right)+\bar{\alpha}_{\theta}(z,y)\delta_{(z,y)}(A)\;.$ $$

We now establish a minorization by separately analyzing each term.

- •
Proposal Component:
Since the proposal density $q$ is assumed to be positive on the compact set $B(0,R)$, for all $z,z^{\prime}\in B(0,R)$,
$q(z^{\prime}\mid z)\geq\varepsilon_{q}:=\inf_{z,z^{\prime}\in B(0,R)}q(z^{\prime}\mid z)>0.$
- •
Encoder Component:
By Assumption [1](#Thmassumption1), $p_{\phi}(\cdot|x)$ is positive, so that by applying a similar argument as above, we obtain the existence of $\varepsilon_{e}>0$ such that
$\frac{p_{\phi}(z^{\prime}|x)}{p_{\phi}(z|x)}\geq\varepsilon_{e}\;.$
- •
Decoder Component:
By Assumption [1](#Thmassumption1), the categorical transition probability satisfies:
$\inf_{z^{\prime}\in\mathsf{Z},y^{\prime}\in\mathsf{Y}}p_{\theta}(y^{\prime}|z^{\prime},x)\geq\varepsilon_{d}>0\;.$
- •
Exponential Weighting:
Since $e^{-\lambda(C(y^{\prime},x)-C(y,x))}$ is always positive and $C$ is bounded (Assumption [1](#Thmassumption1)), there exits $\varepsilon_{w}>0$ such that:
$e^{-\lambda(C(y^{\prime},x)-C(y,x))}\geq\varepsilon_{w}\;.$
(10)

Combining these bounds, for all $A\in\mathcal{Z}\times\mathcal{Y}$, we get:

$$ $\displaystyle P_{\theta}((z,y),A)$ $\displaystyle\quad=\int_{A}q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)\min\left(1,e^{-\lambda(C(y^{\prime},x)-C(y,x))}\frac{p_{\phi}(z^{\prime}|x)}{p_{\phi}(z|x)}\right)+\bar{\alpha}_{\theta}(z,y)\delta_{(z,y)}(A)$ $\displaystyle\quad\geq\mu^{\mathrm{Leb}}(B(0,R))|\mathsf{Y}|\varepsilon_{q}\varepsilon_{d}\min\left(1,\varepsilon_{w}\varepsilon_{e}\right)\int_{A}\nu_{q}(\mathrm{d}z^{\prime})\nu_{d}(\mathrm{d}y^{\prime})\;,$ $$

where $\nu_{C}$ is the uniform probability measure over $\mathsf{Y}$:

$$ $\nu_{d}(\mathrm{d}y^{\prime})=\frac{1}{|\mathsf{Y}|}\sum_{y\in\mathsf{Y}}\delta_{y}(\mathrm{d}y^{\prime})$ $$

and

$$ $\nu_{q}(\mathrm{d}z^{\prime}):=\frac{\mu^{\mathrm{Leb}}(\mathrm{d}z^{\prime})}{\mu^{\mathrm{Leb}}(B(0,R))}1_{B(0,R)}(z^{\prime}),$ $$

where $\mu^{\mathrm{Leb}}$ denotes the Lebesgue measure.
Thus, the transition kernel satisfies the uniform minorization condition:

$$ $P_{\theta}((z,y),A)\geq\varepsilon\nu(A)\;,$ (11) $$

where $\varepsilon=\mu^{\mathrm{Leb}}(B(0,R))|\mathsf{Y}|\varepsilon_{q}\varepsilon_{d}\min\left(1,\varepsilon_{w}\varepsilon_{e}\right)$ and $\nu(\mathrm{d}z^{\prime},\mathrm{d}y^{\prime})=\nu_{q}(\mathrm{d}z^{\prime})\nu_{d}(\mathrm{d}y^{\prime})$ is a probability measure.

Drift condition.

For a fixed $x\in\mathsf{X}$, we define the Lyapunov function for all $0<s\leq\lambda$ as

$$ $V_{x}(z,y)=e^{sC(y,x)}\;.$ $$

For all $\left(z,y\right)\in\mathcal{K}$, applying $P_{\theta}$ to $V_{x}$, we have:

$$ $\displaystyle P_{\theta}V_{x}(z,y)$ $\displaystyle=\int e^{sC(y^{\prime},x)}\left(q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)\alpha(z,y,z^{\prime},y^{\prime})+\bar{\alpha}_{\theta}(z,y)\delta_{(z,y)}(\mathrm{d}z^{\prime},\mathrm{d}y^{\prime})\right)$ $\displaystyle=\int e^{sC(y^{\prime},x)}q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)\alpha(z,y,z^{\prime},y^{\prime})+\bar{\alpha}_{\theta}(z,y)e^{sC(y,x)}$ $\displaystyle=\int\left(e^{sC(y^{\prime},x)}\alpha(z,y,z^{\prime},y^{\prime})+e^{sC(y,x)}-\alpha(z,y,z^{\prime},y^{\prime})e^{sC(y,x)}\right)q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)\;.$ $$

Since $C$ is bounded, we conclude that there exists a constant $b<\infty$ such that:

$$ $P_{\theta}V_{x}(z,y)\leq b\;.$ $$

For all $\left(z,y\right)\notin\mathcal{K}$,

$$ $\displaystyle\frac{P_{\theta}V_{x}(z,y)}{V_{x}(z,y)}$ $\displaystyle=\int\left(e^{s\left(C(y^{\prime},x)-C(y,x)\right)}\alpha(z,y,z^{\prime},y^{\prime})+1-\alpha(z,y,z^{\prime},y^{\prime})\right)q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)$ $\displaystyle\leq\int_{y^{\prime}\in\left\{C(y^{\prime},x)<C(y,x)\right\}}\left(e^{s\left(C(y^{\prime},x)-C(y,x)\right)}\alpha(z,y,z^{\prime},y^{\prime})+1-\alpha(z,y,z^{\prime},y^{\prime})\right)q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)$ $\displaystyle+\int_{y^{\prime}\in\left\{C(y^{\prime},x)=C(y,x)\right\}}\left(e^{s\left(C(y^{\prime},x)-C(y,x)\right)}\alpha(z,y,z^{\prime},y^{\prime})+1-\alpha(z,y,z^{\prime},y^{\prime})\right)q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)\;.$ $$

There exists $\eta\in(0,1)$ such that:

$$ $\displaystyle\int_{y^{\prime}\in\left\{C(y^{\prime},x)<C(y,x)\right\}}\left(e^{s\left(C(y^{\prime},x)-C(y,x)\right)}\alpha(z,y,z^{\prime},y^{\prime})+1-\alpha(z,y,z^{\prime},y^{\prime})\right)q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)$ $\displaystyle\leq\eta\int_{y^{\prime}\in\left\{C(y^{\prime},x)<C(y,x)\right\}}q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)\;.$ $$

Thus,

$$ $\displaystyle I(z,y)$ $\displaystyle\leq\eta+\left(1-\eta\right)\int_{y^{\prime}\in\left\{C(y^{\prime},x)=C(y,x)\right\}}q(\mathrm{d}z^{\prime}|z)p_{\theta}(\mathrm{d}y^{\prime}|z^{\prime},x)\;,$ $$

which establishes the drift condition and completes the proof of geometric ergodicity in total variation distance using Theorem [B.1](#A2.Thmtheorem1).

For $L_{2}$-geometric ergodicity, note that the Markov Chain is reversible with respect to the probability measure $\pi_{\theta}$ (Proposition [5.1](#S5.ThmProposition1)), so the Markov operator $P_{\theta}$ acts as a self-adjoint operator on $L_{2}$. The equivalence of geometric ergodicity and the existence of a spectral gap for $P_{\theta}$ acting on $L_{2}$ was established in (Roberts and Rosenthal, 1997, Theorem 2.1).
Specifically, it is shown that $P_{\theta}$ is $L_{2}$-geometrically ergodic if and only if it is $\pi_{\theta}$-TV geometrically ergodic. As a result, there exist constants $(\rho_{2},\kappa_{2})\in(0,1)\times\mathbb{R}_{+}$ such that for all $\mu\in\mathsf{M}_{1}(\mathsf{Z}\times\mathsf{Y})$, $\theta\in\Theta$, and $m\in\mathbb{N}$,

$$ $\|\mu P_{\theta}^{m}-\pi_{\theta}\|_{\mathrm{L}_{2}}\leq\kappa_{2}\rho_{2}^{m}\|\mu-\pi_{\theta}\|_{\mathrm{L}_{2}}\;.$ $$

∎

###### Proof.

#### C.3 Extension to K > 1 K>1

Here, we show how Theorem [5.1](#S5.Thmtheorem1) can be extended to the general case of $K>1$. Notably, at each iteration, $K$ chains are generated simultaneously and independently, with interactions occurring only among the particles from the previous iteration.

###### Proof.

The Markov kernel for general $K\geq 1$ is defined, for all $\mathsf{A}\in\mathcal{Z}^{K}\times\mathcal{Y}^{K}$, as

$$ $\mathsf{P}_{\theta}\big{(}(z^{1:K},y^{1:K}),\mathsf{A}\big{)}\\ =\int_{\mathsf{A}}\prod_{k=1}^{K}q(\mathrm{d}\tilde{z}^{k}|z)p_{\theta}(\mathrm{d}\tilde{y}^{k}|\tilde{z}^{k},x)\mathsf{\alpha}(z^{k},y^{k},\tilde{z}^{k},\tilde{y}^{k})+\bar{\alpha}_{\theta}(z^{k},y^{k})\delta_{(z^{k},y^{k})}(\mathrm{d}\tilde{z}^{k},\mathrm{d}\tilde{y}^{k})\;,$ $$

where $\alpha$ and $\bar{\alpha}_{\theta}$ are defined in ([8](#A3.E8)).

Using the same arguments as in the case $K=1$ (cf. Section [C.2](#A3.SS2)), we obtain the following minorization condition:

$$ $\mathsf{P}_{\theta}\big{(}(z^{1:K},y^{1:K}),A\big{)}\geq\varepsilon^{K}\nu^{\otimes K}(A)\;,$ $$

where $\varepsilon$ and $\nu$ are defined in ([11](#A3.E11)).
This completes the proof of the minorization condition.
We now turn to the drift condition: for a fixed $x\in\mathsf{X}$, we define the Lyapunov function for all $0<s\leq\lambda$ as

$$ $V^{K}_{x}(z^{1:K},y^{1:K})=\sum_{k=1}^{K}e^{sC(y^{k},x)}\;.$ $$

As in the case $K=1$, for all $\left(z^{1:K},y^{1:K}\right)\in\mathcal{K}^{K}$, applying $\mathsf{P}_{\theta}$ to $V^{K}_{x}$ yields:

$$ $\displaystyle\mathsf{P}_{\theta}V^{K}_{x}(z^{1:K},y^{1:K})$ $\displaystyle=\sum_{k=1}^{K}\int e^{sC(\tilde{y}^{k},x)}\left(q(\mathrm{d}\tilde{z}^{k}|z)p_{\theta}(\mathrm{d}\tilde{y}^{k}|\tilde{z}^{k},x)\alpha(z^{k},y^{k},\tilde{z}^{k},\tilde{y}^{k})+\bar{\alpha}_{\theta}(z^{k},y^{k})\delta_{(z^{k},y^{k})}(\mathrm{d}\tilde{z}^{k},\mathrm{d}\tilde{y}^{k})\right)$ $\displaystyle=\sum_{k=1}^{K}\int\left(\alpha(z^{k},y^{k},\tilde{z}^{k},\tilde{y}^{k})\left(e^{sC(\tilde{y}^{k},x)}-e^{sC(y^{k},x)}\right)+e^{sC(y^{k},x^{k})}\right)q(\mathrm{d}\tilde{z}^{k}|z)p_{\theta}(\mathrm{d}\tilde{y}^{k}|\tilde{z}^{k},x)\;.$ $$

Using the same bounding argument as in the $K=1$ case, we conclude that there exists a constant $b<\infty$ such that:

$$ $P_{\theta}V_{x}(z,y)\leq bK\;.$ $$

The case where the particles $(z^{1:K},y^{1:K})$ lie outside a compact set can be handled similarly, following the same reasoning as in the proof for $K=1$, thereby completing the extension to general $K$ for geometric ergodicity in total variation.

As in the case $K=1$, the reversibility of the $k$-th chain can be verified. Since the Markov kernel decomposes as a product, this structure ensures that reversibility holds component-wise, which in turn implies $L_{2}$-geometric ergodicity for $K\geq 1$.

∎

###### Proof.

### Appendix D Convergence Analysis for Adaptive θ \theta

#### D.1 Convergence Analysis of Time-Inhomogeneous MCMC algorithm with Stochastic Approximation Update

In this section, we present the general results of the time-inhomogeneous MCMC algorithm, where the objective is to sample from $\pi_{\theta}$, which itself depends on the parameter $\theta$ updated via SA. This setting corresponds to Algorithm [1](#alg1) with $K=1$, without imposing any assumptions on the proposal density $q$, in particular without requiring symmetry. We then apply these results to our proposed inference method.
We consider the following time-inhomogeneous MCMC algorithm with Stochastic Approximation update.

Figure: Algorithm 4 Time-Inhomogeneous MCMC with Stochastic Approximation Update

To analyze its convergence, we introduce the following assumptions.

###### Assumption 6

Let $(P_{k})_{k\geq 1}$ be a sequence of Markov transition kernels on $\left(\mathsf{X},\mathcal{X}\right)$. Suppose there exists a function $V:\mathsf{X}\to[1,\infty)$ satisfying $\sup_{x\in\mathsf{X}}V(x)<\infty$, and the following conditions hold:

- 1.
Minorization condition. There exist $\mathcal{K}\in\mathcal{X}$, $\varepsilon>0$ and a probability measure $\nu$ such that $\nu(\mathcal{K})>0$ and, for all
$A\in\mathcal{X}$ and $x\in\mathcal{K}$,
$P_{k}(x,A)\geq\varepsilon\nu(A)\;.$
- 2.
Drift condition. There exist constants $\lambda\in[0,1)$, $b\in(0,\infty)$ satisfying
$P_{k}V(x)\leq\begin{cases}\lambda V(x)&x\notin\mathcal{K}\;,\\
b&x\in\mathcal{K}\;.\end{cases}$

Assumption [6](#Thmassumption6) corresponds to a minorization and drift condition similar to those used in time-homogeneous MCMC, but it holds uniformly for the sequence of kernels.
While similar convergence guarantees can be established under weaker, non-uniform conditions (e.g., allowing the constants to depend on $k$), we focus on the uniform case as it aligns with our setting.

###### Assumption 7

There exists $L\in\mathsf{F}(\mathsf{X})$ such that for all $x\in\mathsf{X}$ and $\theta\in\Theta$,

$$ $\left\|\nabla_{\theta}\log\pi_{\theta}(x)\right\|\leq L(x)\;.$ $$

###### Assumption 8

There exists $\theta_{\infty}\in\Theta$ and a positive sequence $\left(a_{m}\right)_{m\in\mathbb{N}}$, with $a_{m}\to 0$ as $m\to\infty$ such that

$$ $\left\|\theta_{m}-\theta_{\infty}\right\|_{\mathrm{L}_{2}}^{2}=O(a_{m})\;.$ $$

Assumption [7](#Thmassumption7) is similar to the assumptions considered in (Andrieu and Atchadé, 2007; Andrieu and Moulines, 2006).
Notably, Assumption [8](#Thmassumption8) does not require $\theta_{\infty}$ to be a unique minimizer; it can simply be a critical point.

###### Theorem D.1

Let Assumptions [6](#Thmassumption6) - [8](#Thmassumption8) hold. Then, there exist a constant $\rho\in(0,1)$ and positive sequences $\left(b_{m}\right)_{m\in\mathbb{N}}$ such that for all $\mu\in\mathsf{M}_{1}(\mathsf{X})$, and $m\in\mathbb{N}$,

$$ $\displaystyle\mathbb{E}\Big{[}\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{\infty}}\right\|_{\mathrm{TV}}\Big{]}$ $\displaystyle=\mathcal{O}\bigg{(}\rho^{b_{m}}+\sum_{j=m-b_{m}}^{m-1}\gamma_{j+1}+a_{m}\bigg{)}\;.$ $$

Furthermore, if $\limsup_{m\to\infty}\left(b_{m}^{-1}+b_{m}/m+b_{m}\gamma_{m}\right)=0$, then

$$ $\mathbb{E}\Big{[}\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{\infty}}\right\|_{\mathrm{TV}}\Big{]}\xrightarrow[m\to\infty]{}0\;.$ $$

Theorem [D.1](#A4.Thmtheorem1) establishes that the iterates of the time-inhomogeneous MCMC with Stochastic Approximation step converge to the target distribution $\pi_{\theta_{\infty}}$. To establish this result, we first present a decomposition of the error in total variation in [D.1.1](#A4.SS1.SSS1), followed by an upper bound on the mixing error in [D.1.2](#A4.SS1.SSS2). The proof of the theorem is then provided in [D.1.3](#A4.SS1.SSS3).

###### Assumption 6

###### Assumption 7

###### Assumption 8

###### Theorem D.1

##### D.1.1 Error Decomposition

###### Lemma D.1

For all $1\leq s\leq m$, we have:

$$ $\displaystyle\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{\infty}}\right\|_{\mathrm{TV}}$ $\displaystyle\leq\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{s}}P_{\theta_{s+1}}\cdots P_{\theta_{m}}\right\|_{\mathrm{TV}}+\sum_{j=s}^{m-1}\left\|\pi_{\theta_{j+1}}-\pi_{\theta_{j}}\right\|_{\mathrm{TV}}$ $\displaystyle\quad+\left\|\pi_{\theta_{m}}-\pi_{\theta_{\infty}}\right\|_{\mathrm{TV}}\;.$ $$

###### Proof.

For all $m\in\mathbb{N}$, using the triangle inequality, we have:

$$ $\displaystyle\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{\infty}}\right\|_{\mathrm{TV}}$ $\displaystyle\leq\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{m}}\right\|_{\mathrm{TV}}+\left\|\pi_{\theta_{m}}-\pi_{\theta_{\infty}}\right\|_{\mathrm{TV}}\;.$ (12) $$

Using the fact that $P_{\theta_{m}}$ admits $\pi_{\theta_{m}}$ as an invariant measure and applying the triangle inequality, for all $1\leq s\leq m$, we have:

$$ $\displaystyle\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{m}}\right\|_{\mathrm{TV}}$ $\displaystyle\leq\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{s}}P_{\theta_{s+1}}\cdots P_{\theta_{m}}\right\|_{\mathrm{TV}}$ $\displaystyle\quad+\sum_{j=s}^{m-1}\left\|\pi_{\theta_{j}}P_{\theta_{j+1}}P_{\theta_{j+2}}\cdots P_{\theta_{m}}-\pi_{\theta_{j+1}}P_{\theta_{j+1}}P_{\theta_{j+2}}\cdots P_{\theta_{m}}\right\|_{\mathrm{TV}}$ $\displaystyle\leq\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{s}}P_{\theta_{s+1}}\cdots P_{\theta_{m}}\right\|_{\mathrm{TV}}+\sum_{j=s}^{m-1}\left\|\pi_{\theta_{j+1}}-\pi_{\theta_{j}}\right\|_{\mathrm{TV}}\;,$ $$

where the last inequality follows from the fact that, for all $j$, the Markov kernels $P_{\theta_{j}}$ are contractions. Together with ([12](#A4.E12)), this completes the proof.
∎

This bound decomposes the total variation error into three components: $(i)$ the mixing error (first term), which measures how well the Markov chain mixes from an arbitrary initial distribution; $(ii)$ the tracking error (second term), which measures how much the stationary distributions shift over time due to changes in the parameters; and $(iii)$ the optimization error (last term), which measures the difference between the current parameters and their limiting value $\theta_{\infty}$.

###### Lemma D.1

###### Proof.

##### D.1.2 Upper Bound on the Mixing Error

###### Proposition D.1

Let Assumption [6](#Thmassumption6) hold for a function $V:\mathsf{X}\to[1,\infty)$, where $\sup_{x\in\mathsf{X}}V(x)<\infty$. Let $(P_{k})_{k\geq 1}$ be a sequence of Markov transition kernels on $\left(\mathsf{X},\mathcal{X}\right)$. Then, there exists a constant $\rho\in(0,1)$ such that for all $\xi,\xi^{\prime}\in\mathsf{M}_{1}(\mathsf{X})$ and all $m\in\mathbb{N}$,

$$ $\displaystyle\|\xi P_{1}\cdots P_{m}-\xi^{\prime}P_{1}\cdots P_{m}\|_{\mathrm{TV}}$ $\displaystyle\leq\begin{cases}2\rho^{m}&\text{if }\xi,\xi^{\prime}\text{ are supported on }\mathcal{K}\;,\\ \rho^{m}\big{(}\xi(V)+\xi^{\prime}(V)\big{)}&\text{otherwise}\;.\end{cases}$ $$

Proposition [D.1](#A4.ThmProposition1) corresponds to the mixing rate of a time-homogeneous Markov Chain and is analogous to (Douc et al., 2004, Theorem 2), though it differs in both statement and proof. In general, there are various approaches to establish the geometric ergodicity of Markov Chains. Here, we follow a coupling argument. Specifically, we adapt the proof of homogeneous Markov Chains from (Douc et al., 2018, Theorem 19.4.1) to the time-homogeneous setting.

We construct a bivariate Markov Chain $(X_{k},X^{\prime}_{k})_{k\geq 1}$ such that, marginally, $(X_{k})_{k\geq 1}$ and $(X^{\prime}_{k})_{k\geq 1}$ are Markov Chains starting from $X_{1}=x$ and $X^{\prime}_{1}=x^{\prime}$, respectively, and each evolving according to the transition kernels $(P_{k})_{k\geq 1}$.

To achieve this, for all $k\geq 1$, we define the modified kernel $Q_{k}$, defined for all $A\subset\mathcal{X}$ and $x_{k}\in\mathsf{X}$ as:

$$ $Q_{k}\left(x_{k},A\right)=\frac{P_{k}\left(x_{k},A\right)-\varepsilon\nu\left(A\right)}{1-\varepsilon}\;,$ $$

and introduce the coupling kernel $\bar{P}_{k}$ on $\mathsf{X}^{2}$, defined for all $A\times A^{\prime}\subset\mathsf{X}^{2}$ and all $z_{k}=(x_{k},x_{k}^{\prime})\in\mathsf{X}^{2}$ by

$$ $\bar{P}_{k}\left(z_{k},A\times A^{\prime}\right)=\mathbf{1}_{x_{k}=x_{k}^{\prime}}P_{k}\left(x_{k},A\right)\delta_{x_{k+1}}\left(A^{\prime}\right)+\mathbf{1}_{x_{k}\neq x_{k}^{\prime}}\mathbf{1}_{z_{k}\notin\mathcal{K}^{2}}P_{k}\left(x_{k},A\right)P_{k}\left(x_{k}^{\prime},A^{\prime}\right)\\ +\mathbf{1}_{x_{k}\neq x_{k}^{\prime}}\mathbf{1}_{z_{k}\in\mathcal{K}^{2}}\left(\varepsilon\nu\left(A\right)\delta_{x_{k+1}}\left(A^{\prime}\right)+\left(1-\varepsilon\right)Q_{k}\left(x_{k},A\right)Q_{k}\left(x_{k}^{\prime},A^{\prime}\right)\right)\;.$ (13) $$

###### Lemma D.2

Let $\left(\bar{P}_{k}\right)_{k\geq 0}$ be the Markov kernels on $\left(\mathsf{X}^{2},\mathcal{X}^{2}\right)$ defined by ([13](#A4.E13)). Then, for all $n\in\mathbb{N}$ and all $(x,x^{\prime})\in\mathsf{X}^{2}$, we have

$$ $\bar{P}_{1}\cdots\bar{P}_{m}\left(\left(x,x^{\prime}\right),\cdot\right)\in\mathcal{C}\left(P_{1}\cdots P_{m}(x,\cdot),P_{1}\cdots P_{m}\left(x^{\prime},\cdot\right)\right)\;,$ $$

where $\mathcal{C}$ denotes the set of couplings introduced in Definition [B.3](#A2.Thmdefinition3).

###### Proof.

We proceed by induction on $m$. By the definition of $\bar{P}_{1}$, we have by construction

$$ $\bar{P}_{1}\left(\left(x,x^{\prime}\right),\cdot\right)\in\mathcal{C}\left(P_{1}(x,\cdot),P_{1}(x^{\prime},\cdot)\right)\;.$ $$

Suppose that for some $m\geq 1$, we have

$$ $\bar{P}_{1}\cdots\bar{P}_{m}\left(\left(x,x^{\prime}\right),\cdot\right)\in\mathcal{C}\left(P_{1}\cdots P_{m}(x,\cdot),P_{1}\cdots P_{m}\left(x^{\prime},\cdot\right)\right)\;.$ $$

Applying $\bar{P}_{m+1}$ to both sides and using the definition of composition of Markov kernels, we obtain, by definition of $\bar{P}_{m+1}$,

$$ $\displaystyle\bar{P}_{1}\cdots\bar{P}_{m+1}\left(\left(x,x^{\prime}\right),A\times\mathsf{X}\right)$ $\displaystyle=\int_{\mathsf{X}\times\mathsf{X}}\bar{P}_{1}\cdots\bar{P}_{m}\left((x,x^{\prime}),\mathrm{d}y\mathrm{d}y^{\prime}\right)\bar{P}_{m+1}\left((y,y^{\prime}),A\times\mathsf{X}\right)$ $\displaystyle=\int_{\mathsf{X}\times\mathsf{X}}\bar{P}_{1}\cdots\bar{P}_{m}\left((x,x^{\prime}),\mathrm{d}y\mathrm{d}y^{\prime}\right)P_{m+1}\left(y,A\right)$ $\displaystyle=\int_{\mathsf{X}\times\mathsf{X}}P_{1}\cdots P_{m}\left(x,\mathrm{d}y\right)P_{m+1}\left(y,A\right)$ $\displaystyle=P_{1}\cdots P_{m+1}\left(x,A\right)\;,$ $$

where the third equality follows from the inductive hypothesis.
Similarly, we obtain

$$ $\bar{P}_{1}\cdots\bar{P}_{m+1}\left(\left(x,x^{\prime}\right),\mathsf{X}\times A\right)=P_{1}\cdots P_{m+1}\left(x^{\prime},A\right)\;.$ $$

Thus, we conclude that

$$ $\bar{P}_{1}\cdots\bar{P}_{m+1}\left(\left(x,x^{\prime}\right),\cdot\right)\in\mathcal{C}\left(P_{1}\cdots P_{m+1}(x,\cdot),P_{1}\cdots P_{m+1}\left(x^{\prime},\cdot\right)\right)\;.$ $$

This completes the induction and proof.
∎

###### Lemma D.3

For all $(x,x^{\prime})\in\mathsf{X}^{2}$ define $\Delta(x,x^{\prime})=\mathbf{1}_{x\neq x^{\prime}}$ and $\bar{V}(x,x^{\prime})=\left(V(x)+V(x^{\prime})\right)/2$. Then, for all $k\in\mathbb{N}$:

- •
If $(x,x^{\prime})\in\mathcal{K}^{2}$, then
$\bar{P}_{k}\Delta(x,x^{\prime})\leq(1-\varepsilon)\Delta(x,x^{\prime}),\quad\bar{P}_{k}\bar{V}(x,x^{\prime})\leq b\;.$
- •
If $(x,x^{\prime})\notin\mathcal{K}^{2}$, then
$\bar{P}_{k}\Delta(x,x^{\prime})\leq\Delta(x,x^{\prime}),\quad\bar{P}_{k}\bar{V}(x,x^{\prime})\leq\lambda\bar{V}(x,x^{\prime})\;.$

###### Proof.

The inequality for $\bar{P}_{k}\Delta$ in both cases follows immediately from the definition of $\bar{P}_{k}$.
For the second inequality, we have:

$$ $\bar{P}_{k}\bar{V}(x,x^{\prime})=\frac{P_{k}V(x)+P_{k}V(x^{\prime})}{2}\;.$ $$

If $(x,x^{\prime})\in\mathcal{K}^{2}$, then since $P_{k}V(x)\leq b$ and $P_{k}V(x^{\prime})\leq b$, it follows that:

$$ $\frac{P_{k}V(x)+P_{k}V(x^{\prime})}{2}\leq b\;.$ $$

If $(x,x^{\prime})\notin\mathcal{K}^{2}$, using $P_{k}V(x)\leq\lambda V(x)$ and $P_{k}V(x^{\prime})\leq\lambda V(x^{\prime})$, we obtain:

$$ $\bar{P}_{k}\bar{V}(x,x^{\prime})=\frac{P_{k}V(x)+P_{k}V(x^{\prime})}{2}\leq\frac{\lambda V(x)+\lambda V(x^{\prime})}{2}=\lambda\bar{V}(x,x^{\prime})\;.$ $$

This concludes the proof.
∎

###### Proof of Proposition D.1 .

For any $t\in(0,1)$, we define

$$ $\varrho_{t}=\max\left(\left(1-\varepsilon\right)^{1-t}b^{t},\lambda^{t}\right)\;.$ $$

For chosen $t$, we introduce the function:
$W(x,x^{\prime})=\Delta^{1-t}1_{(x,x^{\prime})\in\mathcal{K}^{2}}+\Delta^{1-t}\bar{V}^{t}1_{(x,x^{\prime})\notin\mathcal{K}^{2}}$. Then, using Lemma [D.2](#A4.Thmlemma2), we have:

$$ $\displaystyle\left\|P_{1}\cdots P_{m}(x,\cdot)-P_{1}\cdots P_{m}\left(x^{\prime},\cdot\right)\right\|_{\mathrm{TV}}$ $\displaystyle=2\inf\left\{\zeta(\Delta):\zeta\in\mathcal{C}(P_{1}\cdots P_{m}(x,\cdot),P_{1}\cdots P_{m}\left(x^{\prime},\cdot\right))\right\}$ $\displaystyle\leq 2\bar{P}_{1}\cdots\bar{P}_{m}\Delta\left(x,x^{\prime}\right)$ $\displaystyle\leq 2\bar{P}_{1}\cdots\bar{P}_{m}W(x,x^{\prime})\;.$ $$

where we used $V\geq 1$.
Finally, applying Hölder’s inequality and using Lemma [D.3](#A4.Thmlemma3), we obtain, for all $\left(x,x^{\prime}\right)\in\mathsf{X}^{2}$,

$$ $\displaystyle\bar{P}_{k}W\left(x,x^{\prime}\right)=\bar{P}_{k}\left(\Delta^{1-t}\bar{V}^{t}\right)\left(x,x^{\prime}\right)$ $\displaystyle\leq\left(\bar{P}_{k}\Delta\left(x,x^{\prime}\right)\right)^{1-t}\left(\bar{P}_{k}\bar{V}\left(x,x^{\prime}\right)\right)^{t}$ $\displaystyle\leq\Delta^{1-t}\left(x,x^{\prime}\right)\times$ $\displaystyle\leq\varrho_{t}W\left(x,x^{\prime}\right)\;.$ $$

This implies by induction that for all $m\in\mathbb{N}$ and all $\left(x,x^{\prime}\right)\in\mathsf{X}^{2}$,

$$ $\bar{P}_{1}\cdots\bar{P}_{m}W\left(x,x^{\prime}\right)\leq\varrho_{t}^{m}W\left(x,x^{\prime}\right)\;.$ $$

Then,

$$ $\displaystyle\left\|P_{1}\cdots P_{m}(x,\cdot)-P_{1}\cdots P_{m}\left(x^{\prime},\cdot\right)\right\|_{\mathrm{TV}}$ $\displaystyle\leq 2\varrho_{t}^{m}W\left(x,x^{\prime}\right)$ $\displaystyle\leq\begin{cases}2\varrho_{t}^{m}&\text{ if }x\in\mathcal{K}\;,\\ \varrho_{t}^{m}\left(V(x)+V\left(x^{\prime}\right)\right)&\text{ if }x\notin\mathcal{K}\;.\end{cases}$ $$

This concludes the proof.
∎

###### Proposition D.1

###### Lemma D.2

###### Proof.

###### Lemma D.3

###### Proof.

###### Proof of Proposition D.1 .

##### D.1.3 Proof of Theorem D.1

###### Proof.

Using Lemma [D.1](#A4.Thmlemma1), and taking $s=m-b_{m}$, we have:

$$ $\displaystyle\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{\infty}}\right\|_{\mathrm{TV}}$ $\displaystyle\leq\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{m-b_{m}}}P_{\theta_{m-b_{m}}}\cdots P_{\theta_{m}}\right\|_{\mathrm{TV}}$ $\displaystyle\quad+\sum_{j=m-b_{m}}^{m-1}\left\|\pi_{\theta_{j+1}}-\pi_{\theta_{j}}\right\|_{\mathrm{TV}}+\left\|\pi_{\theta_{m}}-\pi_{\theta_{\infty}}\right\|_{\mathrm{TV}}$ $\displaystyle\leq\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{m-b_{m}}}P_{\theta_{m-b_{m}}}\cdots P_{\theta_{m}}\right\|_{\mathrm{TV}}$ $\displaystyle\quad+L(x)\sum_{j=m-b_{m}}^{m-1}\left\|\theta_{j+1}-\theta_{j}\right\|+L(x)\left\|\theta_{m}-\theta_{\infty}\right\|\;,$ $$

where we used the Lipschitz condition of $\pi_{\theta}$.
For the first term, using Proposition [D.1](#A4.ThmProposition1) with $\xi=\mu P_{\theta_{1}}\cdots P_{\theta_{m-b_{m}-1}}$ and $\xi^{\prime}=\pi_{\theta_{m-b_{m}}}$, we have:

$$ $\left\|\mu P_{\theta_{1}}\cdots P_{\theta_{m}}-\pi_{\theta_{m-b_{m}}}P_{\theta_{m-b_{m}}}\cdots P_{\theta_{m}}\right\|_{\mathrm{TV}}\leq\kappa\rho^{b_{m}}\;.$ $$

For the second term, using the Lipschitz condition (Assumption [7](#Thmassumption7)) and the recursion of $\theta_{j+1}$, we get:

$$ $\displaystyle\sum_{j=m-b_{m}}^{m-1}\mathbb{E}\left[\left\|\theta_{j+1}-\theta_{j}\right\|\right]$ $\displaystyle=\sum_{j=m-b_{m}}^{m-1}\gamma_{j+1}\mathbb{E}\left[\left\|H_{\theta_{m}}\left(x_{m+1}\right)\right\|\right]$ $\displaystyle=\mathbb{E}\left[L(x)\right]\sum_{j=m-b_{m}}^{m-1}\gamma_{j+1}\;.$ $$

For the last term, Using Jensen inequality and Assumption [8](#Thmassumption8), we obtain:

$$ $\mathbb{E}\left[\left\|\theta_{m}-\theta_{\infty}\right\|\right]\leq\left\|\theta_{m}-\theta_{\infty}\right\|_{\mathrm{L}_{2}}^{2}=\mathcal{O}\left(a_{m}\right)\;.$ $$

∎

###### Proof.

#### D.2 Proof of Theorem 5.2

For a fixed $x\in\mathsf{X}$, we define the Lyapunov function for all $0<s\leq\lambda$ as

$$ $V_{x}(z,y)=1+e^{sC(y,x)}\;.$ $$

Following the procedure outlined in the proof of Theorem [5.1](#S5.Thmtheorem1), it is straightforward to verify Assumption [6](#Thmassumption6) (the minoration and drift condition) with $V\geq 1$.
Additionally, using Assumptions [3](#Thmassumption3) and [4](#Thmassumption4), we can verify Assumptions [7](#Thmassumption7) and [8](#Thmassumption8). The proof is then concluded by applying Theorem [D.1](#A4.Thmtheorem1).

#### D.3 Convergence Rate in Stochastic Approximation

Stochastic Approximation can be traced back to Robbins and Monro (1951). Since then, numerous variants have been proposed, including those using adaptive step sizes, such as Kingma and Ba (2015).
The non-asymptotic convergence of biased SA methods has been studied in various settings. For instance, Karimi et al. (2019) analyzes the case without adaptive step sizes, while Surendran et al. (2024) extends the analysis to include adaptive schemes for non-convex smooth objectives.
These works provide convergence guarantees in terms of the squared norm of the gradient of the objective function. Specifically, they show that the iterates converge to a critical point at a rate of $\mathcal{O}(\log m/\sqrt{m}+b)$, where $b$ corresponds to the bias and $m$ to the number of iterations.
The analysis typically relies on standard assumptions, including the smoothness of the objective function, an assumption on the bias and variance of the gradient estimator (see Assumption H3 in Surendran et al. (2024)), and a decreasing step size.

In our setting, given that the cost is bounded, the smoothness of the test objective $\mathcal{L}_{test}$ hinges on the smoothness of the policy $p_{\theta}$, an assumption also used in Papini et al. (2018); Surendran et al. (2025). The stochastic update defined in ([4](#S4.E4)) is bounded under Assumption [3](#Thmassumption3), which allows us to verify the necessary conditions on the bias and variance. In particular, the bias is of order $\mathcal{O}(1/K)$. Therefore, with an additional smoothness assumption on $p_{\theta}$ and a suitable choice of step sizes, such as $\gamma_{m}=1/\sqrt{m}$, Assumption [4](#Thmassumption4) can be satisfied.

### Appendix E Additional Experiments

#### E.1 Training Details

In this section, we provide the details of our model and training procedure. The encoder uses multi-head attention with 8 heads and an embedding dimension of $d_{h}=128$, and consists of 6 layers. The decoder includes a single multi-head attention layer with 8 heads and a key dimension of $d_{k}=16$.

For both the TSP and CVRP, node coordinates $c_{i}$ are sampled uniformly within the unit square. In CVRP instances, customer demands $d_{i}$ are drawn from a uniform distribution $\mathcal{U}([1,10])$. The vehicle capacity $D$ is set based on the number of nodes: $D=50$ for $n=100$, $D=55$ for $n=125$, and $D=60$ for $n=150$, following the setup used in the literature (Hottung et al., 2021).

Training is conducted only for instances with $n=100$ nodes, using $K=100$ latent samples. The latent space is defined as a compact space with diameter $R=40$ and dimension $d_{z}=100$. We use the Adam optimizer with a learning rate of $5\times 10^{-4}$, a batch size of $128$, and train for $2000$ epochs. The momentum parameters are fixed at $\beta_{1}=0.9$ and $\beta_{2}=0.999$, with a weight decay of $1\times 10^{-6}$. The entropic regularization parameter $\beta$ is set to $0.01$, and the weights $\tau$ in the loss ([2](#S4.E2)) are chosen according to an exponential decay schedule.
The training loss and corresponding cost are shown in Figure [4](#A5.F4).

Figure: Figure 4: Training loss and score for our model trained on CVRP instances with nodes $n=100$
Refer to caption: /html/2506.03672/assets/x4.png

#### E.2 Inference Details and Additional Illustrations

##### E.2.1 Inference Details

The inference results typically include the objective value (cost), the optimality gap, and the computation time. For example, in Tables [1](#S6.T1) and [2](#S6.T2), Obj. denotes the value of the cost function defined in ([5](#A1.E5)) for the TSP and CVRP. The Gap indicates the percentage gap to optimality, computed as:

$$ $Gap(y,y^{*})=\left(\frac{C(y,x)}{C(y^{*},x)}-1\right)*100\%$ $$

where $y^{*}$ denotes the optimal solution for TSP and the near-optimal solution provided by LKH-3 for CVRP.

In our experiments, we use a batch size of 200 for TSP and 100 for CVRP with $K=600$ latent samples when the augmentation trick is not applied. When using the augmentation trick, we reduce the batch size to 100 for TSP and 50 for CVRP, and $K=300$ latent samples.

The proposal distribution is a Gaussian with density: for all $m\in\mathbb{N}$ and $1\leq k\leq K$,

$$ $q(z^{k}_{m+1}\mid z^{1:K}_{m})=\mathcal{N}\left(z_{m+1}^{k};z_{m}^{k}+\gamma\left(z^{I_{1}}_{m}-z^{I_{2}}_{m}\right),\sigma^{2}I_{d_{z}}\right)\;,$ $$

where $I_{1},I_{2}\sim\mathcal{U}(\{1,\dots,K\})$. The variance parameter is set to $\sigma^{2}=0.01$ and the scaling factor is $\gamma=0.319$ for TSP and $\gamma=0.379$ for CVRP.
Instead of updating the parameters at every iteration, updates are performed at fixed intervals. The update schedule is described in the next section and illustrated in Figure [7](#A5.F7).

##### E.2.2 Illustration of hyperparameters

Figure [5](#A5.F5) illustrates solutions generated by our method, following a similar visualization style as in Perron and Furnon (2019); Kool et al. (2019). Visually, the solutions appear to be optimal.

Figure: Figure 5: Solution representation produced by our model on the TSP (left) and CVRP (right) with $n=50$. In the TSP plot, the red point denotes the starting node, and arrows indicate the visiting order. In the CVRP plot, the large square represents the depot, and each color corresponds to a distinct vehicle route. The bar illustrates vehicle capacity usage: black segments show the portion used by each customer, while white segments indicate unused capacity. The overall height of each bar reflects the total load on the corresponding route.
Refer to caption: /html/2506.03672/assets/x5.png

To highlight the impact of important hyperparameters, we first focus on the number of latent samples $K$. We observe that increasing $K$ improves the results, as stated in Theorem [5.2](#S5.Thmtheorem2), particularly due to the constant arising from the minoration condition. However, beyond a certain threshold, the improvement becomes insignificant. Choosing an appropriate value of $K$ is crucial to balance faster mixing with computational cost.

Figure: Figure 6: Average gap to optimality for different values of $K$ in the CVRP with $n=100$
Refer to caption: /html/2506.03672/assets/x7.png

Next, we discuss the SA step, focusing on how frequently the parameters should be updated. Our initial motivation for not updating the parameters at every iteration is twofold: to reduce computational cost and to better explore the latent space given the current parameter distribution. Consequently, parameter updates can be performed at regular intervals rather than every iteration.

Since the initial parameters are typically far from optimal, allowing long exploration intervals early in training is often unnecessary. Instead, we begin with short exploration intervals and gradually increase them to enable more thorough exploration as learning progresses. Selecting these intervals is non-trivial; we chose them manually without extensive hyperparameter tuning. The update schedule used in our experiments is $[1,1,5,15,25,100,150]$. Optimizing this schedule remains an open question and presents an interesting direction for future work.

Figure: Figure 7: Performance comparison of various SA step update frequencies on CVRP with $n=100$
Refer to caption: /html/2506.03672/assets/x8.png

Figure [7](#A5.F7) illustrates the average gap to optimality for different choices of parameter update frequency, including both regular intervals and the increasing schedule. Here, $M_{0}$ denotes the update frequency, with $M_{0}=50$ indicating that parameters are updated every 50 iterations. We observe that while regular intervals produce similar results overall, $M_{0}=75$ yields slightly better performance. Notably, the increasing update schedule achieves a clear improvement over the fixed schedules, highlighting the potential benefits of adaptive strategies.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 4
questions: What is the sensitivity to the choice of latent dimension d_z and bounded diameter R? The ablations only vary K and update schedules, not d_z or R.
soundness: 3
strengths: The instance-conditioned latent model eliminates the need for pre-computed solutions or pre-trained policies unlike CVAE-Opt and COMPASS, and demonstrates clear performance improvements. LGS inference combines interacting MCMC with real-time parameter updates, consistently outperforming all alternatives including DE, CMA-ES, active search, and gradient-based finetuning. The convergence proof for time-inhomogeneous Markov chains is a non-trivial theoretical contribution for this problem class. The experiments include multiple datasets, comparisons with and without augmentation, runtime reports, and detailed hyperparameters, demonstrating excellent reproducibility.
summary: LGS-Net is an RL-based latent variable neural combinatorial optimization model that requires no labels. It uses instance-conditioned latent embeddings and trains via a cost-weighted, entropy-regularized objective function. The key innovation is the Latent Guided Sampling (LGS) inference procedure, which updates decoder parameters through stochastic approximation while running interacting MCMC chains in latent space. It provides convergence guarantees for both fixed and adaptive parameters. The model achieves learning-based SOTA on TSP and CVRP benchmarks, and even outperforms LKH3 on some CVRP settings.
weaknesses: The adaptive chain convergence relies on Assumption 4, but the paper only provides high-level justification. There is a lack of sufficient conditions linking Algorithm 1's specific step-size choices and gradient-variance conditions, creating a gap between theory and implementation. 
The empirical study is limited to Euclidean routing (TSP/CVRP), restricting the demonstration of generality. Experiments on non-Euclidean problems or other combinatorial domains would help establish the versatility of the latent formulation. 
Although negative gaps versus LKH3 on CVRP are reported, verification via exact solvers or optimality certificates is absent, making the evaluation protocol ambiguous.

---

confidence: 4
contribution: 2
initial_rating: 4
presentation: 2
questions: It would be valuable to provide stronger motivation for the design choices introduced in the proposed method. For instance, why was the POPPY loss from COMPASS removed? Why was a latent-space encoder introduced? And why are all decoder parameters updated if the learned latent space is already expected to capture much of the variability? I can anticipate some of the rationale behind these decisions, but I would really appreciate hearing the authors’ explicit perspective. Clarifying these points, ideally supported by additional ablation studies, could substantially strengthen the manuscript.

In Section 4.1, it is also worth noting that the baseline AM/POMO architecture already includes an encoder. It would therefore help to clarify what the new encoder adds, and what type of information is expected to be captured in the latent space that is not already represented in the existing embedding.

It would also be helpful to report the dimensionality of the latent space used in the experiments.

Regarding the reported timings: could the authors confirm whether these values are accurate? Given the inconsistencies mentioned earlier, I am concerned about the validity of the experimental results.

Finally, on line 116, the paper states that LGS “removes the need for a pre-trained policy.” Could the authors elaborate on why this is the case?

Overall, I find the work promising and I am genuinely keen to support it. However, several concerns remain that need to be addressed to lean toward acceptance. I sincerely hope that the authors will be able to provide clarifications and additional results during the rebuttal phase to help advance this discussion.
soundness: 2
strengths: Significance: The use of inference-time strategies for combinatorial optimization (CO) problems is highly relevant. The proposed method achieves state-of-the-art performance with statistically significant improvements on reference benchmarks, evaluated on two well-studied problems (TSP and CVRP), both in- and out-of-distribution. I appreciate that the paper provides not only SOTA results but also mathematical justification for the proposed approach.

Originality: Most underlying concepts already exist, and the general idea of latent-space search is not new (which the authors acknowledge clearly), but the combination of latent-space modeling, MCMC sampling, and gradient-based updates is novel and supported by solid mathematical analysis. While there is no major paradigm shift, the proposed approach presents an interesting and meaningful contribution.

Clarity: The paper makes a commendable effort to explain the method, and it presents the mathematical components and proofs clearly. The main ideas are easy to follow, and the inclusion of figures, algorithms, and plots effectively supports the reader’s understanding.

Quality: The paper engages well with the existing literature, being transparent about its close connections to COMPASS, EAS, and CVAE-Opt. It presents the mathematical aspects rigorously and includes results and extended materials in the appendix. Overall, the paper demonstrates good quality, though there remains room for improvement (see weaknesses).
summary: This paper introduces a new inference strategy for neural solvers applied to combinatorial optimization, with a primary focus on routing problems but applicability beyond them. The approach is based on latent-space models coupled with a subsequent latent-space search. It includes a training phase that anticipates the inference-time search budget and presents a coupled training–inference procedure. The method establishes close connections with prior inference strategies in the literature, namely EAS, COMPASS, and CVAE-Opt. The paper provides both mathematical guarantees and empirical results on two widely studied combinatorial optimization problems (TSP & CVRP) evaluated in- and out-of-distribution.
weaknesses: **W1. Weak motivation and unclear justification of contributions.**

The main motivation of the work remains somewhat weak. Lines 51–54 and 57–60 present the limitations of prior work (EAS and COMPASS) that LGS is supposed to address (lines 60–62). This section is crucial, as it motivates the entire paper, yet it is not fully convincing.
(i) The authors claim that EAS “fine-tunes” the policy, which poses computational challenges, but LGS also performs fine-tuning.
 (ii) The statement that “COMPASS enforces independence between the problem instance and the latent space structure” would benefit from clarification, specifically, how limiting this assumption truly is, and why relaxing it is expected to yield substantial performance gains. 

Moreover, the results and ablation studies do not provide evidence that the “instance dependence” introduced in LGS is responsible for the observed performance improvements. The only ablation reported in the main paper concerns the inference-time search/sampling component. One could therefore wonder whether applying the LGS sampling and parameter-update mechanism to a COMPASS checkpoint would suffice. Overall, the paper does not convincingly establish that the proposed modifications address the stated weaknesses of prior methods. Clearer intuitions about why the proposed approach constitutes a conceptual improvement, and corresponding empirical validation through targeted ablations, would strengthen the work considerably.

Another motivation mentioned in the introduction is the removal of reliance on augmentation tricks. However, it is not clear how LGS achieves this, as other methods such as COMPASS and EAS do not fundamentally rely on augmentation either. They may use it, but they function well without it. The authors should clarify what makes LGS different in this regard; if there is no clear distinction, this point should not be used as a central motivation.

**W2. Missing ablation studies.**

Related to the previous point, some important ablations are missing. The only ablation presented concerns the sampling method, whereas one of the main motivations, conditioning on the instance, is not ablated. Ideally, this should be tested explicitly. Alternatively, the authors could revise their narrative to avoid presenting this aspect as a central feature of the method without empirical support.

**W3. Questionable timing results.**

The timing results reported in Tables 1, 4, and 5 appear inconsistent, which raises concerns about the validity of the experimental study. The times reported for POMO, EAS, and COMPASS do not align with expectations, and the reported runtime for LGS is also surprising.

First, the times for POMO and COMPASS should be roughly identical, as both share the same bottlenecks, the forward pass and environment rollouts. Their only difference is the CMA-ES update, which is negligible (see the original COMPASS paper). Therefore, one should expect similar runtimes across the benchmark, which is not the case (e.g., 10M vs. 20M on TSP100; 1h30 vs. 40M on TSP150).

Second, the reported time comparison between POMO and the other methods lacks coherence, as POMO is sometimes faster and sometimes slower, although it should consistently be equal or faster.

Third, EAS and COMPASS are reported to take the same time, which is implausible, since EAS involves backpropagation steps that should make it slower.

Finally, LGS-Net is reported to have the same runtime as COMPASS, even though it includes regular parameter updates. It should therefore be slower (perhaps only slightly if updates are infrequent), but this needs discussion if the times are indeed identical.

Overall, I am confident that there are issues with the reported timings. The authors should carefully verify and correct these results, and provide more detailed information on the time performance of their algorithm. Beyond hindering fair comparison between COMPASS and LGS-Net, these inconsistencies undermine confidence in the experimental section as a whole.

**W4. Clarity and presentation.**

The clarity of the paper could be improved. Figure 1 could include more detail to illustrate key components. It is somewhat confusing how the encoding interacts with the existing encoding that AM and POMO architectures are already using; the authors should be explicit about what differs in LGS, and how the latent space relates (or not) to the original embedding space. Additionally, certain elements, such as the “proposal distribution” (used in Algorithm 1 but not discussed in Section 4), should be better introduced.

Minor suggestions:
- “To learn solution strategies” (line 14) sounds awkward.
- It is acceptable to use “NCO” without specifying “learning-based NCO,” as the term already implies learning; this would simplify some sentences (e.g., lines 24 and 39).
- Line 52: stating that EAS is a SOTA search method is inaccurate, as it is outperformed on several benchmarks by COMPASS, PolyNet, and MEMENTO. The sentence should be softened to “among the leading methods for search-based RL.” MEMENTO [1] (NeurIPS 2025, Spotlight) could also be added to the related work.
- Line 62: “rely on the augmentation trick” is misleading, most methods are agnostic to it, even if some experiments make use of it.
- Line 113: COMPASS does not learn the latent space; it uses a fixed prior, and the policy is trained to exploit it effectively.
- Line 116: “removes the need for a pre-trained policy” should be clarified, why does this hold?
- Line 251. Small suggestion to swap theta and phi to follow the order in which encoder and decoder are mentioned.
- Line 382. I believe the usual dataset used in the NCO literature for TSP100 has 10000 instances, and TSP125 and TSP150 have 1000 instances. 

[1] Memory-Enhanced Neural Solvers for Routing Problems, Neurips 2025

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. How would you think about applying such a training and inference scheme outside of this NCO setup?

2. How realistic with respect to real problems are the assumptions made in the convergence analysis?
soundness: 4
strengths: The paper is overall well-written, gives good credit to related works, and performs decent experiments and ablations.
- Strong theoretical grounding regarding convergence.
- Standard experiments on TSP and CVRP: n in [100, 125, 150].
- SOTA results on these benchmarks, even beaten the solver for CVRP 100-125.
summary: The paper introduces LGS-Net, a latent-space neural model for combinatorial optimization (e.g., TSP, CVRP). Similar to [Chalumeau et al., 2023], it learns instance-conditioned latent representations using reinforcement learning, yet performs inference through a new scheme based on MCMC called Latent Guided Sampling (LGS).

Contributions:
- A new latent-space NCO model (LGS-Net) that conditions on problem instances.
- A provably convergent inference method (LGS) integrating MCMC and stochastic approximation (SA).
- SOTA results on TSP and CVRP benchmarks (in distribution and slightly out of distribution), outperforming prior methods like POMO, COMPASS, and EAS, with strong generalization to unseen problem sizes.
weaknesses: The paper and method suffer from a few weaknesses.
- Additional complexity of the inference scheme: MCMC and SA must bring a significant overhead.
- Limited experimental scope: I assume this would transfer well to other combinatorial problems like job-shop scheduling or graph problems, but the paper only tests on TSP and CVRP.

---

confidence: 4
contribution: 1
initial_rating: 4
presentation: 4
questions: - Why is the latent variable needed?
- Why not use only a decoder model that takes as input the problem instance x instead of the latent z(x)?
- Why is the latent variable stochastic and why couldn't a deterministic representation work as well?
- Can you show the benefit of having an encoder and the stochastic variable?
soundness: 3
strengths: - The paper addresses recognized limitations of prior NCO approaches (e.g., data requirements, lack of inference-time optimization) and positions the work within a clear conceptual framework.

- Principled combination of latent modeling and RL. The unification of stochastic approximation with MCMC inference in a learned latent space is elegant and theoretically grounded.

- Theoretical soundness. The convergence proofs are rigorous and build on established MCMC theory, adapted to a time-inhomogeneous setting.

- Comprehensive experimental comparison. The authors benchmark against numerous recent NCO baselines and demonstrate consistent improvements, confirming the competitiveness of the approach.
summary: The paper proposes LGS-Net, a latent-space model for Neural Combinatorial Optimization (NCO). The approach trains an encoder to map a problem instance into a distribution for a latent variable.  Based on a sample of this distribution, the decoder generates a distribution over the solution space which then is sampled to obtain solutions. At inference time, the authors introduce Latent Guided Sampling (LGS) — an MCMC-based method augmented with stochastic approximation updates to the decoder parameters, allowing test-time adaptation. Theoretical analysis establishes convergence guarantees for both fixed and adaptive parameter settings. Experiments on TSP and CVRP benchmarks show that LGS-Net achieves state-of-the-art or near-optimal performance, sometimes improving upon strong baselines such as COMPASS and EAS by small margins.
weaknesses: - Unclear role of the concept of the latent space: During inference this method employs a local search based on MCMC updates and weight updates of the decoder. This iterative search evaluates the solution quality of intermediate samples in solution space. Hence this method can hardly be considered a latent CO method since it does not learn to solve CO problems via latent space dynamics but relies on generating samples in solution space and evaluating the cost function for them. More generally, the role of the stochastic latent variable is not clear.

- Incremental novelty. The method extends earlier latent-space NCO works such as CVAE-Opt and COMPASS. The architectural and training components are similar, with the main novelty residing in the inference mechanism.

- Limited empirical diversity. All experiments are on Euclidean routing problems (TSP/CVRP ≤ 150 nodes). The approach’s generality to other combinatorial settings (e.g., knapsack, scheduling) remains untested.

- Small empirical margins. Performance gains over prior SOTA are modest (typically < 0.2 % gap). It is unclear whether these differences are statistically significant and practically meaningful given the added complexity of MCMC + SA inference.


Decision: Reject