Search

Random projections or sketching are widely used in many algorithmic and learning contexts. Here we study the performance of iterative Hessian sketch for least-squares problems. By leveraging and extending recent results from random matrix theory on the limiting spectrum of matrices randomly projected with the subsampled randomized Hadamard transform, and truncated Haar matrices, we can study and compare the resulting algorithms to a level of precision that has not been possible before. Our technical contributions include a novel formula for the second moment of the inverse of projected matrices. We also find simple closed-form expressions for asymptotically optimal step-sizes and convergence rates. These show that the convergence rate for Haar and randomized Hadamard matrices are identical, and asymptotically improve upon Gaussian random projections. These techniques may be applied to other algorithms that employ randomized dimension reduction.

Normalization methods such as batch [Ioffe and Szegedy, 2015], weight [Salimansand Kingma, 2016], instance [Ulyanov et al., 2016], and layer normalization [Baet al., 2016] have been widely used in modern machine learning. Here, we study the weight normalization (WN) method [Salimans and Kingma, 2016] and a variant called reparametrized projected gradient descent (rPGD) for overparametrized least-squares regression. WN and rPGD reparametrize the weights with a scale g and a unit vector w and thus the objective function becomes non-convex. We show that this non-convex formulation has beneficial regularization effects compared to gradient descent on the original objective. These methods adaptively regularize the weights and converge close to the minimum l2 norm solution, even for initializations far from zero. For certain stepsizes of g and w , we show that they can converge close to the minimum norm solution. This is different from the behavior of gradient descent, which converges to the minimum norm solution only when started at a point in the range space of the feature matrix, and is thus more sensitive to initialization., Comment: NeurIPS 2020

We study the following three fundamental problems about ridge regression: (1) what is the structure of the estimator? (2) how to correctly use cross-validation to choose the regularization parameter? and (3) how to accelerate computation without losing too much accuracy? We consider the three problems in a unified large-data linear model. We give a precise representation of ridge regression as a covariance matrix-dependent linear combination of the true parameter and the noise. We study the bias of $K$-fold cross-validation for choosing the regularization parameter, and propose a simple bias-correction. We analyze the accuracy of primal and dual sketching for ridge regression, showing they are surprisingly accurate. Our results are illustrated by simulations and by analyzing empirical data., Comment: Published as a conference paper at ICLR 2020

Data augmentation is a widely used trick when training deep neural networks: in addition to the original data, properly transformed data are also added to the training set. However, to the best of our knowledge, a clear mathematical framework to explain the performance benefits of data augmentation is not available. In this paper, we develop such a theoretical framework. We show data augmentation is equivalent to an averaging operation over the orbits of a certain group that keeps the data distribution approximately invariant. We prove that it leads to variance reduction. We study empirical risk minimization, and the examples of exponential families, linear regression, and certain two-layer neural networks. We also discuss how data augmentation could be used in problems with symmetry where other approaches are prevalent, such as in cryo-electron microscopy (cryo-EM)., Comment: To appear in Journal of Machine Learning Research

In many areas, practitioners need to analyze large datasets that challenge conventional single-machine computing. To scale up data analysis, distributed and parallel computing approaches are increasingly needed. Here we study a fundamental and highly important problem in this area: How to do ridge regression in a distributed computing environment? Ridge regression is an extremely popular method for supervised learning, and has several optimality properties, thus it is important to study. We study one-shot methods that construct weighted combinations of ridge regression estimators computed on each machine. By analyzing the mean squared error in a high dimensional random-effects model where each predictor has a small effect, we discover several new phenomena. 1. Infinite-worker limit: The distributed estimator works well for very large numbers of machines, a phenomenon we call "infinite-worker limit". 2. Optimal weights: The optimal weights for combining local estimators sum to more than unity, due to the downward bias of ridge. Thus, all averaging methods are suboptimal. We also propose a new Weighted ONe-shot DistributEd Ridge regression (WONDER) algorithm. We test WONDER in simulation studies and using the Million Song Dataset as an example. There it can save at least 100x in computation time, while nearly preserving test accuracy., Comment: Gave the name "Wonder" to the algorithm, updated title, added algorithm for general non-isotropic design

