Papers

Gradient-variation online learning has drawn increasing attention due to its deep connections to game theory and optimization. It has been studied extensively in the full-information setting, but is underexplored with bandit feedback. In this work, we focus on gradient variation in Bandit Convex Optimization (BCO) with two-point feedback. By proposing a refined analysis of the non-consecutive gradient variation, a fundamental quantity in gradient variation with bandit feedback, we improve the dimension dependence for both convex and strongly convex functions compared with the best known results (Chiang et al., 2013). Our improved analysis of the non-consecutive gradient variation also implies other favorable problem-dependent guarantees, such as gradient-variance and small-loss regret bounds. Beyond the two-point setup, we demonstrate the versatility of our technique by achieving the first gradient-variation bound for one-point bandit linear optimization over hyper-rectangular domains. Finally, we validate the effectiveness of our results in more challenging tasks such as dynamic and universal regret minimization, establishing the first gradient-variation dynamic and universal regret bounds for two-point BCO.

Missing data are ubiquitous in empirical databases, yet statistical analyses typically require complete data matrices. Multiple imputation offers a principled solution for filling these gaps. This study evaluates the performance of several multiple imputation methods, both in the presence and absence of extreme values, using the MICE package in R. Through Monte Carlo simulations, we generated incomplete data sets with three variables and assessed each imputation method within regression models. The results indicate that the linear regression based imputation method showed the best overall predictive performance (CV-MSE), whereas the sparse model approach was generally less efficient. Our findings underscore the relevance of extreme values when selecting an imputation strategy and highlight sample size, proportion of missingness, presence of extremes, and the type of fitted model as key determinants of performance. Despite its limitations, the study offers practical recommendations for researchers, stressing the need to examine the missingness mechanism and the occurrence of extreme values before choosing an imputation method.

A complete understanding of heterogeneous treatment effects involves characterizing the full conditional distribution of potential outcomes. To this end, we propose the Conditional Counterfactual Mean Embeddings (CCME), a framework that embeds conditional distributions of counterfactual outcomes into a reproducing kernel Hilbert space (RKHS). Under this framework, we develop a two-stage meta-estimator for CCME that accommodates any RKHS-valued regression in each stage. Based on this meta-estimator, we develop three practical CCME estimators: (1) Ridge Regression estimator, (2) Deep Feature estimator that parameterizes the feature map by a neural network, and (3) Neural-Kernel estimator that performs RKHS-valued regression, with the coefficients parameterized by a neural network. We provide finite-sample convergence rates for all estimators, establishing that they possess the double robustness property. Our experiments demonstrate that our estimators accurately recover distributional features including multimodal structure of conditional counterfactual distributions.

The wave speed of a stochastic wave equation driven by Riesz noise on the unbounded multidimensional spatial domain is estimated based on discrete measurements. Central limit theorems for second-order variations of the observations in space, time, and space-time are established. Under general assumptions on the spatial and temporal sampling frequencies, the resulting method-of-moments estimators are asymptotically normally distributed. The covariance structure of the discrete increments admits a closed-form representation involving two different Fejér-type kernels, enabling a precise analysis of the interplay between spatial and temporal contributions.

This article proposes a novel estimator for regression coefficients in clustered data that explicitly accounts for within-cluster dependence. We study the asymptotic properties of the proposed estimator under both finite and infinite cluster sizes. The analysis is then extended to a standard random coefficient model, where we derive asymptotic results for the average (common) parameters and develop a Wald-type test for general linear hypotheses. We also investigate the performance of the conventional pooled ordinary least squares (POLS) estimator within the random coefficients framework and show that it can be unreliable across a wide range of empirically relevant settings. Furthermore, we introduce a new test for parameter stability at a higher (superblock; Tier 2, Tier 3,...) level, assuming that parameters are stable across clusters within that level. Extensive simulation studies demonstrate the effectiveness of the proposed tests, and an empirical application illustrates their practical relevance.

