Billiard Walk (BW) Algorithm

• Billiard Walk (BW) is an MCMC sampler that generates uniform samples in bounded, open, convex sets using billiard-style trajectories with specular reflections.
• It combines random direction sampling, exponentially distributed free-flight lengths, and deterministic reflection dynamics to improve mixing and efficiency over traditional methods like Hit-and-Run.
• Empirical studies show that BW reduces reflection counts and autocorrelation in high-dimensional and geometrically complex domains, highlighting its potential for efficient uniform sampling.

The Billiard Walk (BW) algorithm is a Markov chain Monte Carlo (MCMC) sampler for generating uniformly distributed samples from a bounded, open, connected, convex domain $K \subset \mathbb{R}^d$ with piecewise-smooth boundary. Distinguished by its use of billiard-style trajectories with specular reflections, BW aims to accelerate mixing—particularly in high-dimensional or geometrically challenging settings—relative to classical algorithms such as Hit-and-Run. The method’s core innovation lies in combining independent random directions, exponentially distributed free-flight lengths, and deterministic reflection dynamics to traverse the sample space efficiently while maintaining uniformity and ergodicity of the stationary law (Gryazina et al., 2012).

1. Mathematical Formulation

Let $K \subset \mathbb{R}^d$ be a bounded, open, connected, convex set with piecewise-smooth boundary $\partial K$. The sampling target is the uniform measure on $K$ with density

$\pi(x) = \frac{1}{\mathrm{Vol}(K)},\quad x \in K.$

BW operates as follows:

• Direction sampling: From the current state $x_i \in K$, draw direction $d \sim \mathrm{Unif}(S^{d-1})$. Typically, $d = \xi/\|\xi\|$ where $\xi \sim \mathcal{N}(0, I_d)$.
• Random trajectory length: Draw $U \sim \mathrm{Unif}(0,1)$, set $\ell = -\tau \ln U$ for $\tau > 0$, so $\ell$ has density $f_\tau(\ell) = \frac{1}{\tau} e^{-\ell/\tau}$ for $\ell \geq 0$.
• Billiard trajectory with reflections: Starting from $x_i$ in direction $d$, move at unit speed until either (a) distance $\ell$ is exhausted, or (b) $\partial K$ is hit. On collision at $y \in \partial K$ with inward normal $s(y)$, if remaining length allows, the direction is updated via the reflection law:

$d \longmapsto d - 2\langle d, s(y) \rangle s(y)$

and the trajectory continues. The walk terminates when the total distance $\ell$ is achieved; the endpoint becomes $x_{i+1}$.

• Safeguards: If more than $R$ reflections are incurred before exhausting $\ell$ or a nonsmooth corner is encountered, the step is rejected and restarted with a new direction.

2. Algorithmic Specification

The BW algorithm admits the following pseudocode structure:

$$\begin{aligned} &\textbf{Input: } x^0 \in K,\ \tau > 0,\ R \in \mathbb{N}.\ &\textbf{For } i=0,1,2,\dots,N-1\text{:}\ &\quad\text{1. Sample } U \sim \mathrm{Unif}(0,1),\quad \ell \gets -\tau\ln U.\ &\quad\text{2. Sample direction } d \in S^{d-1} \text{ uniformly.}\ &\quad\text{3. Set } y \gets x^i,\, r \gets \ell,\, k \gets 0.\ &\quad\textbf{while } r>0 \text{ and } k<R\textbf{ do}\ &\qquad(a)\ \text{Compute}\ t^* = \min\{ t>0 : y + t\,d \in \partial K \}.\ &\qquad(b)\ \text{If } t^* \geq r \text{ then } y \gets y + r\,d,\, r \gets 0.\ &\qquad\text{else}\ &\qquad\quad y \gets y + t^* d,\, r \gets r - t^*,\, k \gets k+1.\ &\qquad\quad \text{Let } s \text{ be the inward normal at } y.\ &\qquad\quad d \gets d - 2(d\cdot s)s.\ &\quad\textbf{end while}\ &\quad\text{If } r > 0 \text{ (too many reflections or corner), reject and repeat from step 2.}\ &\quad x^{i+1} \gets y. \end{aligned}$$

The transition law in each step is symmetric and reversible, essential for guaranteeing convergence to the uniform distribution.

3. Stationarity and Mixing-Time Analysis