We consider a least squares regression problem where the data has been generated from a linear model, and we are interested to learn the unknown regression parameters. We consider "sketch-and-solve" methods that randomly project the data first, and do regression after. Previous works have analyzed the statistical and computational performance of such methods. However, the existing analysis is not fine-grained enough to show the fundamental differences between various methods, such as the Subsampled Randomized Hadamard Transform (SRHT) and Gaussian projections. In this paper, we make progress on this problem, working in an asymptotic framework where the number of datapoints and dimension of features goes to infinity. We find the limits of the accuracy loss (for estimation and test error) incurred by popular sketching methods. We show separation between different methods, so that SRHT is better than Gaussian projections. Our theoretical results are verified on both real and synthetic data. The analysis of SRHT relies on novel methods from random matrix theory that may be of independent interest.

Distributed statistical learning problems arise commonly when dealing with large datasets. In this setup, datasets are partitioned over machines, which compute locally, and communicate short messages. Communication is often the bottleneck. In this paper, we study one-step and iterative weighted parameter averaging in statistical linear models under data parallelism. We do linear regression on each machine, send the results to a central server, and take a weighted average of the parameters. Optionally, we iterate, sending back the weighted average and doing local ridge regressions centered at it. How does this work compared to doing linear regression on the full data? Here we study the performance loss in estimation, test error, and confidence interval length in high dimensions, where the number of parameters is comparable to the training data size. We find the performance loss in one-step weighted averaging, and also give results for iterative averaging. We also find that different problems are affected differently by the distributed framework. Estimation error and confidence interval length increase a lot, while prediction error increases much less. We rely on recent results from random matrix theory, where we develop a new calculus of deterministic equivalents as a tool of broader interest., Comment: Fixing a typo

Drawing statistical inferences from large datasets in a model-robust way is an important problem in statistics and data science. In this paper, we propose methods that are robust to large and unequal noise in different observational units (i.e., heteroskedasticity) for statistical inference in linear regression. We leverage the Hadamard estimator, which is unbiased for the variances of ordinary least-squares regression. This is in contrast to the popular White's sandwich estimator, which can be substantially biased in high dimensions. We propose to estimate the signal strength, noise level, signal-to-noise ratio, and mean squared error via the Hadamard estimator. We develop a new degrees of freedom adjustment that gives more accurate confidence intervals than variants of White's sandwich estimator. Moreover, we provide conditions ensuring the estimator is well-defined, by studying a new random matrix ensemble in which the entries of a random orthogonal projection matrix are squared. We also show approximate normality, using the second-order Poincare inequality. Our work provides improved statistical theory and methods for linear regression in high dimensions.

Multiple hypothesis testing problems arise naturally in science. In this paper, we introduce the new Fast Closed Testing (FACT) method for multiple testing, controlling the family-wise error rate. This error rate is state of the art in many important application areas, and is preferred to false discovery rate control for many reasons, including that it leads to stronger reproducibility. The closure principle rejects an individual hypothesis if all global nulls of subsets containing it are rejected using some test statistics. It takes exponential time in the worst case. When the tests are symmetric and monotone, our method is an exact algorithm for computing the closure, quadratic in the number of tests, and linear in the number of discoveries. Our framework generalizes most examples of closed testing such as Holm's and the Bonferroni method. As a special case of our method, we propose the Simes-higher criticism fusion test, which is powerful for detecting both a few strong signals, and also many moderate signals., Comment: This version fixes an indexing error in the calculation of adjusted p-values, pointed out in arXiv:1912.06116. To appear in Biometrika