Network data are increasingly common in the social sciences and infectious disease epidemiology. Analyses often link network structure to node-level covariates, but existing methods falter with sparse networks and high-dimensional node features. We propose a joint latent space modeling framework for sparse networks with high-dimensional binary node covariates that performs covariate selection while accounting for uncertainty in estimated latent positions. Building on joint latent space models that couple edges and node variables through shared latent positions, we introduce a group lasso screening step and incorporate a measurement-error-aware stabilization term to mitigate bias from using estimated latent positions as predictors. We establish prediction error rates for the covariate component both when latent positions are treated as observed and when they are estimated with bounded error; under uniform control across $q$ covariates and $n$ nodes, the rate is of order $O(\log q / n)$ up to an additional term due to latent position estimation error. Our method addresses three challenges: (1) incorporating information from isolated nodes, which are common in sparse networks but often ignored; (2) selecting relevant covariates from high-dimensional spaces; and (3) accounting for uncertainty in estimated latent positions. Simulations show predictive performance remains stable as covariate sparsity grows, while naive approaches degrade. We illustrate how the method can support efficient study design using household social networks from 75 Indian villages, where an emulated pilot study screens a large covariate battery and substantially reduces required subsequent data collection without sacrificing network predictive accuracy.

The geography of innovation offers a framework to understand how territorial characteristics shape innovation, often via spatial and cognitive proximity. Empirical research has focused largely on national and regional scales, while urban and sub-regional geographies receive less attention. Local studies typically rely on limited indicators (e.g., firm-level data, patents, basic socioeconomic measures), with few offering a systematic framework integrating urban form, mobility, amenities, and human-capital proxies at the neighborhood scale. Our study investigates innovation at a finer spatial resolution, going beyond proprietary or static indicators. We develop the Local Innovation Determinants (LID) database and framework to identify key enabling factors across regions, combining traditional government data with publicly available data via APIs for a more granular understanding of spatial dynamics shaping innovation capacity. Using exploratory big and geospatial data analytics and random forest models, we examine neighborhoods in New York and Massachusetts across four dimensions: social factors, economic characteristics, land use and mobility, morphology, and environment. Results show that alternative data sources offer significant yet underexplored potential to enhance insights into innovation dynamics. City policymakers should consider neighborhood-specific determinants and characteristics when designing and implementing local innovation strategies.

Stochastic processes are often represented through orthonormal series expansions, a framework originating in the classical works of Loève and Karhunen and widely used for simulation and numerical approximation. While truncation error in such expansions has been extensively studied, practical models frequently involve an additional source of error arising from the approximation of coefficient functions when closed-form expressions are unavailable. The combined effect of these two errors remains insufficiently addressed in the literature. Building on the author's earlier work on reliability and accuracy estimates for $\varphi$-sub-Gaussian processes, this paper extends the methodology to orthonormal polynomial systems that do not possess normalized generating functions in analytical form, including the Legendre, generalized Laguerre, and Gegenbauer families. New bounds are derived for models in $L_p(T)$ space that simultaneously account for truncation and coefficient approximation. The resulting criteria provide practical guidance for selecting the number of series terms required to achieve prescribed levels of reliability and accuracy across a broader class of polynomial-based stochastic process models.

Root-cause analysis in controlled time dependent systems poses a major challenge in applications. Especially energy systems are difficult to handle as they exhibit instantaneous as well as delayed effects and if equipped with storage, do have a memory. In this paper we adapt the causal root-cause analysis method of Budhathoki et al. [2022] to general time-dependent systems, as it can be regarded as a strictly causal definition of the term "root-cause". Particularly, we discuss two truncation approaches to handle the infinite dependency graphs present in time-dependent systems. While one leaves the causal mechanisms intact, the other approximates the mechanisms at the start nodes. The effectiveness of the different approaches is benchmarked using a challenging data generation process inspired by a problem in factory energy management: the avoidance of peaks in the power consumption. We show that given enough lags our extension is able to localize the root-causes in the feature and time domain. Further the effect of mechanism approximation is discussed.

