Beyond Self-Repellent Kernels: History-Driven Target Towards Efficient Nonlinear MCMC on General Graphs

Q1. Downstream benefits in ML algorithms

Our HDT framework enables more efficient sampling of large or complex graphs for arbitrary target measures, reducing estimator variance and cost. Key impacts include:

1. Estimating Global Properties of Large-Scale Graphs: In social networks, e-commerce systems, recommender systems and other domains, practitioners need to estimate graph size, node degree distribution, vertex label distribution [Xie et al. 2023]. These applications include detecting bot populations, identifying high-value user segments or scarce features [Nakajima et al. 2023], and gradient of local loss for token algorithms in distributed optimizations [Sun et al. 2018, Even 2023], especially with limited graph access. However, a standard random walk can be slow to discover underrepresented regions or rare features, forcing more samples to achieve acceptable levels of accuracy.
2. Sampling from Large-Scale Graphs: In statistical physics (e.g., Ising or Potts models) [Grathwohl et al. 2021, Zhang et al. 2022], machine learning (e.g., Restricted Boltzmann machines / energy-based models for image generation, reconstruction or partition, Bayesian variable selection for high-dimensional regression/classification) [Hinton 2010, Zanella 2020, Rhodes et al. 2022, Chang et al. 2022, Sun et al. 2023], and econometrics (e.g., MCMC-MLE method for exponential random graph models) [Bhamidi et al. 2011, Byshkin et al. 2016], the general discrete state space can always be translated into a graph where each node represents a possible configuration and edges link configurations that differ by the Hamming distance. Efficient exploration yields quicker, more accurate system representation and inference.

HDT-MCMC strategically modifies the target distribution at minimal extra cost compared to baseline MCMC and SRRW, achieving equivalent accuracy with fewer samples and computations. While fully deployment in exponentially large configuration state spaces present ongoing memory challenges beyond this paper's scope, our method offers versatile and simple implementation alongside significant efficiency improvements for sampling on general graphs. We plan to detail these applications in a revised appendix if our paper is accepted.

[References not included in our manuscript are listed below for this response.]

G. Hinton. "A practical guide to training restricted Boltzmann machines." Momentum, 2010.

S. Bhamidi et al. "Mixing time of exponential random graphs." Ann. Appl. Probab., 2011.

M. Byshkin et al. "Auxiliary parameter MCMC for exponential random graph models." J. Stat. Phys., 2016.

T. Sun et al. "On markov chain gradient descent." NeurIPS, 2018.

H. Xie et al. "A bootstrapping approach to optimize random walk based statistical estimation over graphs." IEEE TKDE, 2023.

K. Nakajima et al. "Random walk sampling in social networks involving private nodes." ACM TKDD, 2023.

M. Even. "Stochastic gradient descent under markovian sampling schemes." ICML, 2023.

Q2. Impact of graph properties on the gain of HDT-MCMC

The gains of HDT-MCMC relative to baseline MCMC and SRRW are addressed as follows:

1. HDT-MCMC v.s. Baseline MCMC: Our Theorem 3.3 shows HDT-MCMC reduces sampling variance by a factor of $1/(2α+1)$ compared to baseline MCMC:

$V^{base}$ is defined via transition kernel’s spectrum in Eq. (4) in Line 215 and is also related to spectral gap and mean hitting time [Aldous et al. 2014, Remark after Theorem 2.17 & Eq. (11.4)], where graph with higher connectivity usually ensures smaller spectral gap and mean hitting time. However, HDT-MCMC uniformly reduces variance by the $1/(2α+1)$ factor. Corollary 3.4 adds that a more efficient baseline amplifies this gain: the better the baseline, the more HDT lowers the variance, regardless of the graph or target measure. 2. HDT-MCMC v.s. SRRW: Our Lemma 3.6 highlights HDT-MCMC’s cost-based variance reduces by a factor of $\\frac{2}{\\mathbb{E}\_{i∼μ} [|\\bar{\\mathcal{N}}(i)|]}$ compared to SRRW. The average degree $\\mathbb{E}\_{i∼μ} [|\\bar{\\mathcal{N}}(i)|]$ (including self-loops) depends on graph density and the target measure $\mu$. Even for a sparse graph with an average degree of $3$, the factor is smaller than $0.5$, translated into over 50% variance reduction compared to SRRW under limited budgets. Gains are more pronounced in denser graphs (higher average degree under $\mu$), making HDT particularly effective for large, high-degree networks, as shown numerically in Figure 3 and around Lines 418 - 423.

D. Aldous et al. "Reversible Markov chains and random walks on graphs." Unfinished monograph, 2014.

Additional experiments on non-uniform target

