This paper describes a flexible framework for generalized low-rank tensor estimation problems that includes many important instances arising from applications in computational imaging, genomics, and network analysis. The proposed estimator consists of finding a low-rank tensor fit to the data under generalized parametric models. To overcome the difficulty of nonconvexity in these problems, we introduce a unified approach of projected gradient descent that adapts to the underlying low-rank structure. Under mild conditions on the loss function, we establish both an upper bound on statistical error and the linear rate of computational convergence through a general deterministic analysis. Then we further consider a suite of generalized tensor estimation problems, including sub-Gaussian tensor PCA, tensor regression, and Poisson and binomial tensor PCA. We prove that the proposed algorithm achieves the minimax optimal rate of convergence in estimation error. Finally, we demonstrate the superiority of the proposed framework via extensive experiments on both simulated and real data.

The statistical theory of extremes is extended to independent multivariate observations that are non-stationary both over time and across space. The non-stationarity over time and space is controlled via the scedasis (tail scale) in the marginal distributions. Spatial dependence stems from multivariate extreme value theory. We establish asymptotic theory for both the weighted sequential tail empirical process and the weighted tail quantile process based on all observations, taken over time and space. The results yield two statistical tests for homoscedasticity in the tail, one in space and one in time. Further, we show that the common extreme value index can be estimated via a pseudo-maximum likelihood procedure based on pooling all (non-stationary and dependent) observations. Our leading example and application is rainfall in Northern Germany.

A general framework for principal component analysis (PCA) in the presence of heteroskedastic noise is introduced. We propose an algorithm called HeteroPCA, which involves iteratively imputing the diagonal entries of the sample covariance matrix to remove estimation bias due to heteroskedasticity. This procedure is computationally efficient and provably optimal under the generalized spiked covariance model. A key technical step is a deterministic robust perturbation analysis on singular subspaces, which can be of independent interest. The effectiveness of the proposed algorithm is demonstrated in a suite of problems in high-dimensional statistics, including singular value decomposition (SVD) under heteroskedastic noise, Poisson PCA, and SVD for heteroskedastic and incomplete data.

We study predictive density estimation under Kullback–Leibler loss in ${\ell }_0$-sparse Gaussian sequence models. We propose proper Bayes predictive density estimates and establish asymptotic minimaxity in sparse models. Fundamental for this is a new risk decomposition for sparse, or spike-and-slab priors.

A surprise is the existence of a phase transition in the future-to-past variance ratio r. For $\mathit{r}<{\mathit{r}}_0=(\surd 5-1)/4$, the natural discrete prior ceases to be asymptotically optimal. Instead, for subcritical r, a ‘bi-grid’ prior with a central region of reduced grid spacing recovers asymptotic minimaxity. This phenomenon seems to have no analog in the otherwise parallel theory of point estimation of a multivariate normal mean under quadratic loss.

For spike-and-uniform slab priors to have any prospect of minimaxity, we show that the sparse parameter space needs also to be magnitude constrained. Within a substantial range of magnitudes, such spike-and-slab priors can attain asymptotic minimaxity.

We develop a general theory and estimation methods for functional linear sufficient dimension reduction, where both the predictor and the response can be random functions, or even vectors of functions. Unlike the existing dimension reduction methods, our approach does not rely on the estimation of conditional mean and conditional variance. Instead, it is based on a new statistical construction—the weak conditional expectation, which is based on Carleman operators and their inducing functions. Weak conditional expectation is a generalization of conditional expectation. Its key advantage is to replace the projection on to an ${\mathit{L}}_2$-space—which defines conditional expectation—by projection on to an arbitrary Hilbert space, while still maintaining the unbiasedness of the related dimension reduction methods. This flexibility is particularly important for functional data, because attempting to estimate a full-fledged conditional mean or conditional variance by slicing or smoothing over the space of vector-valued functions may be inefficient due to the curse of dimensionality. We evaluated the performances of the our new methods by simulation and in several applied settings.

