2026-10-14 | | Total: 20
Sampling from high-dimensional, heavy-tailed distributions is a fundamental challenge in computational statistics, as many standard Markov chain Monte Carlo (MCMC) methods mix poorly in such settings. Recently, Stereographic MCMC [Yang et al., 2024] and the Sub-Cauchy Projection Sampler [Grazzi et al., 2026] have been shown to perform successfully on such tasks. However, establishing their non-asymptotic convergence properties remains an important open problem. In this work, we fill this gap by establishing non-asymptotic upper bounds on the mixing time of the stereographic projection and sub-Cauchy projection samplers. Our results demonstrate that, under certain conditions on the target and initial distributions, the mixing time is polynomial in dimension for a broad class of distributions, including light- and heavy-tailed cases. Motivated by the theoretical analysis, we further establish a new weighted isoperimetric inequality that extends the classical version for log-concave distributions to the heavy-tailed setting with optimal dimension dependence. The proof techniques provide new insights into the geometric properties of heavy-tailed distributions that govern rapid mixing in high dimensions.
Coordinate descent solvers for regularized generalized linear models (GLMs) differ in how they update the intercept. This seemingly routine choice can substantially affect convergence when the response is imbalanced. For direct coordinate descent on the original GLM loss, we compare a conservative step based on worst-case curvature, a safeguarded Newton step based on current curvature, and repeated Newton steps that solve the intercept subproblem. The conservative update can leave much of the interaction between the intercept and the coefficients unresolved, producing a slowdown when that interaction controls convergence. A single Newton step avoids this attenuation and, in our experiments, follows nearly the same outer trajectory as exact minimization while requiring fewer inner steps. Simulated and real-data experiments, including comparisons across six production solvers and a controlled intervention within skglm, support this explanation. We also show how the same mechanism appears in quadratic-approximation methods and extends to multinomial models with rare classes. We recommend a safeguarded Newton update for direct coordinate descent and full intercept optimization within frozen local quadratic surrogates.
Bridge sampling estimates of marginal likelihoods are routinely reported with a delta-method Monte Carlo standard error (MCSE). Existing theory is asymptotic and cannot say whether a given run is long enough to trust the estimate or its MCSE. We give finite-sample answers. The optimal bridge estimator fits a logistic regression that separates posterior from proposal draws, and its accuracy is governed by their overlap, within a factor of two of the classification's Bayes error. Because the estimating function is a sum of bounded terms, no moment conditions are needed, and all bounds hold uniformly over targets, proposals and dimensions. A number of draws inversely proportional to the overlap suffices and is necessary in the worst case, paralleling the exponential-in-Kullback--Leibler requirement of importance sampling. A Berry--Esseen bound with an explicit constant holds at the optimal rate, and the MCSE is accurate once the effective overlap is large; it is estimated at no cost. For Markov chain draws, the required sample size grows in proportion to the sampler's relaxation time, unavoidably in the worst case. The estimated bridge, the fixed-point iteration and the fitted proposal matter only at second order, and we give free diagnostics for sampler failure and misestimated autocorrelation.
We present CLASSO, a path-following method for the LASSO subject to general linear equality and inequality constraints. The local direction is obtained from a quadratic programme on an augmented active set, which determines both the coefficients that leave zero and the weakly active inequalities that remain binding. For a positive-definite Hessian and a nonempty feasible set, this programme has a unique primal minimiser, even when the binding constraint rows are dependent, and gives the local derivative of the solution path. The construction recovers the LASSO-modified LARS direction without constraints, and over the full cone gives the segment leaving the origin. We prove that the primal path is piecewise affine with finitely many kinks, without a constraint qualification or unique multipliers. A finite path does not imply a finite traversal: re-certifying at each kink can fail to terminate. CLASSO instead follows the minimum-norm certificate along a lexicographic rate; we prove that this rule reaches $λ=0$ in exact arithmetic within an explicit bound and traverses the whole path. These results are machine-checked in Lean~4, and the floating-point implementation agrees with an exact rational one on all $210$ small test problems. We assess the implementation on a $39$-case synthetic suite and six constrained versions of the diabetes benchmark. On the $31$ nontrivial synthetic cases with $n\ge p$, sampled relative objective differences from numerical reference solves are below $10^{-9}$ in magnitude; unnormalised endpoint differences reach $1.5\times10^{-7}$ in magnitude. For $n<p$, matched differences range from $-3.6\times10^{-7}$ to $2\times10^{-9}$ across five ridge strengths. Across the synthetic and diabetes benchmarks, checks of all $580$ recorded segments found no missed interior event and all checked KKT residuals within the stated tolerances.
The multivariate linear mixed model combined with restricted maximum likelihood estimation of (co)variance components and solving the mixed model equations is widely used for the analysis of plant breeding data. Several algorithms for restricted maximum likelihood estimation have been proposed, with the average information algorithm being widely used for its computational efficiency, quadratic rate of convergence, and ability to handle large amounts of parameters. A key reason for the scalability of linear mixed models is the exploitation of sparse matrix techniques. Recently, decreasing costs of high-throughput data collection have led to high-dimensional sets of features characterizing genotypes and environments. A convenient way of using these features is through the use of kernel-based covariance structures. However, while convenient, these kernels typically have dense inverses meaning optimizations based on sparse matrix techniques are not always applicable. For these settings, the introduction of general-purpose computing on graphics processing units brought significant computational advantages for dense linear algebra operations. We present a hybrid GPU/CPU implementation of the average information restricted maximum likelihood algorithm. We eliminate sparse equations using supernodal Cholesky factorization and rely on CUDA GPU acceleration to solve the resulting reduced dense system. We show the significant computational advantages of this approach for fully dense problems, as well as problems with a mix of dense and sparse equations.
Exact reversible locally adaptive step sizes are introduced by Bou-Rabee et al. (2026) in the context of the No-U-Turn Hamiltonian Monte Carlo (HMC) sampler for the purpose of high performance HMC sampling in challenging target distributions. This paper carefully explains and illustrates the locally adaptive step size procedure under more general assumptions. Alternative error measures are introduced with the aim of further enhancing MCMC performance by reducing the prevalence of reversibility failure rejections. In addition, the paper shows how to introduce the locally adaptive step sizes also within standard HMC and Metropolis-corrected kinetic Langevin/generalized HMC (GHMC). Finally, a numerical study of the relative merits of the proposed error measures in the context of different Hamiltonian flow-based MCMC methods is provided.
Split Gibbs sampling (SGS) is a versatile tool to sample from complex posterior distributions, such as those arising in imaging inverse problems. However, existing SGS instances are limited to quadratic variable couplings, restricting their applicability to Gaussian likelihoods. This letter genearalizes SGS to non-Euclidean geometry by introducing a coupling kernel defined through a Jensen-Bregman divergence. This framework not only enables posterior sampling in these settings but also unifies several \emph{ad hoc} sampling schemes from the literature under a principled theoretical foundation.
We consider the problem of estimating the expectation of a function (observable) with respect to a high-dimensional probability distribution, under the assumption that a low-dimensional collective variable map captures all the directions of slow mixing for the distribution. More precisely, under the assumption that the fibers of the collective variable map satisfy a log-Sobolev inequality, we establish non-asymptotic relative-error guarantees for this problem using the eigenvector method for umbrella sampling (EMUS). Our analysis applies to both log-Lipschitz observables and indicator of convex sets, and gives complexity bounds that are polynomial in regularity parameters and inverse relative accuracy, and exponential in the dimension of the collective variable space. As a relative error guarantee, our bounds apply to rare event estimation. We consider Gaussian biasing windows defined on a low-dimensional collective variable space and analyze a modified sparse matrix version of EMUS. The bounds quantify the effect of the map through mixing of the conditional distribution and fluctuation of the observable along the fibers.
We construct an approximate sampler for finite doubly free-boundary Macdonald rail-yard processes with Macdonald, boundary, and edge parameters in $(0,1)$. The construction allows arbitrary finite sequences of additions and removals of horizontal or vertical strips and all four boundary types. It reverses a finite commutation--reflection history from a partition sequence whose entries are equal, with the common partition weighted by both boundaries. We derive local reconstruction kernels from the skew Cauchy, Littlewood, and branching identities. Bulk and boundary updates combine integer increments with conditional shape sampling at a prescribed size; same-sign swaps use finite conditional resampling. For $K$ reflection periods and an initial partition-size cutoff $M$, the algorithm terminates almost surely and samples the resulting approximating law exactly. We prove total variation convergence to the target process as $K,M\to\infty$ and give explicit bounds for both truncation errors, decaying geometrically in the cutoffs for fixed input parameters. Under fixed parameter bounds, we also bound the expected operation count of direct enumeration in a unit-cost model. We further construct an exact rejection sampler for the initial partition law without a size cutoff or evaluation of its normalizing constant. Using this sampler removes the size-truncation error, leaving only the error from the finite reflection history.
Sampling Gaussian measures restricted to nonlinear manifolds in high dimensions requires a proposal kernel that respects both geometric constraints and dimensional scaling requirements. We introduce tangent preconditioned Crank-Nicolson (pCN), which combines a manifold projection with a cancellation of the ambient Gaussian contribution, leaving only a geometric residual. Considering a shifted sphere with $\sqrt d$ scaling of radius and centre as a tractable benchmark embedded manifold, we prove that uncorrected tangent pCN results in a proposal noise scale of order $d^{-1/6}$ (compared to $d^{-1/2}$ for a tangent random-walk). We establish asymptotic optimal tuning for the geodesic expected squared jump distance which is maximised at an acceptance probability of 0.574 for the family of drift-corrected tangent proposals containing tangent pCN and tangent MALA. We also prove a joint diffusion limit for the parallel and orthogonal coordinates to the sphere's centre vector of tangent pCN. The orthogonal diffusion speed matches the limiting geodesic expected squared jump distance and is maximised at acceptance probability 0.574, whereas the parallel speed is maximised at 0.712, demonstrating that optimal tuning depends on the observable. We then identify an intrinsic precision correction characterised by potential Hessian and Ricci curvature that yields a proposal kernel with dimension-independent fixed step size limit. We discuss a generalisation to arbitrary hypersurfaces and present nondegenerate scaling of tangent pCN applied to conditioned diffusions under resolution refinement.
In modern sequential Monte Carlo (SMC) applications, such as diffusion model samplers, propagating a particle is expensive, so only a few particles can be afforded. In this regime, each resampling step decides which regions of the state space survive, so a better resampling scheme is worth extra computation. However, classical resampling methods only replicate existing particles and suffer from sample impoverishment, especially when the particle size is small relative to the dimension. To create new particles, we relax the unbiasedness of resampling and control the mean squared error via a Wasserstein criterion. Specifically, we propose Wasserstein-Dirichlet resampling (WDR), which constructs an equal-weight measure that optimally approximates the weighted particle distribution by solving a free-support Wasserstein barycenter problem. We then randomize the resulting optimal coupling by Dirichlet mixing. Theoretically, we provide non-asymptotic mean squared error bounds that quantify the bias-variance trade-off and show that the rate depends on the intrinsic dimension. A truncated variant preserves high-weight particles while applying WDR to the remaining mass, targeting mild weight degeneracy. Experiments on state-space models and conditional image generation demonstrate that WDR consistently improves estimation accuracy, stability, and sample diversity. The improvements persist at matched wall-clock cost when the propagation cost dominates.
The Hybrid forecast covariance model is an operationally successful way to build a background error covariance matrix. It is constructed as a weighted sum, with scalar weights, of the climatological covariance matrix and a localized ensemble covariance matrix, with a localization matrix typically depending on a localization length. For best performance, the Hybrid needs to be optimized by tuning those weights and, possibly, other parameters, such as the localization length. Here, a systematic and automated method for this tuning operation is introduced. This technique is referred to as hypervariational method, and it is formulated in terms of a cost function, built using a time series of innovations as an essential ingredient. This cost function, which depends on the parameters, is minimized when they have optimal values. Equations for its minimization, based on Newton's method, are derived. Besides the aforementioned, commonly used Hybrid, the hypervariational method is applicable to any known Hybrid (and to any set of parameters on which it depends analytically), including recently developed Hybrids, not necessarily constructed as weighted sums. Here, the hypervariational method is applied to tune the ensemble covariance matrix weight in the seminal-formulation Hybrid (which does not feature localization). This is done in a synthetic experimental setting where a theoretically-derived expression for the optimal value for that weight is known. Moreover, in that experimental setting, a further method, based on a linear fit procedure, to determine that optimal value is possible. Both the hypervariational method and that linear fit procedure are shown to produce an optimal weight value matching (up to a tolerance of $0.1$) the optimal value known from the theory. Furthermore, Newton's minimization in the hypervariational method is shown to converge immediately (at the first recursive step).
This article develops optimal order (up to logarithmic factors) weighted cubature rules for Monte Carlo numerical integration in Sobolev reproducing kernel Hilbert spaces over common matrix manifolds. Specifically, we consider matrix models of the Stiefel, Grassmann, flag, and Lagrangian Grassmann manifolds each over either real numbers $\mathbb{R}$, complex numbers $\mathbb{C}$, or quaternions $\mathbb{H}$. The proposed cubature rules use evaluation points which are independent and identically distributed uniform on the manifold. The weights are chosen to be optimal with respect to the reproducing kernel Hilbert space worst-case error, under the constraint that the weights sum to one. Crucially, this constraint ensures optimal weights may be found without having to evaluate the kernel mean embedding (the integral of the kernel), which is often difficult to compute. We consider computable Matérn kernels of Frobenius distances between matrices which reproduce the Sobolev reproducing kernel Hilbert spaces with equivalent norms. Open-source numerical experiments validate our theoretical development, including an application to simultaneous approximation of intrinsic volumes in the Cauchy--Kubota formula. We also present preliminary empirical results beyond independent points including transformations of low-discrepancy quasi-Monte Carlo points from the unit cube and jointly optimized points and weights.
Existing cumulative regret guarantees for Bayesian optimization (BO) do not necessarily characterize the quality of the algorithm's final output. This gap is critical under noisy evaluations, where the goal is to return a single high-quality recommendation. To address it, we introduce efficient exploration strategies that can be incorporated into any BO procedure by occasionally replacing the acquisition query with an exploration query. To validate this approach, we bound the Gaussian process posterior standard deviation over the search space via the fill distance of the queries, which measures the largest region left unexplored. We then use this result to derive upper bounds on the inference regret, i.e., the gap between the global optimum and the function value at the algorithm's returned point. Empirically, these exploration strategies improve the final recommendations of existing BO algorithms on synthetic benchmarks, hyperparameter optimization, and a control task.
Verifier-guided inference-time scaling uses additional computation to score intermediate states and steer generation toward better outputs. However, imperfect intermediate verifiers can distort sampling, and repeated refinement alone does not ensure recovery of the desired reward-tilted distribution. We introduce Continuous-State Verifier-Guided Backtracking (CS-VGB), an inference-time scaling method for sampling from reward-tilted distributions with continuous generative models through repeated denoising and renoising. With correctly paired denoising and renoising transitions, the stationary distribution of CS-VGB at the clean checkpoint equals the desired reward-tilted target, even when intermediate verifiers are inaccurate. Under stated regularity conditions, CS-VGB converges exponentially to the target despite verifier error, and the expected number of updates needed to reach a clean sample grows only linearly with the number of checkpoints. In analytically tractable settings, we empirically demonstrate that CS-VGB recovers the desired reward-tilted distribution despite verifier error, whereas the compared verifier-based sampling methods fail to recover the target under the same verifier errors. We also instantiate the framework with diffusion models and one-step stochastic posterior maps. Across prompt-guided ImageNet generation, facial-attribute generation on CelebA, and multi-robot planning, recurrent CS-VGB search improves reward and task success as inference budgets grow, with image-generation rewards continuing beyond the observed plateaus of denoising-only baselines.
Generalized neural network mixed models (GNMMs) combine neural networks and mixed-effects models to analyze longitudinal and clustered data. Although GNMMs have shown good predictive performance in practice, their statistical properties are not well understood. Therefore, this paper studies several important properties of GNMM and provides theoretical support for its use. First, we show that the neural-network component can approximate an unknown nonlinear fixed-effect function arbitrarily well and establish consistency of the fitted neural- network function for its population target within a given network class. Next, we establish consistency of random-effect estimation under increasing within-subject information and show that the GNMM predictor converges to the true predictor. We then prove consistency of the conditional mean function and show that the excess mean squared prediction error (MSPE) relative to the oracle prediction error converges to zero. We conduct simulation studies to examine nonlinear fixed-effect recovery, random-effect estimation, predictor and conditional mean estimation, and excess prediction error under different numbers of subjects and repeated observations. We also apply GNMM to a longitudinal schizophrenia dataset to predict continuous and binary depression outcomes using smartphone- based behavioral measurements and compare its predictive performance with generalized linear mixed models and conventional artificial neural networks. Both the theoretical and empirical results demonstrate that GNMM provides a flexible framework for modeling nonlinear longitudinal and clustered data. These findings provide a statistical foundation for GNMM prediction and support future methodological developments.
Clustering is typically performed using only object features, although additional information about relationships between objects may be available and can be incorporated through constraints. Among these, pairwise constraints require two objects to belong to the same cluster or to different clusters. Several extensions of the k-means algorithm accommodate such constraints, but they sometimes fail to find a feasible assignment or incur high computational cost. Motivated by these limitations, we propose an algebraic approach to k-means clustering with pairwise constraints. The algorithm simplifies the constraints using Gröbner bases before optimization, improving the efficiency of cluster assignment. It can also detect whether the constraint set has a feasible assignment before optimization. Comparative experiments with existing algorithms showed that the proposed algorithm achieved comparable or better clustering quality at a lower computational cost, particularly when the numbers of constraints and clusters were small.
Human reference labels used to evaluate medical AI can vary across readers and annotation workflows. We quantified agreement of a fixed pathology AI system with a future qualified reader without treating inferred consensus as error-free biological truth. Linked Evaluation and Unbiased NuCLS multirater data were analyzed using a fixed external HoVer-Net decision vector. The selected frame contained 1,144 anchors, 52 fields of view, five slides/patients, and five primary readers. The target was fixed-AI agreement with a future reader, with descriptive contrast $Δ=θ_E-θ_U$. Analyses included raw agreement and iMRMC benchmarks, reader-only and shared-detectability Bayesian models, structural sensitivities, adversarial simulations, and predictive assessment. Among 5,188 reader-anchor pairs observed in both conditions, 19.62% changed AI-agreement state, while raw agreement was 0.5466 in Unbiased and 0.5473 in Evaluation. The reader-only model gave $Δ=-0.0008$ (95% CrI, $-0.0033$ to $0.0008$), and the shared-detectability model gave $Δ\approx0$ ($-0.0140$ to $0.0127$). Across structural perturbations, absolute agreement ranged from 0.420 to 0.646 in Unbiased and 0.409 to 0.640 in Evaluation. Frozen adversarial simulations met the comparative-stability criterion (median $R=0.295$; 14/16 model-scenario cells had $R<1$). Held-out-AI prediction supported shared predictive information, but varied by slide composition. The Evaluation-minus-Unbiased contrast was comparatively stable, whereas absolute agreement was structurally sensitive. These findings do not establish equivalence, causal assistance effects, biological accuracy, or new-patient transportability. Evaluations using variable human references should report reader population, selected frame, observation mechanism, and structural sensitivity.
The Copula Active Subspaces (CAS) method of Part I estimates a non-Gaussian observation-noise density $π_{\boldsymbol{n}}$ from residual samples. Its KL divergence $\mathcal{D}_{\mathrm{KL}}(π_{\boldsymbol{n}}\,\|\,\widehatπ_{\boldsymbol{n}})$ admits an exact decomposition into marginal, projection and reduced-density terms and bounds the posterior KL divergence averaged over the true predictive distribution. We derive a priori bounds for these terms, together with estimates or bounds from residual samples. Weighted nearest-neighbour entropy estimates are used for the marginal and reduced-density terms, with asymptotic standard errors in exact rank coordinates. For the projection error, the score is fitted on one subsample, and its mean squared magnitude outside the retained subspace is estimated on the other. Corrections for score and rank-transform errors then give a theoretical upper bound with probability at least $1-δ$. Tests on three known noise laws in dimension $d=20$ use the estimated rank transform. In $11$ of $12$ combinations of example and sample size, the sum of the two estimated terms and the trace bound exceeds the reference divergence after averaging over the estimated models. In these tests, score-estimation error is the largest contribution to looseness. The overestimate from the Gaussian log-Sobolev inequality alone is at most a factor of $1.98$.
The cvxoptdes R-package provides a unified framework to perform convex optimal experiment design and online response optimization for linear, generalized linear and additive, nonlinear, and random intercept mixed models, using exclusively open-source solvers. It supports optimization and augmentation of weighted, penalized and convex constrained continuous designs across a discrete grid of candidate points and rounding to integer-valued designs using point-exchange algorithms or mixed-integer programming for orthogonally constrained designs, including the class of orthogonal minimally aliased response surface (OMARS) designs. Models and design constraints are defined and updated in an object-oriented context with the aim to facilitate the process of sequential (penalized) design augmentation analogous to online design procedures known from active learning and Bayesian optimization.