The Annals of Statistics

Distributed learning offers a practical solution for the integrative analysis of multisource datasets, especially under privacy or communication constraints. However, addressing prospective distributional heterogeneity and ensuring communication efficiency pose significant challenges on distributed statistical analysis. In this article we focus on integrative estimation of distributed heterogeneous precision matrices, a crucial task related to joint precision matrix estimation where computation-efficient algorithms and statistical optimality theories are still underdeveloped. To tackle these challenges, we introduce a novel HEterogeneity-adjusted Aggregating and Thresholding (HEAT) approach for distributed integrative estimation. HEAT is designed to be both communication- and computation-efficient, and we demonstrate its statistical optimality by establishing the convergence rates and the corresponding minimax lower bounds under various integrative losses. To enhance the optimality of HEAT, we further propose an iterative HEAT (IteHEAT) approach. By iteratively refining the higher-order errors of HEAT estimators through multi-round communications, IteHEAT achieves geometric contraction rates of convergence. Extensive simulations and real data applications validate the numerical performance of HEAT and IteHEAT methods.

We consider diffusion of independent molecules in an insulated Euclidean domain with unknown diffusivity parameter. At a random time and position, the molecules may bind and stop diffusing in dependence of a given “binding potential.” The binding process can be modeled by an additive random functional corresponding to the canonical construction of a “killed” diffusion Markov process. We study the problem of conducting inference on the infinite-dimensional diffusion parameter from a histogram plot of the “killing” positions of the process. We show first that these positions follow a Poisson point process whose intensity measure is determined by the solution of a Schrödinger-type equation. The inference problem can then be recast as a nonlinear inverse problem for this partial differential equation, which we show to be consistently solvable in a Bayesian way under natural conditions on the initial state of the system, provided the binding potential is not too “aggressive.” In the proofs, we obtain novel posterior contraction rate results for high-dimensional Poisson count data that are of independent interest. A numerical illustration of the estimator by MCMC methods is also provided.

We introduce a novel generative modeling framework based on a discretized parabolic Monge–Ampère PDE, which emerges as a continuous limit of the Sinkhorn algorithm commonly used in optimal transport. Our method performs iterative refinement in the space of Brenier maps using a mirror gradient descent step. We establish theoretical guarantees for generative modeling through the lens of no-regret analysis, demonstrating that the iterates converge to the optimal Brenier map under a variety of step-size schedules. As a technical contribution, we derive a new Evolution Variational Inequality tailored to the parabolic Monge–Ampère PDE, connecting geometry, transportation cost and regret. Our framework accommodates nonlog-concave target distributions, constructs an optimal sampling process via the Brenier map, and integrates favorable learning techniques from generative adversarial networks and score-based diffusion models. As direct applications, we illustrate how our theory paves new pathways for generative modeling and variational inference.

We study the estimation of the invariant density of additive fractional stochastic differential equations with Hurst parameter $\mathit{H}\in (0,1)$. We first focus on continuous observations and develop a kernel-based estimator achieving faster convergence rates than previously available. This result stems from a martingale decomposition combined with new bounds on the (conditional) convergence in total variation to equilibrium of fractional SDEs. For $\mathit{H}<1/2$, we further refine the rates based on recent bounds on the marginal density. We then extend the methodology to discrete observations, showing that the same convergence rates can be attained. Moreover, we establish concentration inequalities for the estimator and introduce a data-driven bandwidth selection procedure that adapts to unknown smoothness. Numerical experiments for the fractional Ornstein–Uhlenbeck process illustrate the estimator’s practical performance. Finally, our results weaken the usual convexity assumptions on the drift component, allowing us to consider settings where strong convexity only holds outside a compact set.

We provide universality results that quantify how data augmentation affects the variance and limiting distribution of estimates through simple surrogates, and analyze several specific models in detail. The results confirm some observations made in machine learning practice, but also lead to unexpected findings: Data augmentation may increase rather than decrease the uncertainty of estimates, such as the empirical prediction risk. It can act as a regularizer, but fails to do so in certain high-dimensional problems, and it may shift the double-descent peak of an empirical risk. Overall, the analysis shows that several properties attributed to data augmentation are not either true or false, but rather depend on a combination of factors—notably the data distribution, the properties of the estimator, and the interplay of sample size, number of augmentations and dimension. As our main theoretical tool, we develop an adaptation of Lindeberg’s technique for block dependence. The resulting universality regime may be Gaussian or non-Gaussian.