We introduce the concept of pattern graphs–directed acyclic graphs representing how response patterns are associated. A pattern graph represents an identifying restriction that is nonparametrically identified/saturated and is often a missing not at random restriction. We introduce a selection model and a pattern mixture model formulations using the pattern graphs and show that they are equivalent. A pattern graph leads to an inverse probability weighting estimator as well as an imputation-based estimator. We also study the semiparametric efficiency theory and derive a multiply-robust estimator using pattern graphs.

Many complex networks in the real world can be formulated as hypergraphs where community detection has been widely used. However, the fundamental question of whether communities exist or not in an observed hypergraph remains unclear. This work aims to tackle this important problem. Specifically, we systematically study when a hypergraph with community structure can be successfully distinguished from its Erdős–Rényi counterpart, and propose concrete test statistics when the models are distinguishable. The main contribution of this paper is threefold. First, we discover a phase transition in the hyperedge probability for distinguishability. Second, in the bounded-degree regime, we derive a sharp signal-to-noise ratio (SNR) threshold for distinguishability in the special two-community 3-uniform hypergraphs, and derive nearly tight SNR thresholds in the general two-community m-uniform hypergraphs. Third, in the dense regime, we propose a computationally feasible test based on sub-hypergraph counts, obtain its asymptotic distribution, and analyze its power. Our results are further extended to nonuniform hypergraphs in which a new test involving both edge and hyperedge information is proposed. The proofs rely on Janson’s contiguity theory (Combin. Probab. Comput. 4 (1995) 369–405), a high-moments driven asymptotic normality result by Gao and Wormald (Probab. Theory Related Fields 130 (2004) 368–376), and a truncation technique for analyzing the likelihood ratio.

In high-dimensional regression, we attempt to estimate a parameter vector ${\mathit{\beta }}_0\in {\mathbb{R}}^{\mathit{p}}$ from $\mathit{n}\lesssim \mathit{p}$ observations ${\{({\mathit{y}}_{\mathit{i}},{\mathit{x}}_{\mathit{i}})\}}_{\mathit{i}\le \mathit{n}}$, where ${\mathit{x}}_{\mathit{i}}\in {\mathbb{R}}^{\mathit{p}}$ is a vector of predictors and ${\mathit{y}}_{\mathit{i}}$ is a response variable. A well-established approach uses convex regularizers to promote specific structures (e.g., sparsity) of the estimate $\hat{\mathit{\beta }}$ while allowing for practical algorithms. Theoretical analysis implies that convex penalization schemes have nearly optimal estimation properties in certain settings. However, in general the gaps between statistically optimal estimation (with unbounded computational resources) and convex methods are poorly understood.

We show that when the statistican has very simple structural information about the distribution of the entries of ${\mathit{\beta }}_0$, a large gap frequently exists between the best performance achieved by any convex regularizer satisfying a mild technical condition and either: (i) the optimal statistical error or (ii) the statistical error achieved by optimal approximate message passing algorithms. Remarkably, a gap occurs at high enough signal-to-noise ratio if and only if the distribution of the coordinates of ${\mathit{\beta }}_0$ is not log-concave. These conclusions follow from an analysis of standard Gaussian designs. Our lower bounds are expected to be generally tight, and we prove tightness under certain conditions.

Approximate Message Passing (AMP) algorithms have seen widespread use across a variety of applications. However, the precise forms for their Onsager corrections and state evolutions depend on properties of the underlying random matrix ensemble, limiting the extent to which AMP algorithms derived for white noise may be applicable to data matrices that arise in practice.

In this work, we study more general AMP algorithms for random matrices W that satisfy orthogonal rotational invariance in law, where W may have a spectral distribution that is different from the semicircle and Marcenko–Pastur laws characteristic of white noise. The Onsager corrections and state evolutions in these algorithms are defined by the free cumulants or rectangular free cumulants of the spectral distribution of W. Their forms were derived previously by Opper, Çakmak and Winther using nonrigorous dynamic functional theory techniques, and we provide rigorous proofs.