While we use the uniform target in our paper, we've added experiments targeting degree distribution on more graphs to address this. Due to space limits, these results are in the PDF with a link in the response to reviewer mrvq.

Thank you for your constructive and practical feedback. We agree that demonstrating robustness to initialization and quantifying uncertainty in our results will enhance the paper. Below is the setting of additional experiments. Due to space constraints, we defer the simulation results on more real-world graphs under TVD and NRMSE metrics in the anonymous linked PDF:

https://www.dropbox.com/scl/fi/ujmodl9e8721nef5dx2k5/ICML_2025_Rebuttal_Additional_Simulations.pdf?rlkey=80cp2rhdvf59cvjukrtrcz5km&st=b96ipxpx&dl=0

Robustness of HDT-MCMC to different initialization strategies and uncertainty quantification.

We run $1000$ independent trials for different MCMC methods on the wikiVote graph with uniform target and show the TVD $\mathbb{E}[||x_n-\mu||_1]$ at $n=3000$ steps with 95% confidence intervals:

1. Varying the Initial State $X_0$:
   Besides random selection, we started from nodes in different regions of the graph:

   • Low-degree group (nodes with below-average degree)
   • High-degree group (nodes with above-average degree)

   ---

   Table 1. Mean (std error) TVD with different $X_0$ initializations. All std error values are in units of $10^{-4}$.

   Method	Low-Deg Group	High-Deg Group
   MHRW	0.573 (7.316)	0.573 (7.290)
   HDT-MHRW	0.417 (3.785)	0.417 (3.804)
   MTM	0.533 (6.818)	0.533 (7.076)
   HDT-MTM	0.331 (5.417)	0.330 (4.154)
   MHDA	0.544 (7.419)	0.543 (7.185)
   HDT-MHDA	0.395 (3.564)	0.395 (3.569)

   ---

2. Varying the Fake Visit Count $x_0$:
   We use three scenarios:

   • 'Deg' scenario: Fake visit counts set by node degrees, heavily favoring high-degree nodes.
   • 'Non-unif' scenario: Fake visit counts randomly drawn from a Dirichlet distribution with default concentration parameter of $0.5$.
   • 'Unif' scenario: Same initial count for all nodes.

   In this part, we start with $X_0$ uniformly at random. Note that baseline methods (MHRW, MTM, MHDA) do not depend on fake visit counts $x_0$, so they give identical TVD values in all three scenarios.

   ---

   Table 2. Mean (std error) TVD with different $x_0$ initializations. All std error values are in units of $10^{-4}$.

   Method	Deg	Non-unif	Unif
   MHRW	0.573 (7.085)	0.573 (7.085)	0.573 (7.085)
   HDT-MHRW	0.416 (3.767)	0.416 (3.877)	0.416 (3.810)
   MTM	0.532 (7.092)	0.532 (7.092)	0.532 (7.092)
   HDT-MTM	0.325 (4.301)	0.332 (4.200)	0.330 (4.130)
   MHDA	0.543 (7.219)	0.543 (7.219)	0.543 (7.219)
   HDT-MHDA	0.395 (3.605)	0.395 (3.657)	0.396 (3.587)

   ---

In both Tables 1 and 2, HDT-MCMC shows robust convergence and consistent variance reduction relative to the baseline. We observe that the initialization effect quickly vanishes during the burn-in period (one third of the total samples) so they do not affect the performance of HDT-MCMC. Even when the initial visitation $x_0$ is heavily skewed (such as degree-based distribution where weight goes mostly to the node with largest degree), the HDT forces rapid exploration of under-sampled regions.

In most cases, the standard error for our HDT-MCMC framework is smaller than baseline algorithms, indicating the consistency of HDT-MCMC across repeated runs. This shows that every sample path generated by our HDT framework performs consistently. These results suggest that the HDT $π[x]$ contributes to stable sampling trajectories that rapidly converge to the target distribution $μ$. We conjecture that this consistent performance of HDT-MCMC is due to the correctness behavior from our design of history-aware scheme in every trajectory to drive $x_n$ closer to the target distribution $μ$.

We will include detailed simulation setups and plots in the revised paper to illustrate these findings if our paper is accepted.

Thank you for your thoughtful and detailed review. We appreciate your positive remarks on our paper’s theoretical depth and are grateful for the questions you raised. Below, we address each of your points in turn.

Response to Q1

We agree that this point warrants further explanation. Here is a streamlined proof by contradiction:

1. Assume that $g(c_1,c_2,x,\\mu)$ depends on $x, \\mu$.
2. According to Line 650, for any coordinate $i \\in [N]$

