Understanding statistical inference under possibly nonsparse high-dimensional models has gained much interest recently. For a given component of the regression coefficient, we show that the difficulty of the problem depends on the sparsity of the corresponding row of the precision matrix of the covariates, not the sparsity of the regression coefficients. We develop new concepts of uniform and essentially uniform nontestability that allow the study of limitations of tests across a broad set of alternatives. Uniform nontestability identifies a collection of alternatives such that the power of any test, against any alternative in the group, is asymptotically at most equal to the nominal size. Implications of the new constructions include new minimax testability results that, in sharp contrast to the current results, do not depend on the sparsity of the regression parameters. We identify new tradeoffs between testability and feature correlation. In particular, we show that, in models with weak feature correlations, minimax lower bound can be attained by a test whose power has the $\sqrt{\mathit{n}}$ rate, regardless of the size of the model sparsity.

A prominent concern of scientific investigators is the presence of unobserved hidden variables in association analysis. Ignoring hidden variables often yields biased statistical results and misleading scientific conclusions. Motivated by this practical issue, this paper studies the multivariate response regression with hidden variables, $\mathit{Y}={({\mathrm{\Psi }}^{\ast })}^{\mathit{T}}\mathit{X}+{({\mathit{B}}^{\ast })}^{\mathit{T}}\mathit{Z}+\mathit{E}$, where $\mathit{Y}\in {\mathbb{R}}^{\mathit{m}}$ is the response vector, $\mathit{X}\in {\mathbb{R}}^{\mathit{p}}$ is the observable feature, $\mathit{Z}\in {\mathbb{R}}^{\mathit{K}}$ represents the vector of unobserved hidden variables, possibly correlated with X, and E is an independent error. The number of hidden variables K is unknown and both m and p are allowed, but not required, to grow with the sample size n.

Though ${\mathrm{\Psi }}^{\ast }$ is shown to be nonidentifiable due to the presence of hidden variables, we propose to identify the projection of ${\mathrm{\Psi }}^{\ast }$ onto the orthogonal complement of the row space of ${\mathit{B}}^{\ast }$, denoted by ${\mathrm{\Theta }}^{\ast }$. The quantity ${({\mathrm{\Theta }}^{\ast })}^{\mathit{T}}\mathit{X}$ measures the effect of X on Y that cannot be explained through the hidden variables, and thus ${\mathrm{\Theta }}^{\ast }$ is treated as the parameter of interest. Motivated by the identifiability proof, we propose a novel estimation algorithm for ${\mathrm{\Theta }}^{\ast }$, called HIVE, under homoscedastic errors. The first step of the algorithm estimates the best linear prediction of Y given X, in which the unknown coefficient matrix exhibits an additive decomposition of ${\mathrm{\Psi }}^{\ast }$ and a dense matrix due to the correlation between X and Z. Under the sparsity assumption on ${\mathrm{\Psi }}^{\ast }$, we propose to minimize a penalized least squares loss by regularizing ${\mathrm{\Psi }}^{\ast }$ and the dense matrix via group-lasso and multivariate ridge, respectively. Nonasymptotic deviation bounds of the in-sample prediction error are established. Our second step estimates the row space of ${\mathit{B}}^{\ast }$ by leveraging the covariance structure of the residual vector from the first step. In the last step, we estimate ${\mathrm{\Theta }}^{\ast }$ via projecting Y onto the orthogonal complement of the estimated row space of ${\mathit{B}}^{\ast }$ to remove the effect of hidden variables. Nonasymptotic error bounds of our final estimator of ${\mathrm{\Theta }}^{\ast }$, which are valid for any $\mathit{m},\mathit{p},\mathit{K}$ and n, are established. We further show that, under mild assumptions, the rate of our estimator matches the best possible rate with known ${\mathit{B}}^{\ast }$ and is adaptive to the unknown sparsity of ${\mathrm{\Theta }}^{\ast }$ induced by the sparsity of ${\mathrm{\Psi }}^{\ast }$. The model identifiability, estimation algorithm and statistical guarantees are further extended to the setting with heteroscedastic errors. Thorough numerical simulations and two real data examples are provided to back up our theoretical results.

