# Latent Guided Sampling for Combinatorial Optimization

Sobihan Surendran*†, Adeline Fermanian†, and Sylvain Le Corff*

*Sorbonne Université and Université Paris Cité,
CNRS, Laboratoire de Probabilités, Statistique et Modélisation,
F-75005 Paris, France

†LOPF, Califrais’ Machine Learning Lab, Paris, France

## 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 $(\mathcal{X}, \mathcal{X})$, where $\mathcal{X}$ is a countably generated $\sigma$-algebra, let $F(\mathcal{X})$ denote the set of all measurable functions defined on $(\mathcal{X}, \mathcal{X})$. Let $M(\mathcal{X})$ be the set of $\sigma$-finite measures on $(\mathcal{X}, \mathcal{X})$, and $M_1(\mathcal{X}) \subset M(\mathcal{X})$ the probability measures. For all $f \in F(\mathcal{X})$ and $\mu \in M(\mathcal{X})$, we write $\mu(f) = \int f(x)\mu(dx)$. For a Markov kernel $P$ on $(\mathcal{X}, \mathcal{X})$ and $\mu \in M_1(\mathcal{X})$, the composition $\mu P$ is defined as $\mu P: \mathcal{X} \ni A \mapsto \int \mu(dx)P(x, dy)1_A(y)$. For probability measures $\mu$ and $\nu$ defined on the same measurable space, the Total Variation (TV) is defined as $\|\mu - \nu\|_{\text{TV}} := \sup_{A \in \mathcal{X}} |\mu(A) - \nu(A)|$. The $L_2$-norm of a random variable $X$ is defined as $\|X\|_{L_2} := (\mathbb{E}[|X|^2])^{1/2}$. The Hadamard product of two vectors $u$ and $v$ is denoted by $u \odot v$. For a sequence $(a_m)_{m \in \mathbb{N}}$, and all $u \le v$, we write $a_{u:v} = \{a_u, \dots, a_v\}$. Table 4 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 = \{x_i\}_{i=1}^n \in \mathcal{X} \subset \mathbb{R}^{n \times d_x}$ consists of a set of $n$ nodes,
each represented by a feature vector $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 \mathcal{Y}}{\operatorname{argmin}} C(y, x), \qquad (1)$$

where $\mathcal{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 $\mathcal{Y}$ consists of all possible node permutations. In CVRP, $x$ additionally includes demands, and $\mathcal{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.

## 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}(y_t|x, y_{1:t-1}),$$

with the convention $p_{\theta}(y_1|x, y_{1:0}) = p_{\theta}(y_1|x)$, where $y_t \in \{0, \dots, n\}$ is the selected node at step $t$, $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}(\cdot|x)} [C(y, x)].$$

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) dz.$$

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}(y_t|y_{1:t-1}, z, x),$$

where $y_t \in \{0, \dots, n\}$ is the selected node at step $t$, and $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><img src="image_1.png" /><figcaption>Figure 1: Our Latent Model Architecture with the Gaussian Reparameterization Trick</figcaption></figure>

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 \le t \le T$, for all $i \in \{0, \dots, n\}$, this layer computes the probability $p_{\theta}(y_t = i|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.

## 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|x)} \left[ \mathbb{E}_{y^k \sim p_{\theta}(\cdot|x, z^k)} [w^k C(y^k, x)] + \beta \mathcal{H}(p_{\theta}(\cdot|x, z^k)) \right]. \quad (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 = \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.
## 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)dz, \quad \text{where} \quad \pi_{\theta}(z, y|x) \propto p_{\phi}(z|x)p_{\theta}(y|z, x)e^{-\lambda C(y,x)}. \quad (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 $\exp(-\lambda C(y, x))$, where $\lambda$ controls the trade-off between likelihood and cost. However, incorporating this reweighting renders the distribution in (3) 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|x)} [C(y, x)] . $$

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}(x, \{(z^k, y^k)\}_{k=1}^K) = \frac{1}{K} \sum_{k=1}^{K} (C(y^k, x) - b(x)) \nabla_{\theta} \log p_{\theta}(y^k|x, z^k), \quad (4) $$

where $b(x)$ is a baseline function that reduces variance. The proposed method is detailed in Algorithm 1. 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, 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$.

# 5 Theoretical Results

Let $Z \subset \mathbb{R}^{d_z}$ be the latent space and $Y$ the solution space. In this section, we present theoretical results on our inference method described in Algorithm 1.

## 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). 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).
**Algorithm 1** Latent Guided Sampling

1: **Input:** Problem instance $x$, pretrained encoder $p_\phi$, pretrained decoder $p_{\theta_0}$, proposal distribution $q$, number of particles $K$, number of iterations $M$, cost function $C$, and temperature $\lambda$.

2: **Initialize:**

3: Sample initial particles: $z_0^k \sim p_\phi(\cdot|x)$ for all $k = 1, \dots, K$.

4: **for** $m = 0, 1, \dots, M-1$ **do**

5: Propagate new particles: $\tilde{z}_{m+1}^k \sim q(\cdot|z_m)$ for all $k = 1, \dots, K$.

6: Generate candidate solutions: $y_{m+1}^k \sim p_{\theta_m}(\cdot|\tilde{z}_{m+1}^k, x)$ for all $k = 1, \dots, K$.

7: Compute acceptance probabilities:

$$ \alpha_{m+1}^k = \min \left( 1, \frac{p_\phi(\tilde{z}_{m+1}^k | x)}{p_\phi(z_m^k | x)} \times e^{-\lambda(C(y_{m+1}^k, x) - C(y_m^k, x))} \right). $$

8: Accept $z_{m+1}^k = \tilde{z}_{m+1}^k$ with probability $\alpha_{m+1}^k$, otherwise $z_{m+1}^k = z_m^k$ for all $k = 1, \dots, K$.

9: Compute the gradient estimate $H_{\theta_m}(x, \{(z_{m+1}^k, y_{m+1}^k)\}_{k=1}^K)$ using (4).

10: Update parameters: $\theta_{m+1} = \theta_m - \gamma_{m+1} H_{\theta_m}(x, \{(z_{m+1}^k, y_{m+1}^k)\}_{k=1}^K)$.

11: **end for**

**Proposition 5.1** The sequence $\{(Z_m, Y_m) : m \in \mathbb{N}\}$ generated by Algorithm 1 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. 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 \le t \le T$, and all $(y_{1:t-1}, x, z)$, $p_\theta(y_t = i|y_{1:t-1}, x, z) > 0$ for all feasible nodes $i \in \{0, \dots, n\}$.

Since $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 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 and 2 hold. Then, the Markov Kernel $P_\theta$ admits a unique invariant probability measure $\pi_\theta$, defined in (3). There exist constants $\rho_1, \rho_2 \in (0, 1)$ and $\kappa_1, \kappa_2 \in \mathbb{R}_+$ such that for all $\mu \in M_1(\mathbb{Z} \times Y)$, $\theta \in \Theta$, and $m \in \mathbb{N}$,

$$ \|\mu P_\theta^m - \pi_\theta\|_{\text{TV}} \le \kappa_1 \rho_1^m \quad \text{and} \quad \|\mu P_\theta^m - \pi_\theta\|_{\text{L}_2} \le \kappa_2 \rho_2^m \|\mu - \pi_\theta\|_{\text{L}_2}. $$