A functional Bernstein–von Mises theorem is proved for posterior measures arising in a data assimilation problem with the two-dimensional Navier–Stokes equation where a Gaussian process prior is assigned to the initial condition of the system. The posterior measure, which provides the update in the space of all trajectories arising from a discrete sample of the (deterministic) dynamics, is shown to be approximated by a Gaussian random vector field obtained as the solution to a linear parabolic PDE with Gaussian initial condition. The approximation holds in the strong sense of the supremum norm on the regression functions, showing that predicting future states of Navier–Stokes systems admits $1/\sqrt{\mathit{N}}$-consistent estimators even for commonly used nonparametric models. Consequences for coverage of credible bands and uncertainty quantification are discussed. A local asymptotic minimax theorem is derived that describes the lower bound for estimating the state of the nonlinear system, which is shown to be attained by the Bayesian data assimilation algorithm.

Given a random sample from a density function supported on a manifold M, a new method for the estimating highest density regions of the underlying population is introduced. The new proposal is based on the empirical version of the opening operator from mathematical morphology combined with a preliminary estimator of the density function. This results in an estimator that is easy-to-compute since it simply consists of a list of centers and a radius r that are adequately selected from the data. The new estimator is shown to be consistent and its convergence rates in terms of the Hausdorff distance are provided. All consistency results are established uniformly on the level of the set and for any Riemannian manifold M satisfying mild assumptions. The applicability of the procedure is shown by some illustrative examples.

The problem of estimating the slope parameter in regression between two spatial processes under confounding by an unmeasured spatial process has received widespread attention in the recent statistical literature. Yet, a fundamental question remains unresolved: when is this slope consistently estimable under spatial confounding, with existing insights being largely empirical or estimator-specific. In this manuscript, we characterize conditions for consistent estimability of the regression slope between Gaussian random fields (GRFs), the common stochastic model for spatial processes, under spatial confounding. Under fixed-domain (infill) asymptotics, we give sufficient conditions for consistent estimability in terms of the smoothness or local behavior of the exposure and confounder processes. When estimability holds, we provide consistent estimators of the slope using local differencing (taking discrete differences or Laplacians of the processes of suitable order). Using functional analysis results on Paley–Wiener spaces, we then provide an easy-to-verify necessary condition for consistent estimability of the slope in terms of the relative spectral tail decays of the confounder and exposure. As a by-product, we establish a novel and general spectral condition on the equivalence of measures on the paths of multivariate GRFs with component fields of varying smoothnesses, a result of independent importance. Our estimability results or estimators do not rely on specific parametric models for the covariance functions. We show that for many covariance classes like the Matérn, power-exponential, generalized Cauchy and coregionalization families, the necessary and sufficient conditions become identical, except for a boundary point, thereby providing a sharp characterization of consistent estimability of the slope for these processes. The results are extended to multivariate slopes and to accommodate measurement error using local-averaging-and-differencing-based estimators. We show that differencing-based estimation remains consistent for popular classes of nonstationary Gaussian random fields, some non-Gaussian random fields and irregular designs. The finite sample behavior of the estimators is explored via numerical experiments.

Bayesian model-based spatial clustering methods are widely used for their flexibility in estimating latent clusters with an unknown number of clusters while accounting for spatial proximity. Many existing methods are designed for clustering finite spatial units, limiting their ability to make predictions or may impose restrictive geometric constraints on the shapes of subregions. Furthermore, the posterior clustering consistency theory of spatial clustering models remains largely unexplored in the literature. In this study, we propose a Spatial Domain Random Partition Model (Spat-RPM) and demonstrate its application for spatially clustered regression, which extends spanning tree-based Bayesian spatial clustering by partitioning the spatial domain into disjoint blocks and using spanning tree cuts to induce contiguous domain partitions. Under an infill-domain asymptotic framework, we introduce a new distance metric to study the posterior concentration of domain partitions. We show that Spat-RPM achieves a consistent estimation of domain partitions, including the number of clusters (which may go to infinity), and derive posterior concentration rates for partition, parameter and prediction. We also establish conditions on the hyperparameters to achieve consistency, offering important practical guidance for hyperparameter selection. Finally, we examine the asymptotic properties of our model through simulation studies and apply it to Atlantic Ocean data.