Let ${\{({\mathit{X}}_{\mathit{i}},{\mathit{Y}}_{\mathit{i}})\}}_{\mathit{i}=1}^{\mathit{n}}$ be a sequence of independent bivariate random vectors. In this paper, we establish a refined Cramér-type moderate deviation theorem for the general self-normalized sum ${\sum }_{\mathit{i}=1}^{\mathit{n}}{\mathit{X}}_{\mathit{i}}/{({\sum }_{\mathit{i}=1}^{\mathit{n}}{\mathit{Y}}_{\mathit{i}}^2)}^{1/2}$, which unifies and extends the classical Cramér (Actual. Sci. Ind. 736 (1938) 5–23) theorem and the self-normalized Cramér-type moderate deviation theorems by Jing, Shao and Wang (Ann. Probab. 31 (2003) 2167–2215) as well as the further refined version by Wang (J. Theoret. Probab. 24 (2011) 307–329). The advantage of our result is evidenced through successful applications to weakly dependent random variables and self-normalized winsorized mean. Specifically, by applying our new framework on general self-normalized sum, we significantly improve Cramér-type moderate deviation theorems for one-dependent random variables, geometrically β-mixing random variables and causal processes under geometrical moment contraction. As an additional application, we also derive the Cramér-type moderate deviation theorems for self-normalized winsorized mean.

Distributed minimax estimation and distributed adaptive estimation under communication constraints for Gaussian sequence model and white noise model are studied. The minimax rate of convergence for distributed estimation over a given Besov class, which serves as a benchmark for the cost of adaptation, is established. We then quantify the exact communication cost for adaptation and construct an optimally adaptive procedure for distributed estimation over a range of Besov classes.

The results demonstrate significant differences between nonparametric function estimation in the distributed setting and the conventional centralized setting. For global estimation, adaptation in general cannot be achieved for free in the distributed setting. The new technical tools to obtain the exact characterization for the cost of adaptation can be of independent interest.

Network method of moments (Ann. Statist. 39 (2011) 2280–2301) is an important tool for nonparametric network inference. However, there has been little investigation on accurate descriptions of the sampling distributions of network moment statistics. In this paper, we present the first higher-order accurate approximation to the sampling CDF of a studentized network moment by Edgeworth expansion. In sharp contrast to classical literature on noiseless U-statistics, we show that the Edgeworth expansion of a network moment statistic as a noisy U-statistic can achieve higher-order accuracy without nonlattice or smoothness assumptions but just requiring weak regularity conditions. Behind this result is our surprising discovery that the two typically-hated factors in network analysis, namely, sparsity and edgewise observational errors, jointly play a blessing role, contributing a crucial self-smoothing effect in the network moment statistic and making it analytically tractable. Our assumptions match the minimum requirements in related literature. For sparse networks, our theory shows that our empirical Edgeworth expansion and a simple normal approximation both achieve the same gradually depreciating Berry–Esseen-type bound as the network becomes sparser. This result also significantly refines the best previous theoretical result.

For practitioners, our empirical Edgeworth expansion is highly accurate and computationally efficient. It is also easy to implement and convenient for parallel computing. We demonstrate the clear advantage of our method by several comprehensive simulation studies. As a byproduct, we also provide a finite-sample analysis of the network jackknife.

We showcase three applications of our results in network inference. We prove, to our knowledge, the first theoretical guarantee of higher-order accuracy for some network bootstrap schemes, and moreover, the first theoretical guidance for selecting the subsample size for network subsampling. We also derive a one-sample test and the Cornish–Fisher confidence interval for a given moment with higher-order accurate controls of confidence level and type I error, respectively.