Theorem 5.1 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 $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.
## 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. To analyze the convergence, we introduce the following additional assumptions.

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

$$ \|\nabla_{\theta} \log p_{\theta}(y|z, x)\| \le L(x, y, z). $$

Assumption 3 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 $(a_m)_{m \in \mathbb{N}}$, with $a_m \to 0$ as $m \to \infty$ such that

$$ \|\theta_m - \theta_{\infty}\|_{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).

**Theorem 5.2** Let Assumptions 1 - 4 hold. Then, there exist a constant $\rho \in (0, 1)$ and a positive sequence $(b_m)_{m \in \mathbb{N}}$ such that for all $\mu \in M_1(Z \times Y)$, and $m \in \mathbb{N}$,

$$ \mathbb{E}[\|\mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_{\infty}}\|_{TV}] = \mathcal{O}\left(\rho^{b_m} + \sum_{j=m-b_m}^{m-1} \gamma_{j+1} + a_m\right). $$

Furthermore, if $\limsup_{m \to \infty} (b_m^{-1} + b_m/m + b_m\gamma_m) = 0$, then:

$$ \mathbb{E}[\|\mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_{\infty}}\|_{TV}] \xrightarrow{m \to \infty} 0. $$

**Theorem 5.2** 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 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} (b_m^{-1} + b_m/m + b_m\gamma_m) = 0$, the expected total variation distance between the Markov Chain and the target distribution tends to zero as $m \to \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 $O(m^{-\gamma} \log m + a_m)$.

# 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<fnref n="1" />, with training and inference details in Appendix E. 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).

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

