The workshop highlighted recent theoretical advances on inference in high-dimensional statistical models based on the interplay of techniques from mathematical statistics, machine learning, theoretical computer science and related areas. The workshop brought together about 50 researchers in order to present new results, exchange ideas and explore open problems.

Motivated by the statistical analysis of the discrete optimal transport problem, we prove distributional limits for the solutions of linear programs with random constraints. Such limits were first obtained by Klatt, Munk, & Zemel (2022), but their expressions for the limits involve a computationally intractable decomposition of $\mathbb{R}^m$ into a possibly exponential number of convex cones. We give a new expression for the limit in terms of auxiliary linear programs, which can be solved in polynomial time. We also leverage tools from random convex geometry to give distributional limits for the entire set of random optimal solutions, when the optimum is not unique. Finally, we describe a simple, data-driven method to construct asymptotically valid confidence sets in polynomial time.

High-dimensional cellular and molecular profiling of human samples highlights the need for analytical approaches that can integrate multi-omic datasets to generate predictive biomarkers that are in turn accompanied with strong causal inferences. Current methodologies are challenged by the high dimensionality of the combined datasets, the differences in distributions across the datasets, and their integration in a plausible causal framework, beyond merely correlative biomarkers. Here we present CausER, a first-in-class two-step interpretable machine learning approach for high-dimensional multi-omic datasets, that addresses these problems by identifying latent factors and their cause-effect relationships with the system-wide outcome/property of interest. The first step consists of Essential Regression (ER), a novel data-distribution-free regression model that integrates multi-omic datasets and identifies latent factors significantly associated with an outcome. The second involves probabilistic graphical modeling of the significant latent factors to infer plausible causal associations between them and mechanisms that affect outcomes, thereby significantly moving beyond predictive associative markers. By analyzing varied human immunological multi-omic datasets, we demonstrate that CausER significantly outperforms a wide range of state-of-the-art approaches. It generates novel cellular and molecular predictions in a range of contexts, including immunosenescence and sustained immune dysregulation associated with pre-term birth, that are corroborated by biological findings in model systems.

Motivated by modern applications in which one constructs graphical models based on a very large number of features, this paper introduces a new class of cluster-based graphical models, in which variable clustering is applied as an initial step for reducing the dimension of the feature space. We employ model assisted clustering, in which the clusters contain features that are similar to the same unobserved latent variable. Two different cluster-based Gaussian graphical models are considered: the latent variable graph, corresponding to the graphical model associated with the unobserved latent variables, and the cluster-average graph, corresponding to the vector of features averaged over clusters. Our study reveals that likelihood based inference for the latent graph, not analyzed previously, is analytically intractable. Our main contribution is the development and analysis of alternative estimation and inference strategies, for the precision matrix of an unobservable latent vector $Z$. We replace the likelihood of the data by an appropriate class of empirical risk functions, that can be specialized to the latent graphical model and to the simpler, but under-analyzed, cluster-average graphical model. The estimators thus derived can be used for inference on the graph structure, for instance on edge strength or pattern recovery. Inference is based on the asymptotic limits of the entry-wise estimates of the precision matrices associated with the conditional independence graphs under consideration. While taking the uncertainty induced by the clustering step into account, we establish Berry-Esseen central limit theorems for the proposed estimators. It is noteworthy that, although the clusters are estimated adaptively from the data, the central limit theorems regarding the entries of the estimated graphs are proved under the same conditions one would use if the clusters were known....

The study of complex relationships among the elements of a large collection of random variables lead to the development of a number of areas in probability and statistics such as probabilistic network analysis or random matrix theory. The aim of the workshop was to address the challenge to develop a coherent mathematical framework within which these areas can be integrated, for a successful analysis of massive and complicated data sets.