HIV transmission within serodiscordant couples remains a significant public health challenge, particularly in sub-Saharan Africa. Estimating the rate of such infection, alongside the rates of introduction of infection from outside the partnership, is a special case of the more general epidemiological challenge of inferring intensities of within- and between-group intensities of transmission. This study presents a stochastic susceptible-infected (SI) pair model for estimating key epidemiological parameters governing HIV transmission within and between couples, which we further extend to account for gender-specific differences in infection dynamics. Using a likelihood-based inference approach, we estimate transmission parameters and associated uncertainty from observed data. These values can be used to inform infection prevention strategies for HIV, and the methodology proposed can be generalised to other epidemiological settings.

The synthetic control method (SCM) estimates causal effects in panel data with a single-treated unit by constructing a counterfactual outcome as a weighted combination of untreated control units that matches the pre-treatment trajectory. In this paper, we introduce the targeted synthetic control (TSC) method, a new two-stage estimator that directly estimates the counterfactual outcome. Specifically, our TSC method (1) yields a targeted debiasing estimator, in the sense that the targeted updating refines the initial weights to produce more stable weights; and (2) ensures that the final counterfactual estimation is a convex combination of observed control outcomes to enable direct interpretation of the synthetic control weights. TSC is flexible and can be instantiated with arbitrary machine learning models. Methodologically, TSC starts from an initial set of synthetic-control weights via a one-dimensional targeted update through the weight-tilting submodel, which calibrates the weights to reduce bias of weights estimation arising from pre-treatment fit. Furthermore, TSC avoids key shortcomings of existing methods (e.g., the augmented SCM), which can produce unbounded counterfactual estimates. Across extensive synthetic and real-world experiments, TSC consistently improves estimation accuracy over state-of-the-art SCM baselines.

Bayes-filtered transformers are transformers meta-learned on sequences from a prior predictive distribution to approximate the corresponding posterior predictive distribution. They output total predictive uncertainty in a single forward pass but never explicitly represent a posterior distribution, making the standard route to separating aleatoric from epistemic uncertainty unavailable. We address this challenge through the lens of Bayesian predictive inference (BPI). Our main result is a predictive Central Limit Theorem (CLT) for supervised settings under conditions that are among the weakest known in the BPI literature. The CLT characterises the posterior of the limiting predictive distribution given an observed context as asymptotically Gaussian; the variance of this Gaussian quantifies epistemic uncertainty. We apply the framework to TabPFN, a Bayes-filtered transformer that is a state-of-the-art foundation model for tabular prediction. The resulting credible bands achieve near-nominal frequentist coverage as context length grows, and the decomposition largely matches standard desiderata: epistemic uncertainty shrinks with context length and is highest in sparsely observed regions within the span of the context data, while aleatoric uncertainty dominates near decision boundaries where classes overlap.

We study high-dimensional rank regression when data are distributed across multiple machines and the loss is a non-additive U-statistic, as in convoluted rank regression (CRR). Classical communication-efficient surrogate likelihood (CSL) methods crucially rely on the additivity of the empirical loss and therefore break down for CRR, whose global loss couples all sample pairs across machines. We propose a distributed convoluted rank regression (DCRR) framework that constructs a similar surrogate loss and demonstrate its validity under the non-additive losses. We show that this surrogate shares the same population minimizer as the full-data CRR loss and yields estimators that are statistically equivalent to centralized CRR. Building on this, we develop a two-stage sparse DCRR procedure – an iterative $\ell_1$-penalized stage followed by a folded-concave refinement – and establish non-asymptotic error bounds, a distributed strong oracle property, and a DHBIC-type criterion for consistent model selection. A scaling result shows that the number of machines may diverge as $M = o({N/(s^2\log p)})$ while achieving centralized oracle rates with only $O(\log N)$ communication rounds. Simulations and a large-scale real data example demonstrate substantial gains over naive divide-and-conquer, particularly under heavy-tailed errors.