We consider a multivariate response regression model where each coordinate is described by a location-scale non- or semiparametric regression and where the dependence structure of the “noise term” is described by a parametric copula. Our goal is to estimate the associated Euclidean copula parameter, given a sample of the response and the covariate. In the absence of the copula sample, the usual oracle ranks are no longer computable. Instead, we study the normal scores estimator for the Gaussian copula and generalized pseudo-likelihood estimation for general parametric copulas, both based on residual ranks calculated from preliminary non- or semiparametric estimators of the location and scale functions. We show that the residual-based estimators are asymptotically equivalent to their oracle counterparts and provide explicit rate of convergence. Partially to serve this objective, we also study weighted convergence of the residual empirical process under the non- or semiparametric regression model.

This article considers change-point testing and estimation for a sequence of high-dimensional data. In the case of testing for a mean shift for high-dimensional independent data, we propose a new test which is based on U-statistic in Chen and Qin (Ann. Statist. 38 (2010) 808–835) and utilizes the self-normalization principle (Shao J. R. Stat. Soc. Ser. B. Stat. Methodol. 72 (2010) 343–366; Shao and Zhang J. Amer. Statist. Assoc. 105 (2010) 1228–1240). Our test targets dense alternatives in the high-dimensional setting and involves no tuning parameters. To extend to change-point testing for high-dimensional time series, we introduce a trimming parameter and formulate a self-normalized test statistic with trimming to accommodate the weak temporal dependence. On the theory front we derive the limiting distributions of self-normalized test statistics under both the null and alternatives for both independent and dependent high-dimensional data. At the core of our asymptotic theory, we obtain weak convergence of a sequential U-statistic based process for high-dimensional independent data, and weak convergence of sequential trimmed U-statistic based processes for high-dimensional linear processes, both of which are of independent interests. Additionally, we illustrate how our tests can be used in combination with wild binary segmentation to estimate the number and location of multiple change points. Numerical simulations demonstrate the competitiveness of our proposed testing and estimation procedures in comparison with several existing methods in the literature.

Large-scale multiple testing is a fundamental problem in high dimensional statistical inference. It is increasingly common that various types of auxiliary information, reflecting the structural relationship among the hypotheses, are available. Exploiting such auxiliary information can boost statistical power. To this end, we propose a framework based on a two-group mixture model with varying probabilities of being null for different hypotheses a priori, where a shape-constrained relationship is imposed between the auxiliary information and the prior probabilities of being null. An optimal rejection rule is designed to maximize the expected number of true positives when average false discovery rate is controlled. Focusing on the ordered structure, we develop a robust EM algorithm to estimate the prior probabilities of being null and the distribution of p-values under the alternative hypothesis simultaneously. We show that the proposed method has better power than state-of-the-art competitors while controlling the false discovery rate, both empirically and theoretically. Extensive simulations demonstrate the advantage of the proposed method. Datasets from genome-wide association studies are used to illustrate the new methodology.

The Hilbert–Schmidt Independence Criterion (HSIC) is a dependence measure based on reproducing kernel Hilbert spaces that is widely used to test independence between two random vectors. Remains the delicate choice of the kernel. In this work, we develop a new HSIC-based aggregated procedure which avoids such a kernel choice, and provide theoretical guarantees for this procedure. To achieve this, on the one hand, we introduce non-asymptotic single tests based on Gaussian kernels with a given bandwidth, which are of prescribed level. Then, we aggregate several single tests with different bandwidths, and prove sharp upper bounds for the uniform separation rate of the aggregated procedure over Sobolev balls. On the other hand, we provide a lower bound for the non-asymptotic minimax separation rate of testing over Sobolev balls, and deduce that the aggregated procedure is adaptive in the minimax sense over such regularity spaces. Finally, from a practical point of view, we perform numerical studies in order to assess the efficiency of our aggregated procedure and compare it to existing tests in the literature.