Our motivating application is a Bayes-AMP algorithm for Principal Components Analysis, when there is prior structure for the principal components (PCs) and possibly nonwhite noise. For sufficiently large signal strengths and any non-Gaussian prior distributions for the PCs, we show that this algorithm provably achieves higher estimation accuracy than the sample PCs.

Permutation tests are widely used in statistics, providing a finite-sample guarantee on the type I error rate whenever the distribution of the samples under the null hypothesis is invariant to some rearrangement. Despite its increasing popularity and empirical success, theoretical properties of the permutation test, especially its power, have not been fully explored beyond simple cases. In this paper, we attempt to partly fill this gap by presenting a general nonasymptotic framework for analyzing the minimax power of the permutation test. The utility of our proposed framework is illustrated in the context of two-sample and independence testing under both discrete and continuous settings. In each setting, we introduce permutation tests based on U-statistics and study their minimax performance. We also develop exponential concentration bounds for permuted U-statistics based on a novel coupling idea, which may be of independent interest. Building on these exponential bounds, we introduce permutation tests, which are adaptive to unknown smoothness parameters without losing much power. The proposed framework is further illustrated using more sophisticated test statistics including weighted U-statistics for multinomial testing and Gaussian kernel-based statistics for density testing. Finally, we provide some simulation results that further justify the permutation approach.

Model-X knockoffs (J. R. Stat. Soc. Ser. B. Stat. Methodol. 80 (2018) 551–577) allows analysts to perform feature selection using almost any machine learning algorithm while provably controlling the expected proportion of false discoveries. This procedure involves constructing synthetic variables, called knockoffs, which effectively act as controls during feature selection. The gold standard for constructing knockoffs has been to minimize the mean absolute correlation (MAC) between features and their knockoffs, but, surprisingly, we prove this procedure can be powerless in extremely easy settings, including Gaussian linear models with correlated exchangeable features. The key problem is that minimizing the MAC creates joint dependencies between the features and knockoffs, which allow machine learning algorithms to reconstruct the effect of the features on the response using the knockoffs. To improve power, we propose generating knockoffs which minimize the reconstructability (MRC) of the features, and we demonstrate our proposal for Gaussian features by showing it is computationally efficient, robust, and powerful. We also prove that certain MRC knockoffs minimize a notion of estimation error in Gaussian linear models. Through extensive simulations, we show MRC knockoffs often dramatically outperform MAC-minimizing knockoffs, and we find no settings in which MAC-minimizing knockoffs outperform MRC knockoffs by more than a slight margin. We implement our methods and many others from the knockoffs literature in a new python package $\mathtt{knockpy}$.

We consider least squares estimation in a general nonparametric regression model where the error is allowed to depend on the covariates. The rate of convergence of the least squares estimator (LSE) for the unknown regression function is well studied when the errors are sub-Gaussian. We find upper bounds on the rates of convergence of the LSE when the error has a uniformly bounded conditional variance and has only finitely many moments. Our upper bound on the rate of convergence of the LSE depends on the moment assumptions on the error, the metric entropy of the class of functions involved and the “local” structure of the function class around the truth. We find sufficient conditions on the error distribution under which the rate of the LSE matches the rate of the LSE under sub-Gaussian error. Our results are finite sample and allow for heteroscedastic and heavy-tailed errors.

This paper considers the deconvolution problem in the case where the target signal is multidimensional and no information is known about the noise distribution. More precisely, no assumption is made on the noise distribution and no samples are available to estimate it: the deconvolution problem is solved based only on observations of the corrupted signal. We establish the identifiability of the model up to translation when the signal has a Laplace transform with an exponential growth ρ smaller than 2 and when it can be decomposed into two dependent components. Then we propose an estimator of the probability density function of the signal, which is consistent for any unknown noise distribution with finite variance. We also prove rates of convergence and, as the estimator depends on ρ which is usually unknown, we propose a model selection procedure to obtain an adaptive estimator with the same rate of convergence as the estimator with a known tail parameter. This rate of convergence is known to be minimax when $\mathit{\rho }=1$.