Identifying differentially methylated regions is an important task in epigenome-wide association studies, where differential signals often arise across groups of neighboring CpG sites. Many existing methods detect differentially methylated regions by aggregating CpG-level test results, which may limit their ability to capture complex regional methylation patterns. In this paper, we introduce the R package mmcmcBayes, which implements a multistage Markov chain Monte Carlo procedure for region-level detection of differentially methylated regions. The method models sample-wise regional methylation summaries using the alpha-skew generalized normal distribution and evaluates evidence for differential methylation between groups through Bayes factors. We use a multistage region-splitting strategy to refine candidate regions based on statistical evidence. We describe the underlying methodology and software implementation, and illustrate its performance through simulation studies and applications to Illumina 450K methylation data. The mmcmcBayes package provides a practical region-level alternative to existing CpG-based differentially methylated regions detection methods and includes supporting functions for summarizing, comparing, and visualizing detected regions.

Standard approaches to quantum statistical inference rely on measurements that induce a collapse of the wave function, effectively consuming the quantum state to extract information. In this work, we investigate the fundamental limits of locally-gentle quantum state certification, where the learning algorithm is constrained to perturb the state by at most $\alpha$ in trace norm, thereby allowing for the reuse of samples. We analyze the hypothesis testing problem of distinguishing whether an unknown state $\rho$ is equal to a reference $\rho_0$ or $\epsilon$-far from it. We derive the minimax sample complexity for this problem, quantifying the information-theoretic price of non-destructive measurements. Specifically, by constructing explicit measurement operators, we show that the constraint of $\alpha$-gentleness imposes a sample size penalty of $\frac{d}{\alpha^2}$, yielding a total sample complexity of $n = \Theta(\frac{d^3}{\epsilon^2 \alpha^2})$. Our results clarify the trade-off between information extraction and state disturbance, and highlight deep connections between physical measurement constraints and privacy mechanisms in quantum learning. Crucially, we find that the sample size penalty incurred by enforcing $\alpha$-gentleness scales linearly with the Hilbert-space dimension $d$ rather than the number of parameters $d^2-1$ typical for high-dimensional private estimation.

We develop a general mathematical framework to analyze scaling regimes and derive explicit analytic solutions for gradient flow (GF) in large learning problems. Our key innovation is a formal power series expansion of the loss evolution, with coefficients encoded by diagrams akin to Feynman diagrams. We show that this expansion has a well-defined large-size limit that can be used to reveal different learning phases and, in some cases, to obtain explicit solutions of the nonlinear GF. We focus on learning Canonical Polyadic (CP) decompositions of high-order tensors, and show that this model has several distinct extreme lazy and rich GF regimes such as free evolution, NTK and under- and over-parameterized mean-field. We show that these regimes depend on the parameter scaling, tensor order, and symmetry of the model in a specific and subtle way. Moreover, we propose a general approach to summing the formal loss expansion by reducing it to a PDE; in a wide range of scenarios, it turns out to be first-order and solvable by the method of characteristics. We observe a very good agreement of our theoretical predictions with experimental results.

In the context of election security, a Risk-Limiting Audit (RLA) is a statistical framework that uses a minimal partial recount of the ballots to guarantee that the results of the election were correctly reported. A generalized RLA framework has remained elusive for algorithmic election rules such as the Single Transferable Vote (STV) rule, because of the dependence of these rules on the chronology of eliminations and elections leading to the outcome of the election. This paper proposes a new graph-based approach to audit these algorithmic election rules, by considering the space of all possible sequences of elections and eliminations. If we fix a subgraph of this universal space ahead of the audit, a sufficient strategy is to verify statistically that the true election sequence does not leave the fixed subgraph. This makes for a flexible framework to audit these elections in a chronology-agnostic way.