We consider a general model for high-dimensional empirical risk minimization whereby the data ${\mathit{x}}_{\mathit{i}}$ are d-dimensional Gaussian vectors, the model is parametrized by $\mathrm{\Theta }\in {\mathbb{R}}^{\mathit{d}\times \mathit{k}}$ and the loss depends on the data via the projection ${\mathrm{\Theta }}^{\mathsf{T}}{\mathit{x}}_{\mathit{i}}$. This setting covers as special cases classical statistics methods (e.g., multinomial regression and other generalized linear models), but also two-layer fully connected neural networks with k hidden neurons.

We use the Kac–Rice formula from Gaussian process theory to derive a bound on the expected number of local minima of this empirical risk, under the proportional asymptotics in which $\mathit{n},\mathit{d}\to \infty$, with $\mathit{n}\asymp \mathit{d}$. Via Markov’s inequality, this bound allows one to determine the positions of these minimizers (with exponential deviation bounds), and hence derive sharp asymptotics on the estimation and prediction error.

As a special case, we apply our characterization to convex losses. We show that our approach is tight and allows us to prove previously conjectured results. In addition, we characterize the spectrum of the Hessian at the minimizer. A companion paper applies our general result to nonconvex examples.

We study the spectral properties of the Jacobian matrix in artificial neural networks and establish both its limiting spectral distribution and the second-order fluctuations of its spectrum. Our analysis uncovers fundamental phenomena, including universality properties with respect to the weight distribution and a phase transition in spectral behavior driven by the choice of activation function. Leveraging tools from random matrix theory, we first analyze a class of nonlinear covariance ensembles, which may be of independent interest in high-dimensional statistics, and subsequently characterize the spectral fluctuations of the Jacobian matrix, highlighting their dependence on weight initialization, layer widths, activation functions and network depth. These findings provide finite-width distributional calibration for Jacobian-based stability diagnostics and quantify initialization-induced uncertainty in signal propagation, gradient stability and local sensitivity of high-dimensional neural networks.

Large-scale network inference with uncertainty quantification has important applications in natural, social and medical sciences. The recent work of Fan et al. (J. R. Stat. Soc. Ser. B. Stat. Methodol. 84 (2022) 630–653) introduced a general framework of statistical inference on membership profiles in large networks (SIMPLE) for testing the sharp null hypothesis that a pair of given nodes share the same membership profile. In real applications, there are often groups of nodes under investigation that may share similar membership profiles in the presence of relatively weaker signals than the setting considered in SIMPLE. To address these practical challenges, in this paper, we propose a SIMPLE method with random coupling (SIMPLE-RC) for testing the nonsharp composite null hypothesis that a group of given nodes shares similar (not necessarily identical) membership profiles under weaker signals. Utilizing the idea of random coupling, we construct our test as the maximum of the SIMPLE tests for subsampled node pairs from the group. Such a technique significantly reduces the correlation among individual SIMPLE tests while largely maintaining the power, enabling delicate analysis of the asymptotic distributions of the SIMPLE-RC test. Our method and theory cover both the cases with and without node degree heterogeneity. These new theoretical developments are empowered by a second-order expansion of spiked eigenvectors under the ${\ell }_{\infty }$-norm, built upon our work for random matrices with weak spikes. Our theoretical results and the practical advantages of the newly proposed method are demonstrated through simulations and real data examples.

We consider two hypothesis testing problems for low-rank and high-dimensional tensor signals, namely the tensor signal alignment and tensor signal matching problems. These problems are challenging due to the high dimension of tensors and the lack of suitable test statistics. By exploiting a recent tensor contraction method, we propose and validate relevant test statistics using eigenvalues of a data matrix resulting from the tensor contraction. The matrix entries exhibit long-range dependence, which makes the analysis of the matrix challenging, involved, and distinct from standard random matrix theory. Our approach provides a novel framework for addressing hypothesis testing problems in the context of high-dimensional tensor signals.