Variable clustering is one of the most important unsupervised learning methods, ubiquitous in most research areas. In the statistics and computer science literature, most of the clustering methods lead to non-overlapping partitions of the variables. However, in many applications, some variables may belong to multiple groups, yielding clusters with overlap. It is still largely unknown how to perform overlapping variable clustering with statistical guarantees. To bridge this gap, we propose a novel Latent model-based OVErlapping clustering method (LOVE) to recover overlapping sub-groups of a potentially very large group of variables. In our model-based formulation, a cluster is given by variables associated with the same latent factor, and can be determined from an allocation matrix A that indexes our proposed latent model. We assume that some of the observed variables are pure, in that they are associated with only one latent factor, whereas the remaining majority has multiple allocations. We prove that the corresponding allocation matrix A, and the induced overlapping clusters, are identifiable, up to label switching. We estimate the clusters with LOVE, our newly developed algorithm, which consists in two steps. The first step estimates the set of pure variables, and the number of clusters. In the second step we estimate the allocation matrix A and determine the overlapping clusters. Under minimal signal strength conditions, our algorithm recovers the population level clusters consistently. Our theoretical results are fully supported by our empirical studies, which include extensive simulation studies that compare LOVE with other existing methods, and the analysis of a RNA-seq dataset.

We introduce a new sparse estimator of the covariance matrix for high-dimensional models in which the variables have a known ordering. Our estimator, which is the solution to a convex optimization problem, is equivalently expressed as an estimator which tapers the sample covariance matrix by a Toeplitz, sparsely-banded, data-adaptive matrix. As a result of this adaptivity, the convex banding estimator enjoys theoretical optimality properties not attained by previous banding or tapered estimators. In particular, our convex banding estimator is minimax rate adaptive in Frobenius and operator norms, up to log factors, over commonly-studied classes of covariance matrices, and over more general classes. Furthermore, it correctly recovers the bandwidth when the true covariance is exactly banded. Our convex formulation admits a simple and efficient algorithm. Empirical studies demonstrate its practical effectiveness and illustrate that our exactly-banded estimator works well even when the true covariance matrix is only close to a banded matrix, confirming our theoretical results. Our method compares favorably with all existing methods, in terms of accuracy and speed. We illustrate the practical merits of the convex banding estimator by showing that it can be used to improve the performance of discriminant analysis for classifying sound recordings.

We introduce a new criterion, the Rank Selection Criterion (RSC), for selecting the optimal reduced rank estimator of the coefficient matrix in multivariate response regression models. The corresponding RSC estimator minimizes the Frobenius norm of the fit plus a regularization term proportional to the number of parameters in the reduced rank model.The rank of the RSC estimator provides a consistent estimator of the rank of the coefficient matrix; in general, the rank of our estimator is a consistent estimate of the effective rank, which we define to be the number of singular values of the target matrix that are appropriately large. The consistency results are valid not only in the classic asymptotic regime, when n, the number of responses, and p, the number of predictors, stay bounded, and m, the number of observations, grows, but also when either, or both, n and p grow, possibly much faster than m.We establish minimax optimal bounds on the mean squared errors of our estimators. Our finite sample performance bounds for the RSC estimator show that it achieves the optimal balance between the approximation error and the penalty term.Furthermore, our procedure has very low computational complexity, linear in the number of candidate models, making it particularly appealing for large scale problems. We contrast our estimator with the nuclear norm penalized least squares (NNP) estimator, which has an inherently higher computational complexity than RSC, for multivariate regression models. We show that NNP has estimation properties similar to those of RSC, albeit under stronger conditions. However, it is not as parsimonious as RSC.We offer a simple correction of the NNP estimator which leads to consistent rank estimation. We verify and illustrate our theoretical findings via an extensive simulation study.

We propose and analyse fully data-driven methods for inference about the mean function of a Gaussian process from a sample of independent trajectories of the process, observed at random time points and corrupted by additive random error. Our methods are based on thresholded least squares estimators relative to an approximating function basis. The variable threshold levels are determined from the data and the resulting estimates adapt to the unknown sparsity of the mean function relative to the approximating basis. These results are obtained via novel oracle inequalities, which are further used to derive the rates of convergence of our mean estimates. In addition, we construct confidence balls that adapt to the unknown regularity of the mean and covariance function of the stochastic process. They are easy to compute since they do not require explicit estimation of the covariance operator of the process. A simulation study shows that the new method performs very well in practice and is robust against large variations that may be introduced by the random-error terms.

