This paper constructs a new estimator for large covariance matrices by drawing a bridge between the classic (Stein (1975)) estimator in finite samples and recent progress under large-dimensional asymptotics. The estimator keeps the eigenvectors of the sample covariance matrix and applies shrinkage to the inverse sample eigenvalues. The corresponding formula is quadratic: it has two shrinkage targets weighted by quadratic functions of the concentration (that is, matrix dimension divided by sample size). The first target dominates mid-level concentrations and the second one higher levels. This extra degree of freedom enables us to outperform linear shrinkage when the optimal shrinkage is not linear, which is the general case. Both of our targets are based on what we term the “Stein shrinker”, a local attraction operator that pulls sample covariance matrix eigenvalues towards their nearest neighbors, but whose force diminishes with distance (like gravitation). We prove that no cubic or higher-order nonlinearities beat quadratic with respect to Frobenius loss under large-dimensional asymptotics. Non-normality and the case where the matrix dimension exceeds the sample size are accommodated. Monte Carlo simulations confirm state-of-the-art performance in terms of accuracy, speed, and scalability.

We establish a Berry–Esseen bound for general multivariate nonlinear statistics by developing a new multivariate-type randomized concentration inequality. The bound is the best possible for many known statistics. As applications, Berry–Esseen bounds for M-estimators and averaged stochastic gradient descent algorithms are obtained.

Discretizations of the Langevin diffusion have been proven very useful for developing and analyzing algorithms for sampling and stochastic optimization. We present an improved non-asymptotic analysis of the Euler-Maruyama discretization of the Langevin diffusion. Our analysis does not require global contractivity, and yields polynomial dependence on the time horizon. Compared to existing approaches, we make an additional smoothness assumption, and improve the existing rate in discretization step size from $O(\mathrm{\eta })$ to $O({\mathrm{\eta }}^2)$ in terms of the KL divergence. This result matches the correct order for numerical SDEs, without suffering from exponential time dependence. When applied to MCMC, this result simultaneously improves on the analyses of a range of sampling algorithms that are based on Dalalyan’s approach.

The branching process in random environment ${\{Z_n\}}_{n\ge 0}$ is a population growth process where individuals reproduce independently of each other with the reproduction law randomly picked at each generation. We focus on the supercritical case, when the process survives with positive probability and grows exponentially fast on the nonextinction set. Using Fourier techniques we improve existing central limit theorem as well as we obtain Edgeworth expansions and the renewal theorem for the sequence ${\{log{Z}_n\}}_{n\ge 0}$. The strategy is to compare $log{Z}_n$ with partial sums of i.i.d. random variables in order to obtain precise estimates.

In certain non-linear regression models, the functional form of the regression function is not explicitly available, but is only described by a set of differential equations. For regression models described by a set of ordinary differential equations (ODEs), both Bayesian and non-Bayesian methods for inference were developed in the literature. In this paper, we consider a Bayesian approach to non-linear regression with respect to a multidimensional predictor variable given by a set of partial differential equations (PDEs). We consider a computationally convenient two-step approach by first representing the functions nonparametrically, constructing a finite random series prior using tensor products of B-splines and directly inducing a posterior distribution on parameter space through an appropriate projection map. By considering three different choices of the projection map, we propose three different approaches with their merits. We allow generalized non-linear regression with the response variable following an exponential family of distributions, extending the method beyond regression with additive normal errors. We establish Bernstein-von Mises type theorems which show $\sqrt{n}$-consistency and asymptotically correct frequentist coverage of Bayesian credible regions. We also conduct a simulation study to evaluate the finite sample performances of the proposed methods.

We study the “periodic homogenization” for a class of nonlocal partial differential equations of parabolic-type with rapidly oscillating coefficients, related to stochastic differential equations driven by multiplicative isotropic α-stable Lévy noise ($1<\mathit{\alpha }<2$) which is nonlinear in the noise component. Our homogenization method is probabilistic. It turns out that, under suitable regularity assumptions, the limit of the solutions satisfies a nonlocal partial differential equation with constant coefficients, which are associated to a symmetric α-stable Lévy process.

Consider semiparametric estimation where a doubly robust estimating function for a low-dimensional parameter is available, depending on two working models. With high-dimensional data, we develop regularized calibrated estimation as a general method for estimating the parameters in the two working models, such that valid Wald confidence intervals can be obtained for the parameter of interest under suitable sparsity conditions if either of the two working models is correctly specified. We propose a computationally tractable two-step algorithm and provide rigorous theoretical analysis which justifies sufficiently fast rates of convergence for the regularized calibrated estimators in spite of sequential construction and establishes a desired asymptotic expansion for the doubly robust estimator. As concrete examples, we discuss applications to partially linear, log-linear, and logistic models and estimation of average treatment effects. Numerical studies in the former three examples demonstrate superior performance of our method, compared with debiased Lasso.