Factor analysis and principal component analysis (PCA) are used in many application areas. The first step, choosing the number of components, remains a serious challenge. Our work proposes improved methods for this important problem. One of the most popular state-of-the-art methods is Parallel Analysis (PA), which compares the observed factor strengths to simulated ones under a noise-only model. This paper proposes improvements to PA. We first de-randomize it, proposing Deterministic Parallel Analysis (DPA), which is faster and more reproducible than PA. Both PA and DPA are prone to a shadowing phenomenon in which a strong factor makes it hard to detect smaller but more interesting factors. We propose deflation to counter shadowing. We also propose to raise the decision threshold to improve estimation accuracy. We prove several consistency results for our methods, and test them in simulations. We also illustrate our methods on data from the Human Genome Diversity Project, where they significantly improve the accuracy., Comment: Made title consistent with published version

Researchers often have datasets measuring features $x_{ij}$ of samples, such as test scores of students. In factor analysis and PCA, these features are thought to be influenced by unobserved factors, such as skills. Can we determine how many components affect the data? This is an important problem, because it has a large impact on all downstream data analysis. Consequently, many approaches have been developed to address it. Parallel Analysis is a popular permutation method. It works by randomly scrambling each feature of the data. It selects components if their singular values are larger than those of the permuted data. Despite widespread use in leading textbooks and scientific publications, as well as empirical evidence for its accuracy, it currently has no theoretical justification. In this paper, we show that the parallel analysis permutation method consistently selects the large components in certain high-dimensional factor models. However, it does not select the smaller components. The intuition is that permutations keep the noise invariant, while "destroying" the low-rank signal. This provides justification for permutation methods in PCA and factor models under some conditions. Our work uncovers drawbacks of permutation methods, and paves the way to improvements., Comment: To appear in the Annals of Statistics

We consider the linearly transformed spiked model, where observations $Y_i$ are noisy linear transforms of unobserved signals of interest $X_i$: \begin{align*} Y_i = A_i X_i + \varepsilon_i, \end{align*} for $i=1,\ldots,n$. The transform matrices $A_i$ are also observed. We model $X_i$ as random vectors lying on an unknown low-dimensional space. How should we predict the unobserved signals (regression coefficients) $X_i$? The naive approach of performing regression for each observation separately is inaccurate due to the large noise. Instead, we develop optimal linear empirical Bayes methods for predicting $X_i$ by "borrowing strength" across the different samples. Our methods are applicable to large datasets and rely on weak moment assumptions. The analysis is based on random matrix theory. We discuss applications to signal processing, deconvolution, cryo-electron microscopy, and missing data in the high-noise regime. For missing data, we show in simulations that our methods are faster, more robust to noise and to unequal sampling than well-known matrix completion methods., Comment: This paper replaces the preprint "PCA from noisy, linearly reduced data: the diagonal case" by Edgar Dobriban, William Leeb, and Amit Singer (arXiv:1611.10333)

Suppose we observe data of the form $Y_i = D_i (S_i + \varepsilon_i) \in \mathbb{R}^p$ or $Y_i = D_i S_i + \varepsilon_i \in \mathbb{R}^p$, $i=1,\ldots,n$, where $D_i \in \mathbb{R}^{p\times p}$ are known diagonal matrices, $\varepsilon_i$ are noise, and we wish to perform principal component analysis (PCA) on the unobserved signals $S_i \in \mathbb{R}^p$. The first model arises in missing data problems, where the $D_i$ are binary. The second model captures noisy deconvolution problems, where the $D_i$ are the Fourier transforms of the convolution kernels. It is often reasonable to assume the $S_i$ lie on an unknown low-dimensional linear space; however, because many coordinates can be suppressed by the $D_i$, this low-dimensional structure can be obscured. We introduce diagonally reduced spiked covariance models to capture this setting. We characterize the behavior of the singular vectors and singular values of the data matrix under high-dimensional asymptotics where $n,p\to\infty$ such that $p/n\to\gamma>0$. Our results have the most general assumptions to date even without diagonal reduction. Using them, we develop optimal eigenvalue shrinkage methods for covariance matrix estimation and optimal singular value shrinkage methods for data denoising. Finally, we characterize the error rates of the empirical Best Linear Predictor (EBLP) denoisers. We show that, perhaps surprisingly, their optimal tuning depends on whether we denoise in-sample or out-of-sample, but the optimally tuned mean squared error is the same in the two cases., Comment: This technical report has been largely superseded by our later paper arXiv:1709.03393. Please cite that one instead of this one. This paper has a slightly different approach, so we want to keep it publicly available