We consider a sparse high-dimensional regression model where the goal is to recover a k-sparse unknown binary vector ${\mathit{\beta }}^{\ast }$ from n noisy linear observations of the form $\mathit{Y}=\mathit{X}{\mathit{\beta }}^{\ast }+\mathit{W}\in {\mathbb{R}}^{\mathit{n}}$ where $\mathit{X}\in {\mathbb{R}}^{\mathit{n}\times \mathit{p}}$ has i.i.d. $\mathit{N}(0,1)$ entries and $\mathit{W}\in {\mathbb{R}}^{\mathit{n}}$ has i.i.d. $\mathit{N}(0,{\mathit{\sigma }}^2)$ entries. In the high signal-to-noise ratio regime and sublinear sparsity regime, while the order of the sample size needed to recover the unknown vector information-theoretically is known to be ${\mathit{n}}^{\ast }:=2\mathit{k}log\mathit{p}/log(\mathit{k}/{\mathit{\sigma }}^2+1)$, no polynomial-time algorithm is known to succeed unless $\mathit{n}>{\mathit{n}}_{\mathrm{alg}}:=(2\mathit{k}+{\mathit{\sigma }}^2)log\mathit{p}$.

In this work, we offer a series of results investigating multiple computational and statistical aspects of the recovery task in the regime $\mathit{n}\in [{\mathit{n}}^{\ast },{\mathit{n}}_{\mathrm{alg}}]$. First, we establish a novel information-theoretic property of the MLE of the problem happening around $\mathit{n}={\mathit{n}}^{\ast }$ samples, which we coin as an “all-or-nothing behavior”: when $\mathit{n}>{\mathit{n}}^{\ast }$ it recovers almost perfectly the support of ${\mathit{\beta }}^{\ast }$, while if $\mathit{n}<{\mathit{n}}^{\ast }$ it fails to recover any fraction of it correctly. Second, at an attempt to understand the computational hardness in the regime $\mathit{n}\in [{\mathit{n}}^{\ast },{\mathit{n}}_{\mathrm{alg}}]$, we prove that at order ${\mathit{n}}_{\mathrm{alg}}$ samples there is an Overlap Gap Property (OGP) phase transition occurring at the landscape of the MLE: for constants $\mathit{c},\mathit{C}>0$ when $\mathit{n}<\mathit{c}{\mathit{n}}_{\mathrm{alg}}$ OGP appears in the landscape of MLE while if $\mathit{n}>\mathit{C}{\mathit{n}}_{\mathrm{alg}}$ OGP disappears. OGP is a geometric “disconnectivity” property, which initially appeared in the theory of spin glasses and is known to suggest algorithmic hardness when it occurs. Finally, using certain technical results obtained to establish the OGP phase transition, we additionally establish various novel positive and negative algorithmic results for the recovery task of interest, including the failure of LASSO with access to $\mathit{n}<\mathit{c}{\mathit{n}}_{\mathrm{alg}}$ samples and the success of a simple local search method with access to $\mathit{n}>\mathit{C}{\mathit{n}}_{\mathrm{alg}}$ samples.

Sufficient dimension reduction (SDR) embodies a family of methods that aim for reduction of dimensionality without loss of information in a regression setting. In this article, we propose a new method for nonparametric function-on-function SDR, where both the response and the predictor are a function. We first develop the notions of functional central mean subspace and functional central subspace, which form the population targets of our functional SDR. We then introduce an average Fréchet derivative estimator, which extends the gradient of the regression function to the operator level and enables us to develop estimators for our functional dimension reduction spaces. We show the resulting functional SDR estimators are unbiased and exhaustive, and more importantly, without imposing any distributional assumptions such as the linearity or the constant variance conditions that are commonly imposed by all existing functional SDR methods. We establish the uniform convergence of the estimators for the functional dimension reduction spaces, while allowing both the number of Karhunen–Loève expansions and the intrinsic dimension to diverge with the sample size. We demonstrate the efficacy of the proposed methods through both simulations and two real data examples.