The goals of this paper are to review the most popular methods of predictor selection in regression models, to explain why some fail when the number P of explanatory variables exceeds the number N of participants, and to discuss alternative statistical methods that can be employed in this case. We focus on penalized least squares methods in regression models, and discuss in detail two such methods that are well established in the statistical literature, the LASSO and Elastic Net. We introduce bootstrap enhancements of these methods, the BE-LASSO and BE-Enet, that allow the user to attach a measure of uncertainty to each variable selected. Our work is motivated by a multimodal neuroimaging dataset that consists of morphometric measures (volumes at several anatomical regions of interest), white matter integrity measures from diffusion weighted data (fractional anisotropy, mean diffusivity, axial diffusivity and radial diffusivity) and clinical and demographic variables (age, education, alcohol and drug history). In this dataset, the number P of explanatory variables exceeds the number N of participants. We use the BE-LASSO and BE-Enet to provide the first statistical analysis that allows the assessment of neurocognitive performance from high dimensional neuroimaging and clinical predictors, including their interactions. The major novelty of this analysis is that biomarker selection and dimension reduction are accomplished with a view towards obtaining good predictions for the outcome of interest (i.e., the neurocognitive indices), unlike principal component analysis that are performed only on the predictors' space independently of the outcome of interest.

This paper studies sparse density estimation via l1 penalization (SPADES). We focus on estimation in high-dimensional mixture models and nonparametric adaptive den- sity estimation. We show, respectively, that SPADES can recover, with high probability, the unknown components of a mixture of probability densities and that it yields minimax adaptive density estimates. These results are based on a general sparsity oracle inequality that the SPADES estimates satisfy. MSC2000 Subject classification: Primary 62G08, Secondary 62C20, 62G05, 62G20

This paper proposes and analyzes fully data driven methods for inference about the mean function of a stochastic process from a sample of independent trajectories of the process, observed at discrete time points and corrupted by additive random error. The proposed method uses thresholded least squares estimators relative to an approximating function basis. The variable threshold levels are estimated from the data and the basis is chosen via cross-validation from a library of bases. The resulting estimates adapt to the unknown sparsity of the mean function relative to the selected approximating basis, both in terms of the mean squared error and supremum norm. These results are based on novel oracle inequalities. In addition, uniform confidence bands for the mean function of the process are constructed. The bands also adapt to the unknown regularity of the mean function, are easy to compute, and do not require explicit estimation of the covariance operator of the process. The simulation study that complements the theoretical results shows that the new method performs very well in practice, and is robust against large variations introduced by the random error terms.

This paper investigates correct variable selection in finite samples via l(1) and l(1) + l(2) type penalization schemes. The asymptotic consistency of variable selection immediately follows from this analysis. We focus on logistic and linear regression models. The following questions are central to our paper: given a level of confidence 1 - delta, under which assumptions on the design matrix, for which strength of the signal and for what values of the tuning parameters can we identify the true model at the given level of confidence? Formally, if (I) over cap is an estimate of the true variable set I*, we study conditions under which P((I) over cap = I*) >= 1 - delta, for a given sample size n, number of parameters 34 and confidence 1 - delta. We show that in identifiable models, both methods can recover coefficients of size 1/root n, up to small multiplicative constants and logarithmic factors in M and 1/delta. The advantage of the l(1) + l(2) penalization over the l(1) is minor for the variable selection problem, for the models we consider here. Whereas the former estimates are unique, and become more stable for highly correlated data matrices as one increases the tuning parameter of the l(2) part, too large an increase in this parameter value may preclude variable selection.

In this article we investigate consistency of selection in regression models via the popular Lasso method. Here we depart from the traditional linear regression assumption and consider approximations of the regression function f with elements of a given dictionary of M functions. The target for consistency is the index set of those functions from this dictionary that realize the most parsimonious approximation to f among all linear combinations belonging to an L_2 ball centered at f and of radius r_n,M^2. In this framework we show that a consistent estimate of this index set can be derived via ℓ_1 penalized least squares, with a data dependent penalty and with tuning sequence r_n,M>√(log(Mn)/n), where n is the sample size. Our results hold for any 1≤ M≤ n^γ, for any γ>0.

In this correspondence, a sequential procedure for aggregating linear combinations of a finite family of regression estimates is described and analyzed. Particular attention is given to linear combinations having coefficients in the generalized simplex. The procedure is based on exponential weighting, and has a computationally tractable approximation. Analysis of the procedure is based in part on techniques from the sequential prediction of nonrandom sequences. Here these techniques are applied in a stochastic setting to obtain cumulative loss bounds for the aggregation procedure. From the cumulative loss bounds we derive an oracle inequality for the aggregate estimator for an unbounded response having a suitable moment-generating function. The inequality shows that the risk of the aggregate estimator is less than the risk of the best candidate linear combination in the generalized simplex, plus a complexity term that depends on the size of the coefficient set. The inequality readily yields convergence rates for aggregation over the unit simplex that are within logarithmic factors of known minimax bounds. Some preliminary results on model selection are also presented.