We study asymmetric rank-one spiked tensor models in the high-dimensional regime, where the noise entries are independent and identically distributed with zero mean, unit variance, and finite fourth moment. This extends the classical Gaussian framework to a substantially broader class of noise distributions. We analyze the maximum-likelihood estimator associated with the best rank-one approximation of an order-$d$ tensor, for $d\ge 3$. Our approach is formulated along an informative, spectrally separated branch of stationary points of the non-convex maximum-likelihood landscape. In the core order-three asymmetric model, we verify locally in the high-signal regime that such an informative branch exists and remains separated from the bulk. Under this branch-selection framework, we show that the empirical spectral distribution of a suitable block-wise tensor contraction converges almost surely to the same deterministic limit as in the Gaussian case. As a consequence, the asymptotic singular value and the mode-wise alignments between the estimated and planted spike directions admit the same explicit characterizations as under Gaussian noise. These results establish a universality principle for asymmetric spiked tensor models: the high-dimensional spectral behavior and statistical limits of the selected maximum-likelihood stationary point are robust beyond the Gaussian setting. Our proof combines resolvent methods from random matrix theory, cumulant expansions under finite fourth-moment assumptions, and Efron–Stein-type variance bounds. A main technical difficulty is to control the statistical dependence between the estimator and the noise, including the associated cross terms in the non-Gaussian setting.

The main contribution of this paper is to develop a hierarchical Bayesian formulation of PINNs for linear inverse problems, which is called BPINN-IP. The proposed methodology extends PINN to account for prior knowledge on the nature of the expected NN output, as well as its weights. Also, as we can have access to the posterior probability distributions, naturally uncertainties can be quantified. Also, variational inference and Monte Carlo dropout are employed to provide predictive means and variances for reconstructed images. Un example of applications to deconvolution and super-resolution is considered, details of the different steps of implementations are given, and some preliminary results are presented.

A/B testing on platforms often faces challenges from network interference, where a unit's outcome depends not only on its own treatment but also on the treatments of its network neighbors. To address this, cluster-level randomization has become standard, enabling the use of network-aware estimators. These estimators typically trim the data to retain only a subset of informative units, achieving low bias under suitable conditions but often suffering from high variance. In this paper, we first demonstrate that the interior nodes - units whose neighbors all lie within the same cluster - constitute the vast majority of the post-trimming subpopulation. In light of this, we propose directly averaging over the interior nodes to construct the mean-in-interior (MII) estimator, which circumvents the delicate reweighting required by existing network-aware estimators and substantially reduces variance in classical settings. However, we show that interior nodes are often not representative of the full population, particularly in terms of network-dependent covariates, leading to notable bias. We then augment the MII estimator with a counterfactual predictor trained on the entire network, allowing us to adjust for covariate distribution shifts between the interior nodes and full population. By rearranging the expression, we reveal that our augmented MII estimator embodies an analytical form of the point estimator within prediction-powered inference framework. This insight motivates a semi-supervised lens, wherein interior nodes are treated as labeled data subject to selection bias. Extensive and challenging simulation studies demonstrate the outstanding performance of our augmented MII estimator across various settings.

We study the Pareto frontier (optimal trade-off) between utility and separation, a fairness criterion requiring predictive independence from sensitive attributes conditional on the true outcome. Through an information-theoretic lens, we prove a characterization of the utility-separation Pareto frontier, establish its concavity, and thereby prove the increasing marginal cost of separation in terms of utility. In addition, we characterize the conditions under which this trade-off becomes strict, providing a guide for trade-off selection in practice. Based on the theoretical characterization, we develop an empirical regularizer based on conditional mutual information (CMI) between predictions and sensitive attributes given the true outcome. The CMI regularizer is compatible with any deep model trained via gradient-based optimization and serves as a scalar monitor of residual separation violations, offering tractable guarantees during training. Finally, numerical experiments support our theoretical findings: across COMPAS, UCI Adult, UCI Bank, and CelebA, the proposed method substantially reduces separation violations while matching or exceeding the utility of established baseline methods. This study thus offers a provable, stable, and flexible approach to enforcing separation in deep learning.