Many applications, such as photon-limited imaging and genomics, involve large datasets with noisy entries from exponential family distributions. It is of interest to estimate the covariance structure and principal components of the noiseless distribution. Principal Component Analysis (PCA), the standard method for this setting, can be inefficient when the noise is non-Gaussian. We develop $e$PCA (exponential family PCA), a new methodology for PCA on exponential family distributions. $e$PCA can be used for dimensionality reduction and denoising of large data matrices. $e$PCA involves the eigendecomposition of a new covariance matrix estimator, constructed in a simple and deterministic way using moment calculations, shrinkage, and random matrix theory. We provide several theoretical justifications for our estimator, including the finite-sample convergence rate, and the Marchenko-Pastur law in high dimensions. $e$PCA compares favorably to PCA and various PCA alternatives for exponential families, in simulations as well as in XFEL and SNP data analysis. An open-source implementation is available., Comment: 24 pages, 9 figures. An open-source implementation can be found at http://github.com/lydiatliu/epca/

Researchers in data-rich disciplines---think of computational genomics and observational cosmology---often wish to mine large bodies of $p$-values looking for significant effects, while controlling the false discovery rate or family-wise error rate. Increasingly, researchers also wish to prioritize certain hypotheses, for example those thought to have larger effect sizes, by upweighting, and to impose constraints on the underlying mining, such as monotonicity along a certain sequence. We introduce Princessp, a principled method for performing weighted multiple testing by constrained convex optimization. Our method elegantly allows one to prioritize certain hypotheses through upweighting and to discount others through downweighting, while constraining the underlying weights involved in the mining process. When the $p$-values derive from monotone likelihood ratio families like the Gaussian means model, the new method allows exact solution of an important optimal weighting problem previously thought to be nonconvex and computationally infeasible. Our method scales to massive dataset sizes. We illustrate the applications of Princessp on a series of standard genomics datasets and offer comparisons with several previous `standard' methods. Princessp offers both ease of operation and the ability to scale to extremely large problem sizes. The method is available as open-source software from http://github.com/dobriban/pvalue_weighting_matlab ., Comment: This is an entirely rewritten version of the paper. The title of the paper, the name of the method, and the introduction have been changed, with the goal of making the paper more accessible and appealing to practitioners. New sections on monotone likelihood ratio families and two-sided tests have been added, which expand the scope of the method

Principal component analysis (PCA) is a widely used method for dimension reduction. In high dimensional data, the "signal" eigenvalues corresponding to weak principal components (PCs) do not necessarily separate from the bulk of the "noise" eigenvalues. Therefore, popular tests based on the largest eigenvalue have little power to detect weak PCs. In the special case of the spiked model, certain tests asymptotically equivalent to linear spectral statistics (LSS)---averaging effects over all eigenvalues---were recently shown to achieve some power. We consider a nonparametric, non-Gaussian generalization of the spiked model to the setting of Marchenko and Pastur (1967). This allows a general bulk of the noise eigenvalues, accomodating correlated variables even under the null hypothesis of no significant PCs. We develop new tests based on LSS to detect weak PCs in this model. We show using the CLT for LSS that the optimal LSS satisfy a Fredholm integral equation of the first kind. We develop algorithms to solve it, building on our recent method for computing the limit empirical spectrum. In contrast to the standard spiked model, we find that under "widely spread" null eigenvalue distributions, the new tests have a lot of power., Comment: 46 pages, 9 figures

We provide a unified analysis of the predictive risk of ridge regression and regularized discriminant analysis in a dense random effects model. We work in a high-dimensional asymptotic regime where $p, n \to \infty$ and $p/n \to \gamma \in (0, \, \infty)$, and allow for arbitrary covariance among the features. For both methods, we provide an explicit and efficiently computable expression for the limiting predictive risk, which depends only on the spectrum of the feature-covariance matrix, the signal strength, and the aspect ratio $\gamma$. Especially in the case of regularized discriminant analysis, we find that predictive accuracy has a nuanced dependence on the eigenvalue distribution of the covariance matrix, suggesting that analyses based on the operator norm of the covariance matrix may not be sharp. Our results also uncover several qualitative insights about both methods: for example, with ridge regression, there is an exact inverse relation between the limiting predictive risk and the limiting estimation risk given a fixed signal strength. Our analysis builds on recent advances in random matrix theory., Comment: Added a section on prediction versus estimation for ridge regression. Rewrote introduction. Other results unchanged

Consider an $n \times p$ data matrix $X$ whose rows are independently sampled from a population with covariance $\Sigma$. When $n,p$ are both large, the eigenvalues of the sample covariance matrix are substantially different from those of the true covariance. Asymptotically, as $n,p \to \infty$ with $p/n \to \gamma$, there is a deterministic mapping from the population spectral distribution (PSD) to the empirical spectral distribution (ESD) of the eigenvalues. The mapping is characterized by a fixed-point equation for the Stieltjes transform. We propose a new method to compute numerically the output ESD from an arbitrary input PSD. Our method, called Spectrode, finds the support and the density of the ESD to high precision; we prove this for finite discrete distributions. In computational experiments it outperforms existing methods by several orders of magnitude in speed and accuracy. We apply Spectrode to compute expectations and contour integrals of the ESD. These quantities are often central in applications of random matrix theory (RMT). We illustrate that Spectrode is directly useful in statistical problems, such as estimation and hypothesis testing for covariance matrices. Our proposal may make it more convenient to use asymptotic RMT in aspects of high-dimensional data analysis.

We develop a new method for frequentist multiple testing with Bayesian prior information. Our procedure finds a new set of optimal p-value weights called the Bayes weights. Prior information is relevant to many multiple testing problems. Existing methods assume fixed, known effect sizes available from previous studies. However, the case of uncertain information is usually the norm. For a Gaussian prior on effect sizes, we show that finding the optimal weights is a non-convex problem. Despite the non-convexity, we give an efficient algorithm that solves this problem nearly exactly. We show that our method can discover new loci in genome-wide association studies. On several data sets it compares favorably to other methods. Open source code is available.

Statistical and machine learning theory has developed several conditions ensuring that popular estimators such as the Lasso or the Dantzig selector perform well in high-dimensional sparse regression, including the restricted eigenvalue, compatibility, and $\ell_q$ sensitivity properties. However, some of the central aspects of these conditions are not well understood. For instance, it is unknown if these conditions can be checked efficiently on any given data set. This is problematic, because they are at the core of the theory of sparse regression. Here we provide a rigorous proof that these conditions are NP-hard to check. This shows that the conditions are computationally infeasible to verify, and raises some questions about their practical applications. However, by taking an average-case perspective instead of the worst-case view of NP-hardness, we show that a particular condition, $\ell_q$ sensitivity, has certain desirable properties. This condition is weaker and more general than the others. We show that it holds with high probability in models where the parent population is well behaved, and that it is robust to certain data processing steps. These results are desirable, as they provide guidance about when the condition, and more generally the theory of sparse regression, may be relevant in the analysis of high-dimensional correlated observational data., Comment: Manuscript shortened and more motivation added. To appear in Communications in Mathematics and Statistics

This paper is concerned with an important matrix condition in compressed sensing known as the restricted isometry property (RIP). We demonstrate that testing whether a matrix satisfies RIP is NP-hard. As a consequence of our result, it is impossible to efficiently test for RIP provided P \neq NP.