Gaussian Approximation in Wasserstein-p: Markovian Martingale Differences and Stochastic Approximation

Abstract

We establish nonasymptotic Gaussian approximation bounds in higher-order Wasserstein distance $W_p$, $p\geq2$, for deterministic matrix-weighted sums of multivariate Poisson martingale differences generated by uniformly ergodic Markov chains. Under an $L^{(2+\eta)p}$-moment condition, $\eta>0$, we obtain the explicit bound $$O\left(p^3\|A\|_4^2+pd^{1/4}\|A\|_2^{1/2}\|A\|_4^2\right),$$ where $A\in\mathbb{R}^n$ collects the moment sizes of the weighted increments. When these sizes are of order $n^{-1/2}$, the bound achieves the first optimal $O(n^{-1/2})$ rate for fixed $p$ and $d$. Consequently, we also obtain the first optimal $O(n^{-1/2})$ $W_p$ Gaussian approximation rate for multivariate additive functionals of uniformly ergodic Markov chains. Our analysis develops two techniques for addressing the interplay between higher-order Wasserstein distance and temporal dependence. First, building on the Ornstein--Uhlenbeck relative-score approach of Fang and Koike (2023), we formulate the bound in terms of antisymmetric Stein couplings while retaining the conditional tensor structure. Second, we develop a refresh-then-maximal coupling that combines an independent first-step resampling, which preserves the desired Stein identity, with a subsequent maximal coupling that provides effective control of the coupling increment. We apply these bounds to nonlinear Markovian stochastic approximation by decomposing normalized errors into matrix-weighted Poisson martingale sums and remainders. This yields the first sharp $W_p$ Gaussian approximation rates for both last iterates and Polyak--Ruppert averages. In particular, the optimized rates are $O(n^{-1/3})$ for polynomial-stepsize last iterates and $O(n^{-1/4})$ for Polyak--Ruppert averages, with an $O(n^{-1/2})$ rate for harmonic-stepsize last iterates under an additional spectral condition.