Motivated by models for multiway comparison data, we consider the problem of estimating a coordinatewise isotonic function on the domain ${[0,1]}^{\mathit{d}}$ from noisy observations collected on a uniform lattice, but where the design points have been permuted along each dimension. While the univariate and bivariate versions of this problem have received significant attention, our focus is on the multivariate case $\mathit{d}\ge 3$. We study both the minimax risk of estimation (in empirical ${\mathit{L}}_2$ loss) and the fundamental limits of adaptation (quantified by the adaptivity index) to a family of piecewise constant functions. We provide a computationally efficient Mirsky partition estimator that is minimax optimal while also achieving the smallest adaptivity index possible for polynomial time procedures. Thus, from a worst-case perspective and in sharp contrast to the bivariate case, the latent permutations in the model do not introduce significant computational difficulties over and above vanilla isotonic regression. On the other hand, the fundamental limits of adaptation are significantly different with and without unknown permutations: Assuming a hardness conjecture from average-case complexity theory, a statistical-computational gap manifests in the former case. In a complementary direction, we show that natural modifications of existing estimators fail to satisfy at least one of the desiderata of optimal worst-case statistical performance, computational efficiency and fast adaptation. Along the way to showing our results, we improve adaptation results in the special case $\mathit{d}=2$ and establish some properties of estimators for vanilla isotonic regression, both of which may be of independent interest.

Methods of merging several p-values into a single p-value are important in their own right and widely used in multiple hypothesis testing. This paper is the first to systematically study the admissibility (in Wald’s sense) of p-merging functions and their domination structure, without any information on the dependence structure of the input p-values. As a technical tool, we use the notion of e-values, which are alternatives to p-values recently promoted by several authors. We obtain several results on the representation of admissible p-merging functions via e-values and on (in)admissibility of existing p-merging functions. By introducing new admissible p-merging functions, we show that some classic merging methods can be strictly improved to enhance power without compromising validity under arbitrary dependence.

In the Gaussian sequence model $\mathit{Y}=\mathit{\mu }+\mathit{\xi }$, we study the likelihood ratio test (LRT) for testing ${\mathit{H}}_0:\mathit{\mu }={\mathit{\mu }}_0$ versus ${\mathit{H}}_1:\mathit{\mu }\in \mathit{K}$, where ${\mathit{\mu }}_0\in \mathit{K}$, and K is a closed convex set in ${\mathbb{R}}^{\mathit{n}}$. In particular, we show that under the null hypothesis, normal approximation holds for the log-likelihood ratio statistic for a general pair $({\mathit{\mu }}_0,\mathit{K})$, in the high-dimensional regime where the estimation error of the associated least squares estimator diverges in an appropriate sense. The normal approximation further leads to a precise characterization of the power behavior of the LRT in the high-dimensional regime. These characterizations show that the power behavior of the LRT is in general nonuniform with respect to the Euclidean metric, and illustrate the conservative nature of existing minimax optimality and suboptimality results for the LRT. A variety of examples, including testing in the orthant/circular cone, isotonic regression, Lasso and testing parametric assumptions versus shape-constrained alternatives, are worked out to demonstrate the versatility of the developed theory.

In this paper, we test whether two data sets measured on the same set of subjects share a common clustering structure. As a leading example, we focus on comparing clustering structures in two independent random samples from two deterministic two-component mixtures of multivariate Gaussian distributions. Mean parameters of these Gaussian distributions are treated as potentially unknown nuisance parameters and are allowed to differ. Assuming knowledge of mean parameters, we first determine the phase diagram of the testing problem over the entire range of signal-to-noise ratios by providing both lower bounds and tests that achieve them. When nuisance parameters are unknown, we propose tests that achieve the detection boundary adaptively as long as ambient dimensions of the data sets grow at a sublinear rate with the sample size.