This paper studies statistical aggregation procedures in the regression setting. A motivating factor is the existence of many different methods of estimation, leading to possibly competing estimators. We consider here three different types of aggregation: model selection (MS) aggregation, convex (C) aggregation and linear (L) aggregation. The objective of (MS) is to select the optimal single estimator from the list; that of (C) is to select the optimal convex combination of the given estimators; and that of (L) is to select the optimal linear combination of the given estimators. we are interested in evaluating the rates of convergence of the excess risks of the estimators obtained by these procedures. Our approach is motivated by recently published minimax results [Nemirovski, A. (2000). Topics in non-parametric statistics. Lectures on Probability Theory and Statistics (Saint-Flour, 1998). Lecture Notes in Math. 1738 85-277. Springer, Berlin; Tsybakov, A. B. (2003). Optimal rates of aggregation. Learning Theory and Kernel Machines. Lecture Notes in Artificial Intelligence 2777 303-313. Springer, Heidelberg]. There exist competing aggregation procedures achieving optimal convergence rates for each of the (MS), (C) and (L) cases separately. Since these procedures are not directly comparable with each other, we suggest an alternative solution. We prove that all three optimal rates, as well as those for the newly introduced (S) aggregation (subset selection), are nearly achieved via a single "universal" aggregation procedure. The procedure consists of mixing the initial estimators with weights obtained by penalized least squares. Two different penalties are considered: one of them is of the BIC type, the second one is a data-dependent l(1)-type penalty.

This paper studies oracle properties of ℓ_1-penalized least squares in nonparametric regression setting with random design. We show that the penalized least squares estimator satisfies sparsity oracle inequalities, i.e., bounds in terms of the number of non-zero components of the oracle vector. The results are valid even when the dimension of the model is (much) larger than the sample size and the regression matrix is not positive definite. They can be applied to high-dimensional linear regression, to nonparametric adaptive regression estimation and to the problem of aggregation of arbitrary estimators.

We develop a statistical method for estimating the spectrum from a data set that consists of several signals, all of which are realizations of a common random process. We first find estimates of the common spectrum using each signal; then we construct $M$ partial aggregates. Each partial aggregate is a linear combination of $M-$ 1 of the spectral estimates. The weights are obtained from the data via a least squares criterion. The final spectral estimate is the average of these $M$ partial aggregates. We show that our final estimator is minimax rate adaptive if at least two of the estimators per signal attain the optimal rate $n^-2alpha /2alpha + 1$ for spectra belonging to a generalized Lipschitz ball with smoothness index $alpha $ . Our simulation study strongly suggests that our procedure works well in practice, and in a large variety of situations is preferable to the simple averaging of the $M$ spectral estimates.

This paper connects consistent variable selection with multiple hypotheses testing procedures in the linear regression model Y=Xβ+ɛ, where the dimension p of the parameter β is allowed to grow with the sample size n. We view the variable selection problem as one of estimating the index set I0⊆{1,…,p} of the non-zero components of β∈Rp. Estimation of I0 can be further reformulated in terms of testing the hypotheses β1=0,…,βp=0. We study here testing via the false discovery rate (FDR) and Bonferroni methods. We show that the set I^⊆{1,…,p} consisting of the indices of rejected hypotheses βi=0 is a consistent estimator of I0, under appropriate conditions on the design matrix X and the control values used in either procedure. This technique can handle situations where p is large at a very low computational cost, as no exhaustive search over the space of the 2p submodels is required.

Let X be a random variable taking values in a separable Hilbert space X, with label Y/spl isin/{0,1}. We establish universal weak consistency of a nearest neighbor-type classifier based on n independent copies (X/sub i/,Y/sub i/) of the pair (X,Y), extending the classical result of Stone to infinite-dimensional Hilbert spaces. Under a mild condition on the distribution of X, we also prove strong consistency. We reduce the infinite dimension of X by considering only the first d coefficients of a Fourier series expansion of each X/sub i/, and then we perform k-nearest neighbor classification in /spl Ropf//sup d/. Both the dimension and the number of neighbors are automatically selected from the data using a simple data-splitting device. An application of this technique to a signal discrimination problem involving speech recordings is presented.