Multiple imputation (MI) is a technique especially designed for handling missing data in public-use datasets. It allows analysts to perform incomplete-data inference straightforwardly by using several already imputed datasets released by the dataset owners. However, the existing MI tests require either a restrictive assumption on the missing-data mechanism, known as equal odds of missing information (EOMI), or an infinite number of imputations. Some of them also require analysts to have access to restrictive or nonstandard computer subroutines. Besides, the existing MI testing procedures cover only Wald’s tests and likelihood ratio tests but not Rao’s score tests, therefore, these MI testing procedures are not general enough. In addition, the MI Wald’s tests and MI likelihood ratio tests are not procedurally identical, so analysts need to resort to distinct algorithms for implementation. In this paper, we propose a general MI procedure, called stacked multiple imputation (SMI), for performing Wald’s tests, likelihood ratio tests and Rao’s score tests by a unified algorithm. SMI requires neither EOMI nor an infinite number of imputations. It is particularly feasible for analysts as they just need to use a complete-data testing device for performing the corresponding incomplete-data test.

Interpolators—estimators that achieve zero training error—have attracted growing attention in machine learning, mainly because state-of-the art neural networks appear to be models of this type. In this paper, we study minimum ${\ell }_2$ norm (“ridgeless”) interpolation least squares regression, focusing on the high-dimensional regime in which the number of unknown parameters p is of the same order as the number of samples n. We consider two different models for the feature distribution: a linear model, where the feature vectors ${\mathit{x}}_{\mathit{i}}\in {\mathbb{R}}^{\mathit{p}}$ are obtained by applying a linear transform to a vector of i.i.d. entries, ${\mathit{x}}_{\mathit{i}}={\mathrm{\Sigma }}^{1/2}{\mathit{z}}_{\mathit{i}}$ (with ${\mathit{z}}_{\mathit{i}}\in {\mathbb{R}}^{\mathit{p}}$); and a nonlinear model, where the feature vectors are obtained by passing the input through a random one-layer neural network, ${\mathit{x}}_{\mathit{i}}=\mathit{\phi }(\mathit{W}{\mathit{z}}_{\mathit{i}})$ (with ${\mathit{z}}_{\mathit{i}}\in {\mathbb{R}}^{\mathit{d}}$, $\mathit{W}\in {\mathbb{R}}^{\mathit{p}\times \mathit{d}}$ a matrix of i.i.d. entries, and φ an activation function acting componentwise on $\mathit{W}{\mathit{z}}_{\mathit{i}}$). We recover—in a precise quantitative way—several phenomena that have been observed in large-scale neural networks and kernel machines, including the “double descent” behavior of the prediction risk, and the potential benefits of overparametrization.

Network sampling is an indispensable tool for understanding features of large complex networks where it is practically impossible to search over the entire graph. In this paper, we develop a framework for statistical inference for counting network motifs, such as edges, triangles and wedges, in the widely used subgraph sampling model, where each vertex is sampled independently, and the subgraph induced by the sampled vertices is observed. We derive necessary and sufficient conditions for the consistency and the asymptotic normality of the natural Horvitz–Thompson (HT) estimator, which can be used for constructing confidence intervals and hypothesis testing for the motif counts based on the sampled graph. In particular, we show that the asymptotic normality of the HT estimator exhibits an interesting fourth-moment phenomenon, which asserts that the HT estimator (appropriately centered and rescaled) converges in distribution to the standard normal whenever its fourth-moment converges to 3 (the fourth-moment of the standard normal distribution). As a consequence, we derive the exact thresholds for consistency and asymptotic normality of the HT estimator in various natural graph ensembles, such as sparse graphs with bounded degree, Erdős–Rényi random graphs, random regular graphs and dense graphons.