Performative predictions influence the very outcomes they aim to forecast. We study performative predictions that affect a sample (e.g., only existing users of an app) and/or the whole population (e.g., all potential app users). This raises the question of how well models generalize under performativity. For example, how well can we draw insights about new app users based on existing users when both of them react to the app's predictions? We address this question by embedding performative predictions into statistical learning theory. We prove generalization bounds under performative effects on the sample, on the population, and on both. A key intuition behind our proofs is that in the worst case, the population negates predictions, while the sample deceptively fulfills them. We cast such self-negating and self-fulfilling predictions as min-max and min-min risk functionals in Wasserstein space, respectively. Our analysis reveals a fundamental trade-off between performatively changing the world and learning from it: the more a model affects data, the less it can learn from it. Moreover, our analysis results in a surprising insight on how to improve generalization guarantees by retraining on performatively distorted samples. We illustrate our bounds in a case study on prediction-informed assignments of unemployed German residents to job trainings, drawing upon administrative labor market records from 1975 to 2017 in Germany.

There is a growing need for flexible statistical distributions that can accurately model data defined on the unit interval. This paper introduces a new unit distribution, termed the unit Shiha (USh) distribution, which is derived from the original Shiha (Sh) distribution through an inverse exponential transformation. The probability density function of the USh distribution is sufficiently flexible to model both left- and right-skewed data, while its hazard rate function is capable of capturing various failure-rate patterns, including increasing, bathtub-shaped, and J-shaped forms. Several statistical properties of the proposed distribution are investigated, including moments and related measures, the quantile function, entropy, and stress-strength reliability. Parameter estimation is carried out using the maximum likelihood method, and its performance is evaluated through a simulation study. The practical usefulness of the USh distribution is demonstrated using four real-life data sets, and its performance is compared with several well-known competing unit distributions. The comparative results indicate that the proposed model fits the data better than the competitive models applied in this study.

Prediction sets provide a means of quantifying the uncertainty in predictive tasks. Using held out calibration data, conformal prediction and risk control can produce prediction sets that exhibit statistically valid error control in a computationally efficient manner. However, in the standard formulations, the error is only controlled on average over many possible calibration datasets of fixed size. In this paper, we extend the control to remain valid with high probability over a cumulatively growing calibration dataset at any time point. We derive such guarantees using quantile-based arguments and illustrate the applicability of the proposed framework to settings involving distribution shift. We further establish a matching lower bound and show that our guarantees are asymptotically tight. Finally, we demonstrate the practical performance of our methods through both simulations and real-world numerical examples.

Suppose some cleverness score parameter is sufficiently interesting to be defined and then measured, perhaps for different strata of specialists or for the broader population. Such phenomena could have Gaussian distributions, when it comes to all players in a stratum, but when interest focuses on the very tails, for the top few percent, those above certain high thresholds, different models are called for, along with the need to analyse such based on the listed top scores only. In this note I develop such models and tools, and apply them to the top-100 and above 2100 points lists for regular chess ratings, for the currently active 14671 men and 753 women, as given by the FIDE, January 2026. It is argued that even when two or more distributions have close to identical expected values, or medians, even smaller differences in variance may explain gaps for the few very best ones.

In recent years, instructional practices in Operations Research (OR), Management Science (MS), and Analytics have increasingly shifted toward digital environments, where large and diverse groups of learners make it difficult to provide practice that adapts to individual needs. This paper introduces a method that generates personalized sequences of exercises by selecting, at each step, the exercise most likely to advance a learner's understanding of a targeted skill. The method uses information about the learner and their past performance to guide these choices, and learning progress is measured as the change in estimated skill level before and after each exercise. Using data from an online mathematics tutoring platform, we find that the approach recommends exercises associated with greater skill improvement and adapts effectively to differences across learners. From an instructional perspective, the framework enables personalized practice at scale, highlights exercises with consistently strong learning value, and helps instructors identify learners who may benefit from additional support.