The algorithm’s design ensures that the uniform law is the unique stationary distribution. The essential properties enabling this are:

• Symmetry: Each step’s direction $d$ is sampled uniformly from $S^{d-1}$.
• Reversibility: The billiard dynamics, including the reflection law, are time-reversible.
• Positivity: The one-step transition density $p(y|x) > 0$ for all $x,y \in K$.

Thus, as $n \to \infty$, the chain’s distribution converges in total variation to $\pi$ [(Gryazina et al., 2012), Theorem 1].

Mixing-time bounds are benchmarked against Hit-and-Run (HR), for which it is known that if $K$ has inradius $r$ and outradius $R$, then the mixing time obeys

$$\tau_{\rm HR}(\varepsilon) = O\bigl(d^3 (R/r)^2 \ln (1/\varepsilon)\bigr).$$

For BW, no explicit polynomial bound is established, but a plausible conjecture, reflecting the efficiency of ballistic segments, is

$$\tau_{\rm BW}(\varepsilon) = O\bigl(d^2 (R/r)^2 \ln(1/\varepsilon)\bigr).$$

This suggests potential for improved scaling in dimension and path efficiency, though rigorous proofs are pending.

4. Empirical Performance and Comparisons

Numerical experiments systematically compare BW and HR under a fixed “boundary-oracle” (BO) call budget, which counts line–$\partial K$ intersection queries:

• Plane angle $\alpha$: BW guarantees exit after at most $\lceil \pi/\alpha \rceil$ reflections with probability $1$; HR’s exit probability after $n$ steps is $1-(1-\alpha/\pi)^n$. Empirically, BW averages $\approx (\pi/2\alpha)$ reflections vs HR’s $\approx (\pi/\alpha)$ steps.
• Nonnegative orthant $\{x>0\} \subset \mathbb{R}^n$: BW exits in at most $d$ reflections, while HR requires approximately $2^{n-1}$ steps to exit with constant probability.
• Strip domain $\{0<x_2<1,\, |x_1|<M\}$: BW moves along $x_1$ about six times faster per BO call than HR.
• Unit cube $[0,1]^n$: For $n=10,25,50$ under equal BO budgets, BW yields markedly lower serial correlation and passes $\chi^2$ fit tests more frequently.
• High-dimensional simplex and toroid: In cases such as $n=50$ simplex and $n=10$ toroid, empirical CDFs and 2D projections show that BW samples adhere much more closely to the uniform law.

The ballistic movement with reflection enhances exploration in domains with narrow passages or sharp corners, areas where HR shows slow progress.

5. Assumptions, Parameterization, and Limitations

BW requires that $K$ be open, bounded, connected, and possess a piecewise-smooth boundary so that normals are defined almost everywhere. The algorithm assumes:

• Well-defined tangent planes at collision points for the reflection law (nonsmooth, “cornery” points have measure zero or trigger rejection and restart).
• Trajectory parameter $\tau$ chosen $\approx \mathrm{diam}(K)$, interpolating between local walk behavior for $\tau \ll\mathrm{diam}(K)$ and excessive reflection for $\tau \gg\mathrm{diam}(K)$.
• Reflection limit $R \approx 10d$ is suggested to prevent rare, pathological infinite-bounce situations.
• No explicit spectral or conductance bounds yet exist for BW; derivation of polynomial mixing-time bounds remains open.

A summary is given in the table below:

Assumption / Parameter	Required Value	Practical Guidance
$K$	Open, bounded, connected, piecewise-smooth boundary
Reflection law	Well-defined normal at collision	Corners handled by rejection
$\tau$	$\sim \mathrm{diam}(K)$	Interpolates locality vs. globality
$R$	$\sim 10d$	Prevents infinite reflection loops

6. Applications and Practical Considerations

BW’s most pronounced advantages appear in high-dimensional, sharply bounded, or geometrically intricate sets—domains where rapid movement along straight-line segments and exploitation of reflections yields superior mixing and lower autocorrelation relative to classical HR. Areas with sharp “corners” or thin “tunnels” see especially notable improvements, as confirmed by empirical studies. The uniform ergodicity and reversibility properties ensure robust integration with established MCMC and convex analysis frameworks. A plausible implication is that, subject to theoretical refinements, the algorithm could set a new benchmark for high-dimensional random sampling, once rigorous mixing-time bounds are established (Gryazina et al., 2012).