In this paper, we study multivariate ranks and quantiles, defined using the theory of optimal transport, and build on the work of Chernozhukov et al. (Ann. Statist. 45 (2017) 223–256) and Hallin et al. (Ann. Statist. 49 (2021) 1139–1165). We study the characterization, computation and properties of the multivariate rank and quantile functions and their empirical counterparts. We derive the uniform consistency of these empirical estimates to their population versions, under certain assumptions. In fact, we prove a Glivenko–Cantelli type theorem that shows the asymptotic stability of the empirical rank map in any direction. Under mild structural assumptions, we provide global and local rates of convergence of the empirical quantile and rank maps. We also provide a sub-Gaussian tail bound for the global ${\mathit{L}}_2$-loss of the empirical quantile function. Further, we propose tuning parameter-free multivariate nonparametric tests—a two-sample test and a test for mutual independence—based on our notion of multivariate quantiles/ranks. Asymptotic consistency of these tests are shown and the rates of convergence of the associated test statistics are derived, both under the null and alternative hypotheses.

Optimal linear prediction (aka. kriging) of a random field ${\{\mathit{Z}(\mathit{x})\}}_{\mathit{x}\in \mathcal{X}}$ indexed by a compact metric space $(\mathcal{X},{\mathit{d}}_{\mathcal{X}})$ can be obtained if the mean value function $\mathit{m}:\mathcal{X}\to \mathbb{R}$ and the covariance function $\mathit{\varrho }:\mathcal{X}\times \mathcal{X}\to \mathbb{R}$ of Z are known. We consider the problem of predicting the value of $\mathit{Z}({\mathit{x}}^{\ast })$ at some location ${\mathit{x}}^{\ast }\in \mathcal{X}$ based on observations at locations ${\{{\mathit{x}}_{\mathit{j}}\}}_{\mathit{j}=1}^{\mathit{n}}$, which accumulate at ${\mathit{x}}^{\ast }$ as $\mathit{n}\to \infty$ (or, more generally, predicting $\mathit{\phi }(\mathit{Z})$ based on ${\{{\mathit{\phi }}_{\mathit{j}}(\mathit{Z})\}}_{\mathit{j}=1}^{\mathit{n}}$ for linear functionals $\mathit{\phi },{\mathit{\phi }}_1,\dots ,{\mathit{\phi }}_{\mathit{n}}$). Our main result characterizes the asymptotic performance of linear predictors (as n increases) based on an incorrect second-order structure $(\tilde{\mathit{m}},\tilde{\mathit{\varrho }})$, without any restrictive assumptions on ϱ, $\tilde{\mathit{\varrho }}$ such as stationarity. We, for the first time, provide necessary and sufficient conditions on $(\tilde{\mathit{m}},\tilde{\mathit{\varrho }})$ for asymptotic optimality of the corresponding linear predictor holding uniformly with respect to φ. These general results are illustrated by weakly stationary random fields on $\mathcal{X}\subset {\mathbb{R}}^{\mathit{d}}$ with Matérn or periodic covariance functions, and on the sphere $\mathcal{X}={\mathbb{S}}^2$ for the case of two isotropic covariance functions.

We propose a general modeling and algorithmic framework for discrete structure recovery that can be applied to a wide range of problems. Under this framework, we are able to study the recovery of clustering labels, ranks of players, signs of regression coefficients, cyclic shifts and even group elements from a unified perspective. A simple iterative algorithm is proposed for discrete structure recovery, which generalizes methods including Lloyd’s algorithm and the power method. A linear convergence result for the proposed algorithm is established in this paper under appropriate abstract conditions on stochastic errors and initialization. We illustrate our general theory by applying it on several representative problems: (1) clustering in Gaussian mixture model, (2) approximate ranking, (3) sign recovery in compressed sensing, (4) multireference alignment and (5) group synchronization, and show that minimax rate is achieved in each case.