The average performance of each method is reported in Table 1 (TSP) and Table 2 (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.

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

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.

Table 3 and Figure 3 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, Paral- lel MCMC, Interacting MCMC, and LGS con- sistently outperform other techniques, particu- larly DE and CMA-ES, which are commonly used inference methods in continuous spaces. In contrast, Single MCMC struggles due to lim- ited exploration, leading to poor performance. EAS also underperforms, as the initial particles are insufficiently effective, and adjusting the pa- rameters does not significantly improve the so- lution. This highlights the critical importance of particle propagation and parameter learning in improving solution quality.

Additionally, the "Sampling" method corre-

sponds to direct sampling from the distribution defined in (3), without the reweighting factor $\exp(-\lambda C(y, x))$. The results highlight the advantage of incorporating this reweighting factor, which

<figure><img src="image_2.png" /><figcaption>Figure 2: Visualization of the 2-dimensional la- tent space ($d_z = 2$) learned by our model on a problem instance. The plotted path illustrates the search trajectory leading to the best-found solu- tion.</figcaption></figure>
<table><caption>Table 3: Comparison of different inference methods with our model on CVRP with $n = 100$</caption><thead><tr><td>Method</td><td>Obj.</td><td>Gap</td></tr></thead><tbody><tr><td>Sampling</td><td>15.652</td><td>0.721%</td></tr><tr><td>DE</td><td>15.561</td><td>0.135%</td></tr><tr><td>CMA-ES</td><td>15.582</td><td>0.271%</td></tr><tr><td>EAS</td><td>15.685</td><td>0.933%</td></tr><tr><td>Single MCMC</td><td>15.649</td><td>0.701%</td></tr><tr><td>Parallel MCMC (ours)</td><td>15.557</td><td>0.109%</td></tr><tr><td>Interacting MCMC (ours)</td><td>15.535</td><td>-0.032%</td></tr><tr><td>LGS (ours)</td><td><strong>15.524</strong></td><td><strong>-0.102%</strong></td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 3: Performance of sampling-based methods on CVRP with $n = 100$, with bold lines indicating the mean over 10 inference runs.</figcaption></figure>

improves the exploration of the latent space, as shown in Figure 2. While interacting MCMC benefits from faster mixing (Theorem 5.1), it is nevertheless outperformed by LGS, which achieves better cost convergence despite potentially slower mixing (Theorem 5.2). 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

C. Andrieu and Y. F. Atchadé. On the efficiency of adaptive mcmc algorithms. *Electronic Communications in Probability*, 12:336-349, 2007.

C. Andrieu and É. Moulines. On the ergodicity properties of some adaptive mcmc algorithms. *Annals of Applied Probability*, 16(3):1462-1505, 2006.

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.

D. Applegate, R. Bixby, V. Chvatal, and W. Cook. Concorde tsp solver, 2006.

P. H. Baxendale. Renewal theory and computable convergence rates for geometrically ergodic markov chains. *Annals of Applied Probability*, 15:700-738, 2005.

I. Bello, H. Pham, Q. V. Le, M. Norouzi, and S. Bengio. Neural combinatorial optimization with reinforcement learning. *arXiv preprint arXiv:1611.09940*, 2016.

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.

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.

X. Chen and Y. Tian. Learning to perform local rewriting for combinatorial optimization. In *Advances in Neural Information Processing Systems*, volume 32, 2019.

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.

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.

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.

J. Devlin. Bert: Pre-training of deep bidirectional transformers for language understanding. *arXiv preprint arXiv:1810.04805*, 2018.

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.

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.
R. Douc, E. Moulines, P. Priouret, P. Soulier, R. Douc, E. Moulines, P. Priouret, and P. Soulier. *Markov chains: Basic definitions*. Springer, 2018.

S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth. Hybrid monte carlo. *Physics letters B*, 195(2):216-222, 1987.

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.

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.

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.

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.

F. Glover. Tabu search—part i. *ORSA Journal on computing*, 1(3):190-206, 1989.

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.

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.

H. Haario, E. Saksman, and J. Tamminen. An adaptive metropolis algorithm. *Bernoulli*, 7(2): 223-242, 2001.

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.

N. Hansen and A. Ostermeier. Completely derandomized self-adaptation in evolution strategies. *Evolutionary computation*, 9(2):159-195, 2001.

W. K. Hastings. Monte carlo sampling methods using markov chains and their applications. *Biometrika*, 1970.

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.

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.

J. J. Hopfield and D. W. Tank. "neural" computation of decisions in optimization problems. *Biological cybernetics*, 52(3):141-152, 1985.
A. Hottung and K. Tierney. Neural large neighborhood search for the capacitated vehicle routing problem. In *ECAI 2020*, pages 443–450. IOS Press, 2020.

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.

A. Hottung, Y.-D. Kwon, and K. Tierney. Efficient active search for combinatorial optimization problems. In *International Conference on Learning Representations*, 2022.

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.

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.

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.

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.

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.

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.

D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In *International Conference on Learning Representations*, 2015.

D. P. Kingma and M. Welling. Auto-encoding variational bayes. In *International Conference on Learning Representations*, 2014.

S. Kirkpatrick, C. D. Gelatt Jr, and M. P. Vecchi. Optimization by simulated annealing. *science*, 220(4598):671–680, 1983.

V. Konda and J. Tsitsiklis. Actor-critic algorithms. In *Advances in Neural Information Processing Systems*, volume 12, 1999.

W. Kool, H. Van Hoof, and M. Welling. Attention, learn to solve routing problems! In *International Conference on Learning Representations*, 2019.

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.

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.

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

R. Liu, X. Li, and K. S. Lam. Combinatorial chemistry in drug discovery. *Current opinion in chemical biology*, 38:117–126, 2017.

H. Lu, X. Zhang, and S. Yang. A learning-based iterative method for solving vehicle routing problems. In *International conference on learning representations*, 2019.

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.

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.

S. P. Meyn and R. L. Tweedie. *Markov chains and stochastic stability*. Springer Science & Business Media, 2012.

N. Mladenović and P. Hansen. Variable neighborhood search. *Computers & operations research*, 24 (11):1097–1100, 1997.

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.

C. H. Papadimitriou and K. Steiglitz. *Combinatorial optimization: algorithms and complexity*. Courier Corporation, 1998.

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.

L. Perron and V. Furnon. OR-Tools. https://developers.google.com/optimization/, 2019.

H. Robbins and S. Monro. A stochastic approximation method. *The Annals of Mathematical Statistics*, pages 400–407, 1951.

G. Roberts and J. Rosenthal. Geometric ergodicity and hybrid markov chains. *Electronic Communications in Probability*, 2(2):13–25, 1997.

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.

V. Steinbiss, B.-H. Tran, and H. Ney. Improvements in beam search. In *ICSLP*, volume 94, pages 2143–2146, 1994.

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.

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

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.

A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin. Attention is all you need. In *Advances in Neural Information Processing Systems*, volume 30, 2017.

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.

O. Vinyals, M. Fortunato, and N. Jaitly. Pointer networks. In *Advances in Neural Information Processing Systems*, volume 28, 2015.

R. J. Williams. Simple statistical gradient-following algorithms for connectionist reinforcement learning. *Machine learning*, 8:229-256, 1992.

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.

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.

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.
# Supplementary Material

## Table of Contents

<table><tr><td>A</td><td>Problem and Model Description</td><td>19</td></tr><tr><td></td><td>A.1 Problem Setting</td><td>19</td></tr><tr><td></td><td>A.2 Model Architecture details</td><td>20</td></tr><tr><td></td><td>A.3 Training</td><td>21</td></tr><tr><td>B</td><td>Preliminaries on Markov Chains</td><td>23</td></tr><tr><td>C</td><td>Convergence Analysis for Fixed θ</td><td>25</td></tr><tr><td></td><td>C.1 Proof of Proposition 5.1</td><td>25</td></tr><tr><td></td><td>C.2 Proof of Theorem 5.1 for K = 1</td><td>27</td></tr><tr><td></td><td>C.3 Extension to K > 1</td><td>29</td></tr><tr><td>D</td><td>Convergence Analysis for Adaptive θ</td><td>31</td></tr><tr><td></td><td>D.1 Convergence Analysis of Time-Inhomogeneous MCMC algorithm with Stochastic Approximation Update</td><td>31</td></tr><tr><td></td><td>D.2 Proof of Theorem 5.2</td><td>36</td></tr><tr><td></td><td>D.3 Convergence Rate in Stochastic Approximation</td><td>36</td></tr><tr><td>E</td><td>Additional Experiments</td><td>37</td></tr><tr><td></td><td>E.1 Training Details</td><td>37</td></tr><tr><td></td><td>E.2 Inference Details and Additional Illustrations</td><td>37</td></tr></table>
# Notation

<table><caption>Table 4: Summary of notation used throughout the paper.</caption><thead><tr><th>Object</th><th>Description</th></tr></thead><tbody><tr><td>x = {x<sub>i</sub>}<sup>n</sup><sub>i=1</sub> ∈ X ⊂ ℝ<sup>n×d<sub>x</sub></sup></td><td>Problem instance</td></tr><tr><td>y = (y<sub>1</sub>, ..., y<sub>T</sub>) ∈ Y ⊂ {0, ..., n}<sup>T</sup></td><td>Solution</td></tr><tr><td>z ∈ Z ⊂ ℝ<sup>d<sub>z</sub></sup></td><td>Latent variable</td></tr><tr><td>(X, X)</td><td>Problem instance space X with Borel σ-algebra X = 𝒩(X)</td></tr><tr><td>(Y, Y)</td><td>Discrete solution space Y with the power set Y = 𝒫(Y)</td></tr><tr><td>(Z, Z)</td><td>Latent space Z with Borel σ-algebra Z = 𝒩(Z)</td></tr><tr><td>ℝ<sub>x</sub></td><td>Distribution over problem instances</td></tr><tr><td>p<sub>φ</sub>(z|x)</td><td>Encoder distribution</td></tr><tr><td>p<sub>θ</sub>(y|x, z)</td><td>Decoder distribution</td></tr><tr><td>C</td><td>Cost function</td></tr><tr><td>n</td><td>Number of nodes in the instance</td></tr><tr><td>t, T</td><td>Decoding step index and horizon</td></tr><tr><td>m, M</td><td>Inference step index and total iterations</td></tr><tr><td>k, K</td><td>Particle index and total number of particles</td></tr><tr><td>B</td><td>Batch size used during training</td></tr></tbody></table>

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.

# A Problem and Model Description

## A.1 Problem Setting

**Traveling Salesman Problem (TSP).** A TSP instance $x = \{x_i\}_{i=1}^n$ consists of a set of $n$ nodes, where the feature $x_i$ corresponds to its coordinates $c_i \in ℝ^2$. The objective is to find a permutation $y = (y_1, ..., y_n)$ of the nodes, where $y_t \in \{1, ..., n\}$ and $y_t \neq y_{t'}$ for all $t \neq t'$, that minimizes the total tour length:

$$C(y, x) = \sum_{i=1}^{n-1} \|x_{y_{i+1}} - x_{y_i}\| + \|x_{y_n} - x_{y_1}\|, \quad (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, ..., n\}$ has a demand $d_i > 0$ and a location $c_i \in ℝ^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 $x_i$, the encoder initially computes $d_h$-dimensional node embeddings $h_i^{(0)}$ through a learned linear projection using parameters $W_0$ and $b_0$:

$$ h_i^{(0)} = W_0 x_i + b_0. $$

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

$$ \hat{h}_{i}^{(l+1)} = \text{InstanceNorm} \left( h_{i}^{(l)} + \text{MHA} \left( h_{1}^{(l)}, \dots, h_{n}^{(l)} \right) \right), $$

$$ h_i^{(l+1)} = \text{InstanceNorm} \left( \hat{h}_i^{(l+1)} + \text{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_i^{(L)}$. 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) = \text{FF}(\bar{h}^{(L)}) \quad \text{and} \quad \log \sigma_{\phi}(x)^2 = \text{FF}(\bar{h}^{(L)}). $$

## 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 = [z, h_{y_{t-1}}^{(L)}, h_{y_0}^{(L)}], $$

where $h_{y_0}^{(L)}$ and $h_{y_{t-1}}^{(L)}$ 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}(y_t = i | y_{0:t-1}, x, z) = \begin{cases} \mathrm{softmax} \left( \omega \tanh \left( \frac{q_{(c_t)}^T k_i}{\sqrt{d_k}} \right) \right) & \text{if } i \neq y_s \quad \forall s < t, \\ 0 & \text{otherwise,} \end{cases}$$

where the query and key are given by $q_{(c_t)} = \mathrm{MHA}(c_t, \{h_i^{(L)}\}_{i=1}^n, \{h_i^{(L)}\}_{i=1}^n)$ and $k_i = W^K h_i^{(L)}$ 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 = [z, h_{y_{t-1}}^{(L)}, \hat{D}_t],$$

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

$$\hat{D}_{t+1} = \begin{cases} \max(\hat{D}_t - d_{y_t, t}, 0) & \text{if } y_t \neq 0, \\ D & \text{if } 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}(y_t = i | y_{0:t-1}, x, z) = \begin{cases} \mathrm{softmax} \left( \omega \tanh \left( \frac{q_{(c_t)}^T k_i}{\sqrt{d_k}} \right) \right) & \text{if } i \neq y_s \quad \forall s < t \text{ and } d_{i,t} \le \hat{D}_t, \\ 0 & \text{otherwise,} \end{cases}$$

where the query and key are given by $q_{(c_t)} = \mathrm{MHA}(c_t, \{h_i^{(L)}\}_{i=1}^n, \{h_i^{(L)}\}_{i=1}^n)$ and $k_i = W^K h_i^{(L)}$ respectively.

## A.3 Training

The estimator of the gradient of the objective defined in (2) is computed using the Monte Carlo method:

$$\begin{aligned} \hat{\nabla}_{\theta} \mathcal{L}(\theta, \phi) = & \frac{1}{B} \sum_{i=1}^{B} \sum_{k=1}^{K} w_{(i)}^{k} (C(y_{(i)}^{k}, x_{(i)}) - b(x_{(i)})) \nabla_{\theta} \log p_{\theta}(y_{(i)}^{k} | x_{(i)}, z_{(i)}^{k}) \\ & - \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}), \end{aligned} \quad (6)$$
$$ \hat{\nabla}_{\phi} \mathcal{L}(\theta, \phi) = \frac{1}{B} \sum_{i=1}^{B} \sum_{k=1}^{K} w_{(i)}^{k} (C(y_{(i)}^{k}, x_{(i)}) - b(x_{(i)})) \nabla_{\phi} \log p_{\phi}(z_{(i)}^{k} | x_{(i)}), \quad (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. The update step ADAM for the parameters $(\theta, \phi)$ corresponds to a single Adam update step (Kingma and Ba, 2015).

Algorithm 2 REINFORCE training

**Input:** Distribution over problem instances $\mathbb{P}_x$, number of training steps $N$, batch size $B$, and number of latent samples $K$.

1. Initialize the model parameters $\theta$ and $\phi$.

2. **for** epoch = 1 to $N$ **do**

3. Sample problem instances $x_{(i)} \sim \mathbb{P}_x$ for $i \in \{1, \dots, B\}$.

4. Generate latent samples $z_{(i)}^1, \dots, z_{(i)}^K \sim p_{\phi}^{\otimes K}(\cdot | x_{(i)})$ using the reparameterization trick.

5. Sample solutions $y_{(i)}^k \sim p_{\theta}(\cdot | x_{(i)}, z_{(i)}^k)$ for all $k = 1, \dots, K$.

6. Compute the gradient estimates $\hat{\nabla}_{\theta,\phi} \mathcal{L}(\theta, \phi)$ using (6) and (7).

7. Update parameters: $(\theta, \phi) \leftarrow \text{ADAM}((\theta, \phi), \hat{\nabla}_{\theta,\phi} \mathcal{L}(\theta, \phi))$.

8. **end for**

**Output:** Optimized parameters $\theta$ and $\phi$.
# B Preliminaries on Markov Chains

In this section, we use the following definitions.

* A sequence of random variables {$X_n, n \in \mathbb{N}$} is a Markov Chain with respect to the filtration ($\mathcal{F}_n$)$_{n \ge 0}$ with Markov kernel $P: X \times \mathcal{X} \to \mathbb{R}_+$ if for any bounded measurable function $f: X \to \mathbb{R}$,

$$ \mathbb{E}[f(X_{n+1}) | \mathcal{F}_n] = P f(X_n) = \int f(x) P(X_n, dx) . $$

* Furthermore, the sequence {$X_n, n \in \mathbb{N}$} is a state-dependent Markov Chain if for any bounded measurable function $f: X \to \mathbb{R}$,

$$ \mathbb{E}[f(X_{n+1}) | \mathcal{F}_n] = P_{\theta_n} f(X_n) = \int f(x) P_{\theta_n}(X_n, dx) , $$

where $P_{\theta_n}: X \times \mathcal{X} \to \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 $(X, \mathcal{X})$ is called invariant for the Markov kernel $P$ if it satisfies $\pi P = \pi$. If {$X_n: n \in \mathbb{N}$} is a Markov Chain with Markov kernel $P$ and $X_0$ is distributed according to an invariant probability measure $\pi$, then for all $n \ge 1$, we have $X_n \sim \pi$.

**Definition B.2 (Reversibility)**

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

$$ \pi(dx)P(x, dy) = \pi(dy)P(y, dx). $$

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

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

$$ \mathcal{C}(\mu, \nu) = \{ \zeta \in \mathbb{M}_1(\mathcal{X}^2) : \forall A \in \mathcal{X}, \zeta(A \times X) = \mu(A), \zeta(X \times A) = \nu(A) \} . $$

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

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

$$ \begin{aligned} \|\mu - \nu\|_{\text{TV}} &= 2 \sup \{|\mu(f) - \nu(f)| : f \in \mathcal{F}(X), 0 \le f \le 1\} \\ &= 2 \inf \{\zeta(\Delta) : \zeta \in \mathcal{C}(\mu, \nu)\}, \end{aligned} $$

where $\Delta(x, x') = \mathbf{1}_{x \ne x'}$ for all $(x, x') \in \mathcal{X}^2$.

**Assumption 5** Let $P$ be a Markov transition kernel on $(X, \mathcal{X})$. Suppose there exists a function $V: X \to [0, \infty)$ satisfying $\sup_{x \in 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) \ge \varepsilon\nu(A) . $$
2. Drift Condition (Foster-Lyapunov Condition). There exist constants $\lambda \in [0, 1)$, $b \in (0, \infty)$ satisfying

$$PV(x) \le \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 hold for a function $V: X \to [0, \infty)$, where $\sup_{x \in 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 M_1(X)$, and $m \in \mathbb{N}$,

$$\|\mu P^m - \pi\|_{\text{TV}} \le \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).
C Convergence Analysis for Fixed θ

In this section, we prove that the iterates of Algorithm 1 with K = 1 and fixed θ (the exact algorithm is given in Algorithm 3) form a reversible Markov Chain (Proposition 5.1) and exhibit geometric ergodicity toward the joint target distribution (Theorem 5.1). We then extend these results to the general case of arbitrary K in Section C.3.

**Algorithm 3** Latent Guided Sampling (K = 1)

1: **Input:** Problem instance $x$, pretrained encoder $p_{\phi}$, pretrained decoder $p_{\theta}$, proposal distribution $q$, number of iterations $M$, and cost function $C$.

2: **Initialize:**

3: Sample initial particle: $z_0 \sim p_{\phi}(\cdot|x)$.

4: **for** $m = 0, 1, \dots, M-1$ **do**

5: Propagate new particle: $\tilde{z}_{m+1} \sim q(\cdot|z_m)$.

6: Generate new solution: $y_{m+1} \sim p_{\theta}(\cdot|\tilde{z}_{m+1}, x)$.

7: Compute the acceptation probability:

$$ \alpha_{m+1} = \min \left( 1, e^{-\lambda(C(y_{m+1},x)-C(y_m,x))} \frac{p_{\phi}(\tilde{z}_{m+1}|x)}{p_{\phi}(z_m|x)} \right). $$

8: Accept $z_{m+1} = \tilde{z}_{m+1}$ with probability $\alpha_{m+1}$.

9: **end for**

C.1 Proof of Proposition 5.1

**Proposition C.1** The sequence {($Z_m, Y_m$) : $m \in \mathbb{N}$} generated by Algorithm 3 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 \mathbb{Z}, y \in \mathbb{Y}$, and $A \in \mathcal{Z} \times \mathcal{Y}$,

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

where

$$ \alpha(z, y, z', y') = \min \left( 1, e^{-\lambda(C(y',x)-C(y,x))} \frac{p_{\phi}(z'|x)}{p_{\phi}(z|x)} \right), $$

$$ \bar{\alpha}_{\theta}(z, y) = 1 - \int_{\mathbb{Z} \times \mathbb{Y}} q(dz'|z) p_{\theta}(dy'|z', x) \alpha(z, y, z', y'). $$

*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} = (U_1, \dots, U_m)$ denotes the sequence of uniform random variables. For all bounded or
non-negative measurable function $h$ on $Z \times Y$ and all $m \in \mathbb{N}$,

$$
\begin{align*}
\mathbb{E}[h(Z_{m+1}, Y_{m+1}) | \mathcal{F}_m] \\
&= \mathbb{E} \left[ 1_{\{U_{m+1} < \alpha(Z_m, Y_m, Z'_{m+1}, Y'_{m+1})\}} h(Z'_{m+1}, Y'_{m+1}) \mid \mathcal{F}_m \right] \\
&\quad + \mathbb{E} \left[ 1_{\{U_{m+1} \ge \alpha(Z_m, Y_m, Z'_{m+1}, Y'_{m+1})\}} h(Z_m, Y_m) \mid \mathcal{F}_m \right] \\
&= \int_{Z \times Y} q(dz'|Z_m) p_\theta(dy'|z', x) \alpha(Z_m, Y_m, z', y') h(z, y) + \bar{\alpha}_\theta(Z_m, Y_m) h(Z_m, Y_m),
\end{align*}
$$

where

$$
\bar{\alpha}_{\theta}(Z_{m}, Y_{m}) = 1 - \int_{Z \times Y} q(dz'|Z_{m}) p_{\theta}(dy'|z', x) \alpha(Z_{m}, Y_{m}, z', y').
$$

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

$$
P_{\theta}((z, y), A) = \int_{A} q(dz'|z) p_{\theta}(dy'|z', x) \alpha(z, y, z', y') + \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}((z, y), (z', y')) = \pi_{\theta}(z', y') P_{\theta}((z', y'), (z, y)). \quad (9)
$$

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

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

We separate the analysis in two cases depending on the value of $r(z, y, z', y')$.

**Case 1.** If $r(z, y, z', y') \le 1$, then $\alpha(z, y, z', y') = r(z, y, z', y')$ and $\alpha(z', y', z, y) = 1$. Thus,

$$
\begin{align*}
\pi_{\theta}(z, y) P_{\theta}((z, y), (z', y')) &= p_{\phi}(z|x) p_{\theta}(y|z, x) e^{-\lambda C(y,x)} P_{\theta}((z, y), (z', y')) \\
&= p_{\phi}(z|x) p_{\theta}(y|z, x) e^{-\lambda C(y,x)} \\
&\quad \times q(z'|z) p_{\theta}(y'|z', x) e^{-\lambda(C(y',x)-C(y,x))} \frac{p_{\phi}(z'|x)}{p_{\phi}(z|x)} \\
&= p_{\phi}(z'|x) p_{\theta}(y'|z', x) e^{-\lambda C(y',x)} q(z'|z) p_{\theta}(y|z, x) \\
&= \pi_{\theta}(z', y') P_{\theta}((z', y'), (z, y)).
\end{align*}
$$

**Case 2.**
If $r(z, y, z', y') > 1$, then $\alpha(z, y, z', y') = 1$ and $\alpha(z', y', z, y) = r(z', y', z, y)$. Similarly,

$$
\begin{align*}
\pi_{\theta}(z', y') P_{\theta}((z', y'), (z, y)) &= p_{\phi}(z'|x) p_{\theta}(y'|z', x) e^{-\lambda C(y', x)} P_{\theta}((z', y'), (z, y)) \\
&= p_{\phi}(z'|x) p_{\theta}(y'|z', x) e^{-\lambda C(y', x)} \\
&\quad \times q(z|z') p_{\theta}(y|z, x) e^{-\lambda(C(y,x)-C(y',x))} \frac{p_{\phi}(z|x)}{p_{\phi}(z'|x)} \\
&= p_{\phi}(z|x) p_{\theta}(y|z, x) e^{-\lambda C(y,x)} q(z|z') p_{\theta}(y'|z', x) \\
&= \pi_{\theta}(z, y) P_{\theta}((z, y), (z', y')).
\end{align*}
$$

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

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

*Proof.* This follows from Theorem B.1, provided that we verify the minorization and drift conditions of Assumption 5. We consider the set

$$ \mathcal{K} = \left\{ (z, y) \in \mathbb{Z} \times \mathcal{Y} \mid \|z\|^2 \le R^2, C(y) \ne \max_{y' \in \mathcal{Y}} C(y') \right\} $$

which is compact since $\mathcal{Y}$ is finite. We denote by $B(c, R) = \{z \in \mathbb{Z} \mid \|z - c\|^2 \le 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}$,

$$
\begin{align*}
P_{\theta}((z, y), A) \\
&= \int_{A} q(dz'|z) p_{\theta}(dy'|z', x) \min \left( 1, e^{-\lambda(C(y', x)-C(y,x))} \frac{p_{\phi}(z'|x)}{p_{\phi}(z|x)} \right) + \bar{\alpha}_{\theta}(z, y) \delta_{(z,y)}(A) \\
&\geq \int_{A \cap \mathcal{K}} q(dz'|z) p_{\theta}(dy'|z', x) \min \left( 1, e^{-\lambda(C(y', x)-C(y,x))} \frac{p_{\phi}(z'|x)}{p_{\phi}(z|x)} \right) + \bar{\alpha}_{\theta}(z, y) \delta_{(z,y)}(A) .
\end{align*}
$$

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' \in B(0, R)$,

$$ q(z' | z) \ge \varepsilon_q := \inf_{z, z' \in B(0, R)} q(z' | z) > 0. $$

* Encoder Component: By Assumption 1, $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'|x)}{p_{\phi}(z|x)} \ge \varepsilon_e . $$

* Decoder Component: By Assumption 1, the categorical transition probability satisfies:

$$ \inf_{z' \in \mathcal{Z}, y' \in \mathcal{Y}} p_{\theta}(y'|z', x) \ge \varepsilon_d > 0. $$
• Exponential Weighting: Since $e^{-\lambda(C(y',x)-C(y,x))}$ is always positive and C is bounded (Assumption 1), there exists $\varepsilon_w > 0$ such that:

$$e^{-\lambda(C(y',x)-C(y,x))} \geq \varepsilon_w. \qquad (10)$$

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

$$ \begin{aligned} P_{\theta}((z, y), A) &= \int_{A} q(dz'|z) p_{\theta}(dy'|z', x) \min \left( 1, e^{-\lambda(C(y',x)-C(y,x))} \frac{p_{\phi}(z'|x)}{p_{\phi}(z|x)} \right) + \bar{\alpha}_{\theta}(z, y) \delta_{(z,y)}(A) \\ &\geq \mu^{\text{Leb}}(B(0, R)) |\mathcal{Y}| \varepsilon_q \varepsilon_d \min(1, \varepsilon_w \varepsilon_e) \int_{A} \nu_q(dz') \nu_d(dy'), \end{aligned} $$

where $\nu_C$ is the uniform probability measure over $\mathcal{Y}$:

$$\nu_d(dy') = \frac{1}{|\mathcal{Y}|} \sum_{y \in \mathcal{Y}} \delta_y(dy')$$

and

$$\nu_q(dz') := \frac{\mu^{\text{Leb}}(dz')}{\mu^{\text{Leb}}(B(0, R))} 1_{B(0, R)}(z'),$$

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

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

where $\varepsilon = \mu^{\text{Leb}}(B(0, R))|\mathcal{Y}|\varepsilon_q\varepsilon_d \min(1, \varepsilon_w\varepsilon_e)$ and $\nu(dz', dy') = \nu_q(dz')\nu_d(dy')$ is a probability measure.

### Drift condition.

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

$$V_x(z, y) = e^{sC(y,x)}.$$

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

$$ \begin{aligned} P_{\theta}V_x(z, y) &= \int e^{sC(y', x)} (q(dz'|z) p_{\theta}(dy'|z', x) \alpha(z, y, z', y') + \bar{\alpha}_{\theta}(z, y) \delta_{(z,y)}(dz', dy')) \\ &= \int e^{sC(y', x)} q(dz'|z) p_{\theta}(dy'|z', x) \alpha(z, y, z', y') + \bar{\alpha}_{\theta}(z, y) e^{sC(y,x)} \\ &= \int \left( e^{sC(y', x)} \alpha(z, y, z', y') + e^{sC(y,x)} - \alpha(z, y, z', y') e^{sC(y,x)} \right) q(dz'|z) p_{\theta}(dy'|z', x). \end{aligned} $$

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 $(z, y) \notin \mathcal{K}$,

$$ 
\begin{aligned}
& \frac{P_{\theta}V_x(z, y)}{V_x(z, y)} \\
&= \int \left( e^{s(C(y',x)-C(y,x))} \alpha(z, y, z', y') + 1 - \alpha(z, y, z', y') \right) q(dz'|z) p_{\theta}(dy'|z', x) \\
&\le \int_{y' \in \{C(y',x)<C(y,x)\}} \left( e^{s(C(y',x)-C(y,x))} \alpha(z, y, z', y') + 1 - \alpha(z, y, z', y') \right) q(dz'|z) p_{\theta}(dy'|z', x) \\
&\quad + \int_{y' \in \{C(y',x)=C(y,x)\}} \left( e^{s(C(y',x)-C(y,x))} \alpha(z, y, z', y') + 1 - \alpha(z, y, z', y') \right) q(dz'|z) p_{\theta}(dy'|z', x).
\end{aligned}
$$

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

$$ 
\begin{aligned}
& \int_{y' \in \{C(y',x)<C(y,x)\}} \left( e^{s(C(y',x)-C(y,x))} \alpha(z, y, z', y') + 1 - \alpha(z, y, z', y') \right) q(dz'|z) p_{\theta}(dy'|z', x) \\
&\le \eta \int_{y' \in \{C(y',x)<C(y,x)\}} q(dz'|z) p_{\theta}(dy'|z', x).
\end{aligned}
$$

Thus,

$$ I(z, y) \le \eta + (1 - \eta) \int_{y' \in \{C(y',x)=C(y,x)\}} q(dz'|z) p_{\theta}(dy'|z', x), $$

which establishes the drift condition and completes the proof of geometric ergodicity in total variation distance using Theorem B.1.

For $L_2$-geometric ergodicity, note that the Markov Chain is reversible with respect to the probability measure $\pi_{\theta}$ (Proposition 5.1), 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 M_1(\mathbb{Z} \times \mathcal{Y})$, $\theta \in \Theta$, and $m \in \mathbb{N}$,

$$ \|\mu P_{\theta}^m - \pi_{\theta}\|_{L_2} \le \kappa_2 \rho_2^m \|\mu - \pi_{\theta}\|_{L_2}. $$

□

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

Here, we show how Theorem 5.1 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 \ge 1$ is defined, for all $A \in \mathcal{Z}^K \times \mathcal{Y}^K$, as

$$ P_{\theta}((z^{1:K}, y^{1:K}), A) = \int_{A} \prod_{k=1}^{K} q(d\tilde{z}^k|z) p_{\theta}(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)}(d\tilde{z}^k, d\tilde{y}^k), $$
where $\alpha$ and $\bar{\alpha}_{\theta}$ are defined in (8).

Using the same arguments as in the case $K=1$ (cf. Section C.2), we obtain the following minoriza-
tion condition:

$$
P_{\theta}((z^{1:K}, y^{1:K}), A) \geq \varepsilon^{K} \nu^{\otimes K}(A),
$$

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

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

As in the case $K=1$, for all $(z^{1:K}, y^{1:K}) \in \mathcal{K}^K$, applying $P_\theta$ to $V_x^K$ yields:

$$
\begin{align*}
& P_{\theta} V_x^K (z^{1:K}, y^{1:K}) \\
&= \sum_{k=1}^{K} \int e^{sC(\tilde{y}^k, x)} (q(d\tilde{z}^k|z) p_{\theta}(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)}(d\tilde{z}^k, d\tilde{y}^k)) \\
&= \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(d\tilde{z}^k|z) p_{\theta}(d\tilde{y}^k|\tilde{z}^k, x) .
\end{align*}
$$

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 \ge 1$.

<figure><img src="image_4.png" /></figure>
# 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 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.

<caption>**Algorithm 4** Time-Inhomogeneous MCMC with Stochastic Approximation Update</caption>

1: **Input:** Number of iterations $M$, initial parameter $\theta_0$, step sizes $(\gamma_m)_{m \ge 1}$, initial distribution $\mu$, proposal distribution $q$, and target distribution $\pi_{\theta}$.
2: **Initialize:**
3: Sample initial particle: $x_0 \sim \mu$.
4: **for** $m = 0, 1, \dots, M-1$ **do**
5: Propose a new sample: $\tilde{x}_{m+1} \sim q(\cdot|x_m)$
6: Compute acceptance probability:
$$ \alpha_{m+1} = \min \left( 1, \frac{\pi_{\theta_m}(\tilde{x}_{m+1})q(x_m|\tilde{x}_{m+1})}{\pi_{\theta_m}(x_m)q(\tilde{x}_{m+1}|x_m)} \right). $$
7: Accept $x_{m+1} = \tilde{x}_{m+1}$ with probability $\alpha_{m+1}$.
8: Compute the gradient estimate $H_{\theta_m}(x_{m+1})$ using previous samples.
9: Update parameters: $\theta_{m+1} = \theta_m - \gamma_{m+1}H_{\theta_m}(x_{m+1})$.
10: **end for**

To analyze its convergence, we introduce the following assumptions.

**Assumption 6** Let $(P_k)_{k \ge 1}$ be a sequence of Markov transition kernels on $(X, \mathcal{X})$. Suppose there exists a function $V: X \to [1, \infty)$ satisfying $\sup_{x \in 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) \ge \varepsilon\nu(A). $$

2. **Drift condition.** There exist constants $\lambda \in [0, 1)$, $b \in (0, \infty)$ satisfying

$$ P_k V(x) \le \begin{cases} \lambda V(x) & x \notin \mathcal{K}, \\ b & x \in \mathcal{K}. \end{cases} $$

Assumption 6 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 F(X)$ such that for all $x \in X$ and $\theta \in \Theta$,

$$ \|\nabla_{\theta} \log \pi_{\theta}(x)\| \le L(x). $$

**Assumption 8** There exists $\theta_{\infty} \in \Theta$ and a positive sequence $(a_m)_{m \in \mathbb{N}}$, with $a_m \to 0$ as $m \to \infty$ such that

$$ \|\theta_m - \theta_{\infty}\|_{L_2}^2 = O(a_m). $$

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

**Theorem D.1** Let Assumptions 6 - 8 hold. Then, there exist a constant $\rho \in (0, 1)$ and positive sequences $(b_m)_{m \in \mathbb{N}}$ such that for all $\mu \in M_1(X)$, and $m \in \mathbb{N}$,

$$ \mathbb{E} \left[ \| \mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_{\infty}} \|_{\text{TV}} \right] = \mathcal{O} \left( \rho^{b_m} + \sum_{j=m-b_m}^{m-1} \gamma_{j+1} + a_m \right). $$

Furthermore, if $\limsup_{m \to \infty} (b_m^{-1} + b_m/m + b_m\gamma_m) = 0$, then

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

Theorem D.1 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, followed by an upper bound on the mixing error in D.1.2. The proof of the theorem is then provided in D.1.3.

### D.1.1 Error Decomposition

**Lemma D.1** For all $1 \le s \le m$, we have:

$$ \|\mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_{\infty}}\|_{\text{TV}} \le \|\mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_s} P_{\theta_{s+1}} \cdots P_{\theta_m}\|_{\text{TV}} + \sum_{j=s}^{m-1} \|\pi_{\theta_{j+1}} - \pi_{\theta_j}\|_{\text{TV}} \\ + \|\pi_{\theta_m} - \pi_{\theta_{\infty}}\|_{\text{TV}}. $$

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

$$ \|\mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_{\infty}}\|_{\text{TV}} \le \|\mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_m}\|_{\text{TV}} + \|\pi_{\theta_m} - \pi_{\theta_{\infty}}\|_{\text{TV}}. \quad (12) $$

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

$$ \begin{aligned} \|\mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_m}\|_{\text{TV}} &\le \|\mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_s} P_{\theta_{s+1}} \cdots P_{\theta_m}\|_{\text{TV}} \\ &\quad + \sum_{j=s}^{m-1} \|\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}\|_{\text{TV}} \\ &\le \|\mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_s} P_{\theta_{s+1}} \cdots P_{\theta_m}\|_{\text{TV}} + \sum_{j=s}^{m-1} \|\pi_{\theta_{j+1}} - \pi_{\theta_j}\|_{\text{TV}}, \end{aligned} $$

where the last inequality follows from the fact that, for all $j$, the Markov kernels $P_{\theta_j}$ are contractions. Together with (12), 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}$.

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

**Proposition D.1** Let Assumption 6 hold for a function $V: X \to [1, \infty)$, where $\sup_{x \in X} V(x) < \infty$. Let $(P_k)_{k \ge 1}$ be a sequence of Markov transition kernels on $(X, \mathcal{X})$. Then, there exists a constant $\rho \in (0, 1)$ such that for all $\xi, \xi' \in M_1(X)$ and all $m \in \mathbb{N}$,

$$ \|\xi P_1 \cdots P_m - \xi' P_1 \cdots P_m\|_{TV} \le \begin{cases} 2\rho^m & \text{if } \xi, \xi' \text{ are supported on } \mathcal{K}, \\ \rho^m (\xi(V) + \xi'(V)) & \text{otherwise.} \end{cases} $$

Proposition D.1 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'_k)_{k \ge 1}$ such that, marginally, $(X_k)_{k \ge 1}$ and $(X'_k)_{k \ge 1}$ are Markov Chains starting from $X_1 = x$ and $X'_1 = x'$, respectively, and each evolving according to the transition kernels $(P_k)_{k \ge 1}$.

To achieve this, for all $k \ge 1$, we define the modified kernel $Q_k$, defined for all $A \subset \mathcal{X}$ and $x_k \in X$ as:

$$ Q_k(x_k, A) = \frac{P_k(x_k, A) - \varepsilon\nu(A)}{1 - \varepsilon}, $$

and introduce the coupling kernel $\bar{P}_k$ on $X^2$, defined for all $A \times A' \subset X^2$ and all $z_k = (x_k, x'_k) \in X^2$ by

$$ \bar{P}_k(z_k, A \times A') = \mathbf{1}_{x_k=x'_k} P_k(x_k, A) \delta_{x_{k+1}}(A') + \mathbf{1}_{x_k \neq x'_k} \mathbf{1}_{z_k \notin \mathcal{K}^2} P_k(x_k, A) P_k(x'_k, A') \\ + \mathbf{1}_{x_k \neq x'_k} \mathbf{1}_{z_k \in \mathcal{K}^2} (\varepsilon\nu(A) \delta_{x_{k+1}}(A') + (1-\varepsilon) Q_k(x_k, A) Q_k(x'_k, A')). \quad (13) $$

**Lemma D.2** Let $(\bar{P}_k)_{k \ge 0}$ be the Markov kernels on $(X^2, \mathcal{X}^2)$ defined by (13). Then, for all $n \in \mathbb{N}$ and all $(x, x') \in X^2$, we have

$$ \bar{P}_1 \cdots \bar{P}_m((x, x'), \cdot) \in \mathcal{C}(P_1 \cdots P_m(x, \cdot), P_1 \cdots P_m(x', \cdot)), $$

where $\mathcal{C}$ denotes the set of couplings introduced in Definition B.3.

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

$$ \bar{P}_1((x, x'), \cdot) \in \mathcal{C}(P_1(x, \cdot), P_1(x', \cdot)). $$

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

$$ \bar{P}_1 \cdots \bar{P}_m((x, x'), \cdot) \in \mathcal{C}(P_1 \cdots P_m(x, \cdot), P_1 \cdots P_m(x', \cdot)). $$
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}$,

$$ 
\begin{align*}
\bar{P}_1 \cdots \bar{P}_{m+1}((x, x'), A \times X) &= \int_{X \times X} \bar{P}_1 \cdots \bar{P}_m((x, x'), dy dy') \bar{P}_{m+1}((y, y'), A \times X) \\
&= \int_{X \times X} \bar{P}_1 \cdots \bar{P}_m((x, x'), dy dy') P_{m+1}(y, A) \\
&= \int_{X \times X} P_1 \cdots P_m(x, dy) P_{m+1}(y, A) \\
&= P_1 \cdots P_{m+1}(x, A),
\end{align*}
$$

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

$$ \bar{P}_1 \cdots \bar{P}_{m+1}((x, x'), X \times A) = P_1 \cdots P_{m+1}(x', A). $$

Thus, we conclude that

$$ \bar{P}_1 \cdots \bar{P}_{m+1}((x, x'), \cdot) \in \mathcal{C}(P_1 \cdots P_{m+1}(x, \cdot), P_1 \cdots P_{m+1}(x', \cdot)). $$

This completes the induction and proof. □

**Lemma D.3** For all $(x, x') \in X^2$ define $\Delta(x, x') = \mathbf{1}_{x \neq x'}$ and $\bar{V}(x, x') = (V(x) + V(x'))/2$. Then, for all $k \in \mathbb{N}$:

* If $(x, x') \in \mathcal{K}^2$, then

$$ \bar{P}_k \Delta(x, x') \leq (1 - \varepsilon) \Delta(x, x'), \quad \bar{P}_k \bar{V}(x, x') \leq b. $$

* If $(x, x') \notin \mathcal{K}^2$, then

$$ \bar{P}_k \Delta(x, x') \leq \Delta(x, x'), \quad \bar{P}_k \bar{V}(x, x') \leq \lambda \bar{V}(x, x'). $$

*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') = \frac{P_k V(x) + P_k V(x')}{2}. $$

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

$$ \frac{P_k V(x) + P_k V(x')}{2} \leq b. $$

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

$$ \bar{P}_k \bar{V}(x, x') = \frac{P_k V(x) + P_k V(x')}{2} \leq \frac{\lambda V(x) + \lambda V(x')}{2} = \lambda \bar{V}(x, x'). $$

This concludes the proof. □
Proof of Proposition D.1. For any $t \in (0, 1)$, we define

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

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

$$ \begin{aligned} \|P_1 \cdots P_m(x, \cdot) - P_1 \cdots P_m(x', \cdot)\|_{\text{TV}} &= 2 \inf \{\zeta(\Delta) : \zeta \in \mathcal{C}(P_1 \cdots P_m(x, \cdot), P_1 \cdots P_m(x', \cdot))\} \\ &\le 2\bar{P}_1 \cdots \bar{P}_m \Delta(x, x') \\ &\le 2\bar{P}_1 \cdots \bar{P}_m W(x, x'). \end{aligned} $$

where we used $V \ge 1$. Finally, applying Hölder's inequality and using Lemma D.3, we obtain, for all $(x, x') \in \mathcal{X}^2$,

$$ \begin{aligned} \bar{P}_k W(x, x') &= \bar{P}_k (\Delta^{1-t} \bar{V}^t)(x, x') \le (\bar{P}_k \Delta(x, x'))^{1-t} (\bar{P}_k \bar{V}(x, x'))^t \\ &\le \Delta^{1-t}(x, x') \times \begin{cases} (1-\varepsilon)^{1-t} b^t & \text{if } (x, x') \in \mathcal{K}^2 \\ \lambda^t \bar{V}^t(x, x') & \text{if } (x, x') \notin \mathcal{K}^2 \end{cases} \\ &\le \varrho_t W(x, x'). \end{aligned} $$

This implies by induction that for all $m \in \mathbb{N}$ and all $(x, x') \in \mathcal{X}^2$,

$$ \bar{P}_1 \cdots \bar{P}_m W(x, x') \le \varrho_t^m W(x, x'). $$

Then,

$$ \begin{aligned} \|P_1 \cdots P_m(x, \cdot) - P_1 \cdots P_m(x', \cdot)\|_{\text{TV}} &\le 2\varrho_t^m W(x, x') \\ &\le \begin{cases} 2\varrho_t^m & \text{if } x \in \mathcal{K}, \\ \varrho_t^m (V(x) + V(x')) & \text{if } x \notin \mathcal{K}. \end{cases} \end{aligned} $$

This concludes the proof. $\square$

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

*Proof.* Using Lemma D.1, and taking $s = m - b_m$, we have:

$$ \begin{aligned} \|\mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_{\infty}}\|_{\text{TV}} &\le \|\mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_{m-b_m}} P_{\theta_{m-b_m}} \cdots P_{\theta_m}\|_{\text{TV}} \\ &\quad + \sum_{j=m-b_m}^{m-1} \|\pi_{\theta_{j+1}} - \pi_{\theta_j}\|_{\text{TV}} + \|\pi_{\theta_m} - \pi_{\theta_{\infty}}\|_{\text{TV}} \\ &\le \|\mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_{m-b_m}} P_{\theta_{m-b_m}} \cdots P_{\theta_m}\|_{\text{TV}} \\ &\quad + L(x) \sum_{j=m-b_m}^{m-1} \|\theta_{j+1} - \theta_j\| + L(x) \|\theta_m - \theta_{\infty}\|, \end{aligned} $$

where we used the Lipschitz condition of $\pi_{\theta}$. For the first term, using Proposition D.1 with $\xi = \mu P_{\theta_1} \cdots P_{\theta_{m-b_m-1}}$ and $\xi' = \pi_{\theta_{m-b_m}}$, we have:

$$ \|\mu P_{\theta_1} \cdots P_{\theta_m} - \pi_{\theta_{m-b_m}} P_{\theta_{m-b_m}} \cdots P_{\theta_m}\|_{\text{TV}} \le \kappa \rho^{b_m}. $$
For the second term, using the Lipschitz condition (Assumption 7) and the recursion of $\theta_{j+1}$, we get:

$$ 
\begin{aligned}
\sum_{j=m-b_m}^{m-1} \mathbb{E}[\|\theta_{j+1} - \theta_j\|] &= \sum_{j=m-b_m}^{m-1} \gamma_{j+1} \mathbb{E}[\|H_{\theta_m}(x_{m+1})\|] \\
&= \mathbb{E}[L(x)] \sum_{j=m-b_m}^{m-1} \gamma_{j+1}.
\end{aligned}
$$

For the last term, Using Jensen inequality and Assumption 8, we obtain:

$$ \mathbb{E}[\|\theta_m - \theta_\infty\|] \le \|\theta_m - \theta_\infty\|_{L_2}^2 = \mathcal{O}(a_m). $$

□

## D.2 Proof of Theorem 5.2

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

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

Following the procedure outlined in the proof of Theorem 5.1, it is straightforward to verify Assumption 6 (the minoration and drift condition) with $V \ge 1$. Additionally, using Assumptions 3 and 4, we can verify Assumptions 7 and 8. The proof is then concluded by applying Theorem D.1.

## 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) is bounded under Assumption 3, 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 can be satisfied.
# 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 $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) are chosen according to an exponential decay schedule. The training loss and corresponding cost are shown in Figure 4.

<figure><img src="image_5.png" /><figcaption>Figure 4: Training loss and score for our model trained on CVRP instances with nodes $n = 100$</figcaption></figure>

## 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 and 2, Obj. denotes the value of the cost function defined in (5) for the TSP and CVRP. The Gap indicates the percentage gap to optimality, computed as:

$$ \text{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 \le k \le K$,

$$q(z_{m+1}^k | z_m^{1:K}) = \mathcal{N}(z_{m+1}^k; z_m^k + \gamma(z_m^{I_1} - z_m^{I_2}), \sigma^2 I_{d_z})$$

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.

## E.2.2 Illustration of hyperparameters

Figure 5 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><img src="image_6.png" /><figcaption>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.</figcaption></figure>

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, 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><img src="image_7.png" /><figcaption>Figure 6: Average gap to optimality for different values of *K* in the CVRP with *n* = 100</figcaption></figure>

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><img src="image_8.png" /><figcaption>Figure 7: Performance comparison of various SA step update frequencies on CVRP with $n = 100$</figcaption></figure>

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