where the function $f :\\mathbb{R}\_{>0}×\\mathbb{R}\_{>0}→\\mathbb{R}\_{>0}$. 3. Consider two pairs $(x,μ)$ and $(\\hat x,\\hat μ)$ whose $i$-th coordinate are identical (i.e., $x_i=\\hat x_i,μ_i=\\hat μ_i$). Then, $f(c_1 x_i,c_2 μ_i )=f(c_1 \\hat x_i,c_2 \\hat μ_i )$ and $f(x_i,μ_i )=f(\\hat x_i,\\hat μ_i )$. By the equation in step 2, it follows that

4. Similarly, construct another pair ($\\tilde x,\\tilde μ$) (i.e., $\\tilde x_j=\\hat x_j,\\tilde μ_j=\\hat μ_j$ for $j$-th coordinate) such that $g(c_1,c_2,\\tilde x, \\tilde μ)=g(c_1,c_2,\\hat x,\\hat μ)$. Then, $(\\hat x, \\hat μ)$ serves as a ‘bridge’ to $(x,μ)$ and $(\\tilde x,\\tilde μ)$, leading to

For a given pair $(x,μ)$, we can always construct a distinct pair $(\\tilde x,\\tilde μ)$ such that $g(c_1,c_2,x,μ)=g(c_1,c_2,\\tilde x, \\tilde μ)$, indicating that function $g$ must be identical for any $(x,μ)$. Hence $g$ cannot actually depend on $x$ or $μ$; it must be a function of $(c_1,c_2)$ alone.

We will revise Lines 653–655 in the manuscript with this explanation to clarify why $g$ must be independent of $(x,μ)$.

Response to Q2

Thank you for catching this. We will correct the mentioned typo and any others in the revision.

Response to Q3

This is an excellent question. You are correct that Appendix E.1 and E.2 do not assume reversibility, allowing non-reversible chains to operate directly on the state space $\\mathcal{X}$ of interest. Let $\\mathcal{Y}$ be the auxiliary state space. Below, we explain why analyzing non-reversible chains on the augmented state space $\\mathcal{Z}≜\\mathcal{X}×\\mathcal{Y}$ is more than $P[x]$ preserving $π[x]$, and discuss the technical challenges in directly applying the stochastic approximation (SA) framework, as well as how Appendix E.3 addresses them.

1. State Space Inconsistency in SA Analysis
   The SA iteration

requires $\\{X_n\\}$ to be controlled Markovian dynamics. That is, $X_{n+1}$ is drawn from a transition kernel $P_{(X_n,⋅)} [x_n]$ parameterized by empirical measure $x_n$ (in Line 891). If the chain instead operates on an augmented space $\\mathcal{Z}$, then $\\{X_n\\}$ in isolation cannot be viewed as a controlled Markov chain, since $X_{n+1}$ depends not only on $X_n, x_n$, but also on the auxiliary state $Y_n\\in \\mathcal{Y}$.
More concretely, mixing these spaces (using a transition matrix on $\\mathcal{Z}$ but a function $H(x,⋅)$ only defined on $\\mathcal{X}$) causes a dimensional mismatch in the SA analysis. For example, in Eq. (32) around Line 952, the solution

to the Poisson equation is not well-defined because $P[x]$ is specified on $\\mathcal{Z}$, whereas $H(x,⋅)$ is only on $\\mathcal{X}$. From the viewpoint of SA iteration, a direct redefinition of $H$ via

does not resolve the core mismatch because $δ_Z$ (a canonical vector in $\\mathcal{Z}$) does not align with the dimension of $x$ (an empirical measure in $\\mathcal{X}$).

2. Dimension and Memory Overhead
   A natural solution to the above mismatch is to record visit frequencies for the augmented space $\\mathcal{Z}$, thereby resolving dimensional mismatch but inflating the dimension of $x$. For instance, this would increase the dimension of $x$ from $|\\mathcal{X}|$ to $|\\mathcal{X}|^2$ for MHDA (Lee et al., 2012) and to $2|\\mathcal{X}|$ for the 2-cycle Markov chain (Maire et al., 2014). This is especially problematic for large state space $\\mathcal{X}$.

Our “Trick” in Appendix E.3
We address these issues by replacing $H$ with

so that the augmented state $Z∈\\mathcal{Z}$ is 'accepted' by $\phi$, but we only store visit frequencies for $\\mathcal{X}$ (i.e., we project $\\mathcal{Z}$ onto $\\mathcal{X}$). This reconciles the SA analysis with the non-reversible Markov chain on $\\mathcal{Z}$ (via the function $\phi$) without inflating the dimension of $x$.

We will revise Line 1043 in the manuscript to clearly emphasize how this approach extends our SA analysis to non-reversible chains on augmented state spaces, mitigating both dimensional mismatch and memory concerns.