Classical multiple testing theory prescribes the null distribution, which is often too stringent an assumption for nowadays large scale experiments. This paper presents theoretical foundations to understand the limitations caused by ignoring the null distribution, and how it can be properly learned from the same data set, when possible. We explore this issue in the setting where the null distributions are Gaussian with unknown rescaling parameters (mean and variance) whereas the alternative distributions are let arbitrary. In that case, an oracle procedure is the Benjamini–Hochberg procedure applied with the true (unknown) null distribution and we aim at building a procedure that asymptotically mimics the performances of the oracle (AMO in short). Our main result establishes a phase transition at the sparsity boundary $\mathit{n}/log(\mathit{n})$: an AMO procedure exists if and only if the number of false nulls is of order less than $\mathit{n}/log(\mathit{n})$, where n is the total number of tests. Further sparsity boundaries are derived for general location models where the shape of the null distribution is not necessarily Gaussian. In light of our impossibility results, we also pursue the less stringent aim of building a nonparametric confidence region for the null distribution. From a practical perspective, this provides goodness-of-fit tests for the null distribution and allows to assess the reliability of empirical null procedures via novel diagnostic graphs. Our results are illustrated on numerical experiments and real data sets, as detailed in a companion vignette (Roquain and Verzelen (2021)).

We consider a testing problem for cross-sectional independence for high-dimensional panel data, where the number of cross-sectional units is potentially much larger than the number of observations. The cross-sectional independence is described through linear regression models. We study three tests named the sum, the max and the max-sum tests, where the latter two are new. The sum test is initially proposed by Breusch and Pagan (1980). We design the max and sum tests for sparse and nonsparse correlation coefficients of random errors between the linear regression models, respectively. And the max-sum test is devised to compromise both situations on the correlation coefficients. Indeed, our simulation shows that the max-sum test outperforms the previous two tests. This makes the max-sum test very useful in practice where sparsity or not for a set of numbers is usually vague. Toward the theoretical analysis of the three tests, we have settled two conjectures regarding the sum of squares of sample correlation coefficients asked by Pesaran (2004 and 2008). In addition, we establish the asymptotic theory for maxima of sample correlation coefficients appeared in the linear regression model for panel data, which is also the first successful attempt to our knowledge. To study the max-sum test, we create a novel method to show asymptotic independence between maxima and sums of dependent random variables. We expect the method itself is useful for other problems of this nature. Finally, an extensive simulation study as well as a case study are carried out. They demonstrate advantages of our proposed methods in terms of both empirical powers and robustness for correlation coefficients of residuals regardless of sparsity or not.

In this paper, we study the asymptotic behavior of the extreme eigenvalues and eigenvectors of the high-dimensional spiked sample covariance matrices, in the supercritical case when a reliable detection of spikes is possible. In particular, we derive the joint distribution of the extreme eigenvalues and the generalized components of the associated eigenvectors, that is, the projections of the eigenvectors onto arbitrary given direction, assuming that the dimension and sample size are comparably large. In general, the joint distribution is given in terms of linear combinations of finitely many Gaussian and Chi-square variables, with parameters depending on the projection direction and the spikes. Our assumption on the spikes is fully general. First, the strengths of spikes are only required to be slightly above the critical threshold and no upper bound on the strengths is needed. Second, multiple spikes, that is, spikes with the same strength, are allowed. Third, no structural assumption is imposed on the spikes. Thanks to the general setting, we can then apply the results to various high dimensional statistical hypothesis testing problems involving both the eigenvalues and eigenvectors. Specifically, we propose accurate and powerful statistics to conduct hypothesis testing on the principal components. These statistics are data-dependent and adaptive to the underlying true spikes. Numerical simulations also confirm the accuracy and powerfulness of our proposed statistics and illustrate significantly better performance compared to the existing methods in the literature. In particular, our methods are accurate and powerful even when either the spikes are small or the dimension is large.