In this paper, we study the problem of testing the equality of two multivariate distributions. One class of tests used for this purpose utilizes geometric graphs constructed using interpoint distances. So far, the asymptotic theory of these tests applies only to graphs, which fall under the stabilizing graphs framework of Penrose and Yukich (Ann. Appl. Probab. 13 (2003) 277–303). We study the case of the K-nearest neighbors graph where $\mathit{K}={\mathit{k}}_{\mathit{N}}$ increases with the sample size, which does not fall under the stabilizing graphs framework. Our main result gives detection thresholds for this test in parametrized families when ${\mathit{k}}_{\mathit{N}}=\mathit{o}({\mathit{N}}^{1/4})$, thus extending the family of graphs where the theoretical behavior is known. We propose a 2-sided version of the test, which removes an exponent gap that plagues the 1-sided test. Our result also shows that increasing the number of nearest neighbors improves the power of the test. This provides a theoretical justification for using denser graphs in testing the equality of two distributions.

Using i.i.d. data to estimate a high-dimensional distribution in Wasserstein distance is a fundamental instance of the curse of dimensionality. We explore how structural knowledge about the data-generating process, which gives rise to the distribution can be used to overcome this curse. More precisely, we work with the set of distributions of probabilistic graphical models for a given directed acyclic graph. It turns out that this knowledge is only helpful if it can be quantified, which we formalize via smoothness conditions on the transition kernels in the disintegration corresponding to the graph. In this case, we prove that the rate of estimation is governed by the local structure of the graph, more precisely by dimensions corresponding to single nodes together with their parent nodes. The precise rate depends on the exact notion of smoothness assumed for the kernels, where either weak (Wasserstein–Lipschitz) or strong (bidirectional total variation Lipschitz) conditions lead to different results. We prove sharpness under the strong condition and show that this condition covers, as a special case, distributions having a positive Lipschitz density.

Crossover designs randomly assign each unit to receive a sequence of treatments. By comparing outcomes within the same unit, these designs can effectively eliminate between-unit variation and facilitate the identification of both instantaneous effects of current treatments and carryover effects from past treatments. They are widely used in traditional biomedical studies and are increasingly adopted in modern digital platforms. However, standard analyses of crossover designs often rely on strong parametric models, making inference vulnerable to model misspecification. This paper adopts the design-based framework to analyze general crossover designs. We make two main contributions. First, we use potential outcomes to formally define the causal estimands and assumptions on the data-generating process. For any given type of crossover design and assumptions on potential outcomes, we outline a procedure for identification and estimation, emphasizing the central role of the treatment assignment mechanism in the design-based inference. Second, we unify the analysis of crossover designs using least squares, with restrictions on the coefficients and weights on the units. Based on the theory, we recommend the specification of the regression function, weighting scheme, and coefficient restrictions to assess identifiability, construct efficient estimators, and conduct inference in a unified fashion. Crucially, the least squares procedure is simple to implement, and yields not only consistent and efficient point estimates but also valid inference even when the working regression model is misspecified.

Multiview datasets are common in scientific and engineering applications, yet existing fusion methods offer limited theoretical guarantees, particularly in the presence of heterogeneous and high-dimensional noise. We propose Generalized Robust Adaptive-Bandwidth Multiview Diffusion Maps (GRAB-MDM), a new kernel-based diffusion geometry framework for integrating multiple noisy data sources. The key innovation of GRAB-MDM is a view-dependent bandwidth selection strategy that adapts to the geometry and noise level of each view, enabling a stable and principled construction of multiview diffusion operators. Under a common-manifold model, we establish asymptotic convergence results and show that the adaptive bandwidths lead to provably robust recovery of the shared intrinsic structure, even when noise levels and sensor dimensions differ across views. Numerical experiments and two real data analyses demonstrate that GRAB-MDM significantly improves robustness and embedding quality compared with existing algorithms. The proposed framework offers a practical and theoretically grounded solution for multiview sensor fusion in high-dimensional noisy environments.