In classical statistics, Pinsker’s theorem provides an exact asymptotic minimax bound in nonparametric estimation, improving upon optimal rates of convergence results. We obtain a quantum version of the theorem by establishing asymptotic minimax results for estimation of the displacement vector in a quantum Gaussian white noise model, given by a sequence of shifted vacuum states. Analogous results are then obtained for estimation of a general pure state from an ensemble of identically prepared, independent quantum systems, using the recently established local asymptotic equivalence to the quantum Gaussian white noise model. Optimality holds with respect to the most fundamental distance measure for quantum states, that is, trace norm distance, and in a true quantum sense, allowing for all possible measurements. Adaptive estimators are also obtained for the above cases. As an application, we obtain asymptotic minimax adaptive estimators for Wigner functions of pure states.

We consider a high-dimensional linear regression problem. Unlike many papers on the topic, we do not require sparsity of the regression coefficients; instead, our main structural assumption is a decay of eigenvalues of the covariance matrix of the data. We propose a new family of estimators, called the canonical thresholding estimators, which pick largest regression coefficients in the canonical form. The estimators admit an explicit form and can be linked to LASSO and Principal Component Regression (PCR). A theoretical analysis for both fixed design and random design settings is provided. Obtained bounds on the mean squared error and the prediction error of a specific estimator from the family allow to clearly state sufficient conditions on the decay of eigenvalues to ensure convergence. In addition, we promote the use of the relative errors, strongly linked with the out-of-sample ${\mathit{R}}^2$. The study of these relative errors leads to a new concept of joint effective dimension, which incorporates the covariance of the data and the regression coefficients simultaneously, and describes the complexity of a linear regression problem. Some minimax lower bounds are established to showcase the optimality of our procedure. Numerical simulations confirm good performance of the proposed estimators compared to the previously developed methods.

In long-term follow-up studies, data are often collected on repeated measures of multivariate response variables as well as on time to the occurrence of a certain event. To jointly analyze such longitudinal data and survival time, we propose a general class of semiparametric latent-class models that accommodates a heterogeneous study population with flexible dependence structures between the longitudinal and survival outcomes. We combine nonparametric maximum likelihood estimation with sieve estimation and devise an efficient EM algorithm to implement the proposed approach. We establish the asymptotic properties of the proposed estimators through novel use of modern empirical process theory, sieve estimation theory and semiparametric efficiency theory. Finally, we demonstrate the advantages of the proposed methods through extensive simulation studies and provide an application to the Atherosclerosis Risk in Communities study.

We construct an algorithm for estimating the mean of a heavy-tailed random variable when given an adversarial corrupted sample of N independent observations. The only assumption we make on the distribution of the noncorrupted (or informative) data is the existence of a covariance matrix Σ, unknown to the statistician. Our algorithm outputs $\hat{\mathit{\mu }}$, which is robust to the presence of $|\mathcal{O}|$ adversarial outliers and satisfies

$$(1)\Vert \hat{\mathit{\mu }}-\mathit{\mu }{\Vert }_2\lesssim \sqrt{\frac{Tr(\mathrm{\Sigma })}{\mathit{N}}}+\sqrt{\frac{\Vert \mathrm{\Sigma }{\Vert }_{\mathrm{op}}\mathit{K}}{\mathit{N}}}$$

with probability at least $1-exp(-{\mathit{c}}_0\mathit{K})-exp(-{\mathit{c}}_1\mathit{u})$, and runtime $\tilde{\mathcal{O}}(\mathit{N}\mathit{d}+\mathit{u}\mathit{K}\mathit{d})$ where $\mathit{K}\in \{600|\mathcal{O}|,\dots ,\mathit{N}\}$ and $\mathit{u}\in {\mathbb{N}}^{\ast }$ are two parameters of the algorithm. The algorithm is fully data-dependent and does not use (1) in its construction, which combines recently developed tools for median-of-means estimators and covering semidefinite programming. We also show that this algorithm can automatically adapt to the number of outliers (adaptive choice of K) and that it satisfies the same bound in expectation.