We design confidence sequences—sequences of confidence intervals which are valid uniformly over time—for quantiles of any distribution over a complete, fully-ordered set, based on a stream of i.i.d. observations. We give methods both for tracking a fixed quantile and for tracking all quantiles simultaneously. Specifically, we provide explicit expressions with small constants for intervals whose widths shrink at the fastest possible $\sqrt{t^{-1}loglogt}$ rate, along with a nonasymptotic concentration inequality for the empirical distribution function which holds uniformly over time with the same rate. The latter strengthens Smirnov’s empirical process law of the iterated logarithm and extends the Dvoretzky-Kiefer-Wolfowitz inequality to hold uniformly over time. We give a new algorithm and sample complexity bound for selecting an arm with an approximately best quantile in a multi-armed bandit framework. In simulations, our method needs fewer samples than existing methods by a factor of five to fifty.

In this paper, we establish some new central limit theorems for certain spectral statistics of a high-dimensional sample covariance matrix under a divergent spectral norm population model. This model covers the divergent spiked population model as a special case. Meanwhile, the number of the spiked eigenvalues can either be fixed or grow to infinity. It is seen from our theorems that the divergence of population spectral norm affects the fluctuations of the linear spectral statistics in a fickle way, depending on the divergence rate.

We solve multidimensional SDEs with distributional drift driven by symmetric, α-stable Lévy processes for $\mathit{\alpha }\in (1,2]$ by studying the associated (singular) martingale problem and by solving the Kolmogorov backward equation. We allow for drifts of regularity $(2-2\mathit{\alpha })∕3$, and in particular we go beyond the by now well understood “Young regime”, where the drift must have better regularity than $(1-\mathit{\alpha })∕2$. The analysis of the Kolmogorov backward equation in the low regularity regime is based on paracontrolled distributions. As an application of our results we construct a Brox diffusion with Lévy noise.

We consider a one-dimensional classical Wigner jellium, not necessarily charge neutral, for which the electrons are allowed to exist beyond the support of the background charge. The model can be seen as a one-dimensional Coulomb gas in which the external field is generated by a smeared background on an interval. It is a true one-dimensional Coulomb gas and not a one-dimensional log-gas. The system exists if and only if the total background charge is greater than the number of electrons minus one. For various backgrounds, we show convergence to point processes, at the edge of the support of the background. In particular, this provides asymptotic analysis of the fluctuations of the right-most particle. Our analysis reveals that these fluctuations are not universal, in the sense that depending on the background, the tails range anywhere from exponential to Gaussian-like behavior, including for instance Tracy – Widom-like behavior. We also obtain a Rényi-type probabilistic representation for the order statistics of the particle system beyond the support of the background.

Discretizations of Langevin diffusions provide a powerful method for sampling and Bayesian inference. However, such discretizations require evaluation of the gradient of the potential function. In several real-world scenarios, obtaining gradient evaluations might either be computationally expensive, or simply impossible. In this work, we propose and analyze stochastic zeroth-order sampling algorithms for discretizing overdamped and underdamped Langevin diffusions. Our approach is based on estimating the gradients, based on Gaussian Stein’s identities, widely used in the stochastic optimization literature. We provide a comprehensive oracle complexity analysis – number noisy function evaluations to be made to obtain an ϵ-approximate sample in Wasserstein distance – of stochastic zeroth-order discretizations of both overdamped and underdamped Langevin diffusions, under various noise models. Our theoretical contributions extend the applicability of sampling algorithms to the noisy black-box settings arising in practice.

We study the distribution of the maximum likelihood estimate (MLE) in high-dimensional logistic models, where covariates are Gaussian with an arbitrary covariance structure. We prove that in the limit of large problems holding the ratio between the number p of covariates and the sample size n constant, every finite list of MLE coordinates follows a multivariate normal distribution. Concretely, the jth coordinate ${\hat{\mathit{\beta }}}_j$ of the MLE is asymptotically normally distributed with mean ${\mathit{\alpha }}_{\star }{\mathit{\beta }}_j$ and standard deviation ${\mathit{\sigma }}_{\star }∕{\mathit{\tau }}_j$; here, ${\mathit{\beta }}_j$ is the value of the true regression coefficient, and ${\mathit{\tau }}_j$ the standard deviation of the jth predictor conditional on all the others. The numerical parameters ${\mathit{\alpha }}_{\star }>1$ and ${\mathit{\sigma }}_{\star }$ only depend upon the problem dimensionality $p∕n$ and the overall signal strength, and can be accurately estimated. Our results imply that the MLE’s magnitude is biased upwards and that the MLE’s standard deviation is greater than that predicted by classical theory. We present a series of experiments on simulated and real data showing excellent agreement with the theory.