Solving large scale Optimal Transport (OT) in machine learning typically relies on sampling measures to obtain a tractable discrete problem. While the discrete solver's accuracy is controllable, the rate of convergence of the discretization error is governed by the intrinsic dimension of our data. Therefore, the true bottleneck is the knowledge and control of the sampling error. In this work, we tackle this issue by introducing novel estimators for both sampling error and intrinsic dimension. The key finding is a simple, tuning-free estimator of $\text{OT}_c(\rho, \hat\rho)$ that utilizes the semi-dual OT functional and, remarkably, requires no OT solver. Furthermore, we derive a fast intrinsic dimension estimator from the multi-scale decay of our sampling error estimator. This framework unlocks significant computational and statistical advantages in practice, enabling us to (i) quantify the convergence rate of the discretization error, (ii) calibrate the entropic regularization of Sinkhorn divergences to the data's intrinsic geometry, and (iii) introduce a novel, intrinsic-dimension-based Richardson extrapolation estimator that strongly debiases Wasserstein distance estimation. Numerical experiments demonstrate that our geometry-aware pipeline effectively mitigates the discretization error bottleneck while maintaining computational efficiency.

We introduce smallest valid partitioning (SVP), a segmentation method for multiple change-point detection in time-series. SVP relies on a local notion of segment validity: a candidate segment is retained only if it passes a user-chosen validity test (e.g., a single change-point test). From the collection of valid segments, we propose a coherent aggregation procedure that constructs a global segmentation which is the exact solution of an optimization problem. Our main contribution is the use of a lexicographic order for the optimization problem that prioritizes parsimony. We analyze the computational complexity of the resulting procedure, which ranges from linear to cubic time depending on the chosen cost and validity functions, the data regime and the number of detected changes. Finally, we assess the quality of SVP through comparisons with standard optimal partitioning algorithms, showing that SVP yields competitive segmentations while explicitly enforcing segment validity. The flexibility of SVP makes it applicable to a broad class of problems; as an illustration, we demonstrate robust change-point detection by encoding robustness in the validity criterion.

Rao's spacing test is a widely used nonparametric method for assessing uniformity on the circle. However, its broader applicability in practical settings has been limited because the null distribution is not easily calculated. As a result, practitioners have traditionally depended on pre-tabulated critical values computed for a limited set of sample sizes, which restricts the flexibility and generality of the method. In this paper, we address this limitation by recursively computing higher-order moments of the Rao's spacing test statistic and employing the Gram-Charlier expansion to derive an accurate approximation to its null distribution. This approach allows for the efficient and direct computation of p-values for arbitrary sample sizes, thereby eliminating the dependency on existing critical value tables. Moreover, we confirm that our method remains accurate and effective even for large sample sizes that are not represented in current tables, thus overcoming a significant practical limitation. Comparative evaluations with published critical values and saddlepoint approximations demonstrate that our method achieves a high degree of accuracy across a wide range of sample sizes. These findings greatly improve the practicality and usability of Rao's spacing test in both theoretical investigations and applied statistical analyses.

We study the asymptotic spectral behavior of high-dimensional random Gram matrices with sparsity and a variance profile, motivated by applications in wireless communications. Specifically, we consider the Gram matrices $\mathbf S_n=\mathbf Y_n\mathbf Y_n^*$, where the entries of $\mathbf Y_n$ are independent, centered, heteroscedastic, and sparse through Bernoulli masking. The sparsity level is parameterized as $s=q^2/n$, where $q$ ranges from polynomial order up to order $n^{1/2}$. We investigate two asymptotic regimes: a moderate-sparsity regime with fixed $s\in(0,1]$, and a high-sparsity regime where $s\to0$. In both regimes, we establish the convergence of the empirical spectral distribution of $\mathbf S_n$ to a deterministic limit, and further derive central limit theorems for linear spectral statistics using resolvent techniques and martingale difference arguments. Our analysis reveals a phase transition in the fluctuation behavior across the two regimes. In the high-sparsity regime, the asymptotic fluctuations are entirely governed by fourth-moment effects, with sparsity-scaled contributions being suppressed. Moreover, the leading deterministic term and the variance of the linear spectral statistic scale at different rates in $q$, causing the standard centering to fail and necessitating an explicit correction to recover a valid CLT. The results apply to both Gaussian and non-Gaussian entries and are illustrated through applications to hypothesis testing and outage probability analysis in large-scale MIMO systems.