We consider the problem of estimating the distribution of time-to-event data that is subject to censoring and for which the event of interest might never occur, that is, some subjects are cured. To model this kind of data in the presence of covariates, one of the leading semiparametric models is the promotion time cure model (Stochastic Models of Tumor Latency and Their Biostatistical Applications (1996) World Scientific), which adapts the Cox model to the presence of cured subjects. Estimating the conditional distribution results in a complicated constrained optimization problem, and inference is difficult as no closed-formula for the variance is available. We propose a new model, inspired by the Cox model, that leads to a simple estimation procedure and that presents a closed formula for the variance. In this paper, we show (i) that the new model contains as a special case the promotion time cure model with an exponential link, and hence we have a simpler way to estimate the latter model than what is done so far in the literature; (ii) that in the latter special case, both estimators are equal to the partial likelihood estimator under the usual Cox model; (iii) that the estimators under the new model have certain asymptotic properties when the model is correct and when it is misspecified; (iv) that the error of LASSO type estimators is of order $\sqrt{log(\mathit{n}\mathit{d})/\mathit{n}}$ in the case of high-dimensional covariates with dimension d and sample size n. We also study the practical behaviour of our estimation procedure by means of simulations, and we apply our model and estimation method to a breast cancer data set.

Regression models with crossed random effect errors can be very expensive to compute. The cost of both generalized least squares and Gibbs sampling can easily grow as ${\mathit{N}}^{3/2}$ (or worse) for N observations. Papaspiliopoulos, Roberts and Zanella (Biometrika 107 (2020) 25–40) present a collapsed Gibbs sampler that costs $\mathit{O}(\mathit{N})$, but under an extremely stringent sampling model. We propose a backfitting algorithm to compute a generalized least squares estimate and prove that it costs $\mathit{O}(\mathit{N})$. A critical part of the proof is in ensuring that the number of iterations required is $\mathit{O}(1)$, which follows from keeping a certain matrix norm below $1-\mathit{\delta }$ for some $\mathit{\delta }>0$. Our conditions are greatly relaxed compared to those for the collapsed Gibbs sampler, though still strict. Empirically, the backfitting algorithm has a norm below $1-\mathit{\delta }$ under conditions that are less strict than those in our assumptions. We illustrate the new algorithm on a ratings data set from Stitch Fix.

This paper studies the statistical and computational limits of high-order clustering with planted structures. We focus on two clustering models, constant high-order clustering (CHC) and rank-one higher-order clustering (ROHC), and study the methods and theory for testing whether a cluster exists (detection) and identifying the support of cluster (recovery).

Specifically, we identify the sharp boundaries of signal-to-noise ratio for which CHC and ROHC detection/recovery are statistically possible. We also develop the tight computational thresholds: when the signal-to-noise ratio is below these thresholds, we prove that polynomial-time algorithms cannot solve these problems under the computational hardness conjectures of hypergraphic planted clique (HPC) detection and hypergraphic planted dense subgraph (HPDS) recovery. We also propose polynomial-time tensor algorithms that achieve reliable detection and recovery when the signal-to-noise ratio is above these thresholds. Both sparsity and tensor structures yield the computational barriers in high-order tensor clustering. The interplay between them results in significant differences between high-order tensor clustering and matrix clustering in literature in aspects of statistical and computational phase transition diagrams, algorithmic approaches, hardness conjecture, and proof techniques. To our best knowledge, we are the first to give a thorough characterization of the statistical and computational trade-off for such a double computational-barrier problem. Finally, we provide evidence for the computational hardness conjectures of HPC detection (via low-degree polynomial and Metropolis methods) and HPDS recovery (via low-degree polynomial method).