Conditional Variance Estimation (CVE) is a novel sufficient dimension reduction (SDR) method for additive error regressions with continuous predictors and link function. It operates under the assumption that the predictors can be replaced by a lower dimensional projection without loss of information. Conditional Variance Estimation is fully data driven, does not require the restrictive linearity and constant variance conditions, and is not based on inverse regression as the majority of moment and likelihood based sufficient dimension reduction methods. CVE is shown to be consistent and its objective function to be uniformly convergent. CVE outperforms the mean average variance estimation, (MAVE), its main competitor, in several simulation settings, remains on par under others, while it always outperforms inverse regression based linear SDR methods, such as Sliced Inverse Regression.

We prove an exponential decay concentration inequality to bound the tail probability of the difference between the log-likelihood of discrete random variables on a finite alphabet and the negative entropy. The concentration bound we derive holds uniformly over all parameter values. The new result improves the convergence rate in an earlier result of Zhao (2020), from $(K^{2}logK)∕n=o(1)$ to ${(logK)}^2∕n=o(1)$, where n is the sample size and K is the size of the alphabet. We further prove that the rate ${(logK)}^2∕n=o(1)$ is optimal. The result is extended to misspecified log-likelihoods for grouped random variables. We give applications of the new result in information theory.

By substituting an $L_p$ loss function for the $L_1$ loss function in the optimization problem defining quantiles, one obtains $L_p$-quantiles that, as shown recently, dominate their classical $L_1$-counterparts in financial risk assessment. In this work, we propose a concept of multivariate $L_p$-quantiles generalizing the spatial ($L_1$-)quantiles introduced by Probal Chaudhuri (J. Amer. Statist. Assoc. 91 (1996) 862–872). Rather than restricting to power loss functions, we actually allow for a large class of convex loss functions ρ. We carefully study existence and uniqueness of the resulting ρ-quantiles, both for a general probability measure over ${\mathbb{R}}^d$ and for a spherically symmetric one. Interestingly, the results crucially depend on ρ and on the nature of the underlying probability measure. Building on an investigation of the differentiability properties of the objective function defining ρ-quantiles, we introduce a companion concept of spatial ρ-depth, that generalizes the classical spatial depth. We study extreme ρ-quantiles and show in particular that extreme $L_p$-quantiles behave in fundamentally different ways for $p\le 2$ and $p>2$. Finally, we establish Bahadur representation results for sample ρ-quantiles and derive their asymptotic distributions. Throughout, we impose only very mild assumptions on the underlying probability measure, and in particular we never assume absolute continuity with respect to the Lebesgue measure.

The exchange algorithm is one of the most popular extensions of the Metropolis–Hastings algorithm to sample from doubly-intractable distributions. However, the theoretical exploration of the exchange algorithm is very limited. For example, natural questions like ‘Does exchange algorithm converge at a geometric rate?’ or ‘Does the exchange algorithm admit a Central Limit Theorem?’ have not been answered yet. In this paper, we study the theoretical properties of the exchange algorithm, in terms of asymptotic variance and convergence speed. We compare the exchange algorithm with the original Metropolis–Hastings algorithm and provide both necessary and sufficient conditions for the geometric ergodicity of the exchange algorithm. Moreover, we prove that our results can be applied to various practical applications such as location models, Gaussian models, Poisson models, and a large class of exponential families, which includes most of the practical applications of the exchange algorithm. A central limit theorem for the exchange algorithm is also established. Our results justify the theoretical usefulness of the exchange algorithm.

We consider a family of continuous processes ${\{X^{\mathit{\epsilon }}\}}_{\mathit{\epsilon }>0}$ which are measurable with respect to a white noise measure, take values in the space of continuous functions $C({[0,1]}^d:\mathbb{R})$, and have the Wiener chaos expansion

$$X^{\mathit{\epsilon }}=\sum _{n=0}^{\mathrm{\infty }}{\mathit{\epsilon }}^n{I}_n\left(f_n^{\mathit{\epsilon }}\right).$$