The split-plot design arose from agricultural science with experimental units, also known as the subplots, nested within groups known as the whole plots. It assigns different interventions at the whole-plot and subplot levels, respectively, providing a convenient way to accommodate hard-to-change factors. By design, subplots within the same whole plot receive the same level of the whole-plot intervention, and thereby induce a group structure on the final treatment assignments. A common strategy is to run an ordinary least squares (ols) regression of the outcome on the treatment indicators coupled with the robust standard errors clustered at the whole-plot level. It does not give consistent estimators for the treatment effects of interest when the whole-plot sizes vary. Another common strategy is to fit a linear mixed-effects model of the outcome with normal random effects and errors. It is a purely model-based approach and can be sensitive to violations of the parametric assumptions. In contrast, design-based inference assumes no outcome models and relies solely on the controllable randomization mechanism determined by the physical experiment. We first extend the existing design-based inference based on the Horvitz–Thompson estimator to the Hajek estimator, and establish the finite-population central limit theorem for both under split-plot randomization. We then reconcile the results with those under the model-based approach, and propose two regression strategies, namely (i) the weighted least squares (wls) fit of the unit-level data based on the inverse probability weighting and (ii) the ols fit of the aggregate data based on whole-plot total outcomes, to reproduce the Hajek and Horvitz–Thompson estimators, respectively. This, together with the asymptotic conservativeness of the corresponding cluster-robust covariances for estimating the true design-based covariances as we establish in the process, justifies the validity of the regression estimators for design-based inference. In light of the flexibility of regression formulation for covariate adjustment, we further extend the theory to the case with covariates, and demonstrate the efficiency gain by regression-based covariate adjustment via both asymptotic theory and simulation. Importantly, all our theories are either numeric or design-based, and hold regardless of how well the regression equations represent the true data generating process.

The goal of this paper is to show that a single robust estimator of the mean of a multivariate Gaussian distribution can enjoy five desirable properties. First, it is computationally tractable in the sense that it can be computed in a time, which is at most polynomial in dimension, sample size and the logarithm of the inverse of the contamination rate. Second, it is equivariant by translations, uniform scaling and orthogonal transformations. Third, it has a high breakdown point equal to 0.5, and a nearly minimax rate breakdown point approximately equal to 0.28. Fourth, it is minimax rate optimal, up to a logarithmic factor, when data consists of independent observations corrupted by adversarially chosen outliers. Fifth, it is asymptotically efficient when the rate of contamination tends to zero. The estimator is obtained by an iterative reweighting approach. Each sample point is assigned a weight that is iteratively updated by solving a convex optimization problem. We also establish a dimension-free nonasymptotic risk bound for the expected error of the proposed estimator. It is the first result of this kind in the literature and involves only the effective rank of the covariance matrix. Finally, we show that the obtained results can be extended to sub-Gaussian distributions, as well as to the cases of unknown rate of contamination or unknown covariance matrix.

In this paper, we consider the statistical inference for several low-rank tensor models. Specifically, in the Tucker low-rank tensor PCA or regression model, provided with any estimates achieving some attainable error rate, we develop the data-driven confidence regions for the singular subspace of the parameter tensor based on the asymptotic distribution of an updated estimate by two-iteration alternating minimization. The asymptotic distributions are established under some essential conditions on the signal-to-noise ratio (in PCA model) or sample size (in regression model). If the parameter tensor is further orthogonally decomposable, we develop the methods and nonasymptotic theory for inference on each individual singular vector. For the rank-one tensor PCA model, we establish the asymptotic distribution for general linear forms of principal components and confidence interval for each entry of the parameter tensor. Finally, numerical simulations are presented to corroborate our theoretical discoveries.

In all of these models, we observe that different from many matrix/vector settings in existing work, debiasing is not required to establish the asymptotic distribution of estimates or to make statistical inference on low-rank tensors. In fact, due to the widely observed statistical-computational-gap for low-rank tensor estimation, one usually requires stronger conditions than the statistical (or information-theoretic) limit to ensure the computationally feasible estimation is achievable. Surprisingly, such conditions “incidentally” render a feasible low-rank tensor inference without debiasing.