We provide sufficient conditions for the large deviations principle of ${\{X^{\mathit{\epsilon }}\}}_{\mathit{\epsilon }>0}$ to hold in $C({[0,1]}^d:\mathbb{R})$, thereby refreshing a problem left open by Pérez–Abreu (1993) in the Brownian motion case. The proof is based on the weak convergence approach to large deviations: it involves demonstrating the convergence in distribution of certain perturbations of the original process, and thus the main difficulties lie in analysing and controlling the perturbed multiple stochastic integrals. Moreover, adopting this representation offers a new perspective on pathwise large deviations and induces a variety of applications thereof.

Markov chain Monte Carlo (MCMC) methods, specifically samplers based on random walks, often have difficulty handling target distributions with complex geometry such as multi-modality. We propose an adaptive multiple-try Metropolis algorithm designed to tackle such problems by combining the flexibility of multiple-proposal samplers with the user-friendliness and optimality of adaptive algorithms. We prove the ergodicity of the resulting Markov chain with respect to the target distribution using common techniques in the adaptive MCMC literature. In a Bayesian model for loss of heterozygosity in cancer cells, we find that our method outperforms traditional adaptive samplers, non-adaptive multiple-try Metropolis samplers, and various more sophisticated competing methods.

A well-known difficult problem regarding Metropolis-Hastings algorithms is to get sharp bounds on their convergence rates. Moreover, a fundamental but often overlooked problem in Markov chain theory is to study the convergence rates for different initializations. In this paper, we study the two issues mentioned above of the Independent Metropolis-Hastings (IMH) algorithms on both general and discrete state spaces. We derive the exact convergence rate and prove that the IMH algorithm’s different deterministic initializations have the same convergence rate. We get the exact convergence speed for IMH algorithms on general state spaces.

In this paper a new additive regression technique is developed for response variables that take values in general Hilbert spaces. The proposed method is based on the idea of smooth backfitting that has been developed mainly for real-valued responses. The local polynomial smoothing device is adopted, which renders various advantages of the technique evidenced in the classical univariate kernel regression with real-valued responses. It is demonstrated that the new technique eliminates many limitations which existing methods are subject to. In contrast to the existing techniques, the proposed approach is equipped with the estimation of the derivatives as well as the regression function itself, and provides options to make the estimated regression function free from boundary effects and possess oracle properties. A comprehensive theory is presented for the proposed method, which includes the rates of convergence in various modes and the asymptotic distributions of the estimators. The efficiency of the proposed method is also demonstrated via simulation study and is illustrated through real data applications.

We study systems of three interacting particles, in which drifts and variances are assigned by rank. These systems are degenerate: the variances corresponding to one or two ranks can vanish, so the corresponding ranked motions become ballistic rather than diffusive. Depending on which ranks are allowed to “go ballistic” the systems exhibit markedly different behavior, which we study here in some detail. Also studied are stability properties for the resulting planar process of gaps between successive ranks.

The degrees of freedom quantify the number of parameters of an estimator and is central to various risk minimization procedures. Its computation and properties are well understood when dealing with finite dimensional linear models, possibly regularized using sparsity through the celebrated Lasso estimator. However, for some applications (such as super-resolution in imaging or training multi-layer perceptrons with a single hidden layer) it makes sense to rather consider “continuous” methods. For these applications, these methods avoid the discretization of the parameter space and lead to more efficient numerical solvers. Training these continuous models with a sparsity inducing prior can be achieved by solving a convex optimization problem over the infinite dimensional space of measures, which is often called the Beurling Lasso (Blasso), and is the continuous counterpart of the Lasso. Previous works (Ann. Statist. 35 (2007) 2173–2192; Statist. Sinica 23 (2013) 809–828) show that the size of the smallest solution support is an unbiased estimator of the degrees of freedom of the Lasso. Our main contribution is a proof of a continuous counterpart to this result for the Blasso. In contrast to the Lasso, our new formula shows that the empirical degrees of freedom of the Blasso is not proportional to the size of the support. Our second contribution is a detailed study of the case of sampling Fourier coefficients in 1D, which corresponds to a super-resolution problem. We show that our formula for the unbiased estimation of the degrees of freedom is valid outside of a set of measure zero of observations, which in turn justifies its use to compute an unbiased estimator of the prediction risk using the Stein Unbiased Risk Estimator (SURE). We also report numerical results for both the case of Fourier sampling and the learning of a multilayers perceptron with a single hidden layer.