# What kind of random matrices have rapidly decaying singular values?

I've been told that in machine learning it's common to compute the singular value decomposition of matrices in order to throw out all information in the matrix except that corresponding to, say, the $k$ largest singular values for some reasonable $k$; moreover, in practice it happens that the singular values decay rapidly and so all but the $k$ largest singular values contain very little information anyway.

What kind of distribution on random matrices has this kind of distribution of singular values?

For example, classical results (Theorems 5.31 and 5.32) imply that for an $N \times n$ matrix ($N \ge n$) whose entries are, say, i.i.d. Gaussian with zero mean and unit variance, the largest singular value is around $\sqrt{N} + \sqrt{n}$ with high probability and the smallest singular value is around $\sqrt{N} - \sqrt{n}$ with high probability; in particular, for $N \gg n$ all of the singular values are pretty similar. This is a very different distribution from the above, so it seems reasonable to infer that the matrices people encounter in machine learning don't behave like they have i.i.d. entries.

Probably the independence assumption is the least reasonable one - for example, in an incidence matrix describing which users of Netflix have watched which movies, one expects lots of correlations between similar movies as well as between similar users; that is, one expects the incidence matrix to be well approximated by a low-rank matrix (which is precisely what we're doing above). I guess one also expects power-law (log-Gaussian?) distributions on, say, both the statistics of how many movies users have watched and how many users have watched movies. Have random models with this kind of behavior been studied, and do their singular values behave as expected?

• A standard model is of Wishart matrices where $X\Sigma X^T$ where $X$ are iid Gaussians and $\Sigma$ is diagonal positive; look up spiked models'' for the case where most entries of $\Sigma$ equal to $1$ except for a few larger one. There are many results for such matrices, both on convergence of empirical measures as well as for the behavior of a few large ones. – ofer zeitouni Feb 25 '14 at 2:31
• Power law graphs are typical in machine learning (google finds links to spectra of such graphs); also, the most heavily used gaussian kernel $\exp(-\gamma\|x-y\|^2$ usually leads to kernel matrices with rapidly decaying spectrum; some other kernel functions that are used also have this property. – Suvrit Feb 25 '14 at 4:26
• @Suvrit: interesting. It looks like one way to define these power law random graphs is to impose a power law distribution on degrees. I wonder if there's a general result to the effect that imposing a particular distribution on degrees gives a similar distribution on eigenvalues or singular values... – Qiaochu Yuan Feb 25 '14 at 23:13
• @Qiaochu: I did not get time to chase the details yet, but the following talk seems promising: aladdin.cs.cmu.edu/workshops/wsa/papers/spectrum.pdf – Suvrit Feb 26 '14 at 15:33
• Parameter sensitivity matrices (Hessian of the goodness-of-fit metric) in many models of physical and biological systems exhibit such a an eigenvalue distribution. This phenomenon goes by the name "sloppy models" and has a substantial literature explaining the phenomenon and developing methods to work with sloppy models (see, e.g., science.sciencemag.org/content/342/6158/604) – Yoav Kallus Apr 7 '17 at 19:15

I hope I understood the OP correctly, in case not please let me know. And I will discuss the case of eigenvalue instead of singular value without much loss of generality.

In case you are only interested in the distribution of eigenvalues of an i.i.d. sample covariance matrix(which is also a random matrix but underlying distribution is always of dimension one), the correct reference should be [Anderson] instead of [Johnstone]. [Johnstone] dealt with a much more general case where the underlying distribution is of multi-dimension and especially with double-Wishart situation and its associated concentration bound. On the other hand [Anderson] did not consider $n(p),m(p)$ but only $m,n$ as separated parameters that tend to $\infty$, which is not suitable in "machine learning era".

$\blacksquare$1.Theorem 5.31 in [Vershynin] and Theorem 5.32 as corollary.

Unlike what the OP claimed, Bai's result also requires finite fourth moment [link1] which is not quite easy to verify in some applications, say the censored data. Assuming for the moment that the fourth moment of the underlying distribution is finite and consider the simple Gaussian regression for example, we usually study $N\approx n$ where $N$ is the sample size and $n$ is the number of variables in the model. You are right on the point if $N\gg n$ then all eigenvalues concentrate on a certain neighborhood of $\sqrt{N}$. However, a challenge in the field is that when $n$ is large and $N\approx n$. Intuitively, when the dimension of the regression model is high and the data is sparse. In this situation, some "machine learning" methods were proposed to select variables [Hastie&Tibshirani]. Among these methods, one is called sparse PCA which penalized sparse input and "forces the eigenvalues to deviate from $N$ with positive probability" if you want to view in Bai's law (Theorem 5.31).

So it is not that

"so it seems reasonable to infer that the (random) matrices people encounter in machine learning don't behave like they have i.i.d. entries."

instead researchers will treat i.i.d. in companion with consideration of the dimension of the space that the working model lied in. Dependent data will look almost like i.i.d data if the dimension is high. [Pastur et.al] The consideration of dimension is not unusual even in classical statistic, some admissible models will become inadmissible if the dimension increase, for example Stein's effect [link2].

$\blacksquare$2. What kind of distribution on random matrices has this kind of distribution of singular values(eigenvalues)?

2.1 Behavior of eigenvalues of dependent Gaussian random matrices

The independent assumption is not necessary as you already noted. However, even if we dropped the independence we can still obtain a concentration bound for the largest eigenvalue. Consider $N(\boldsymbol{0},\Sigma)$ instead of $N(\boldsymbol{0},I)$. Again we consider regression setting, the following result tells you that without independence between variables we can still impose a bound on the covariance matrix which is required to proceed PCA or other variable selection techniques using eigenvalue decomposition. Notice that $Wishart(I,\bullet)$ is not a serious constraint as commented by [Johnstone] p.3.

Thm 1 in [Johnstone](modified)

Assume $p$ is even and that $p, m = m(p)$ and $n= n(p)\rightarrow\infty$ together in such a way that $\lim_{p\rightarrow\infty}\frac{\min(n,p)}{(m+n)}>0,\lim_{p\rightarrow\infty}\frac{p}{m}<1$. For each $s_{0}\in\mathbb{R}$, there exists $C > 0$ such that for $s\geq s_0$, $$|Pr\{W_p\leq\mu_p+\sigma_p\cdot s\}-F_1(s)|<Cp^{-2/3}e^{-s/2}$$ where $W_p=\log\left(\frac{\theta_{1,p}}{1-\theta_{1,p}}\right)$ and $\theta_{1,p}$ is the largest eigenvalue of $(A + B)^{-1}B$ where $A\sim Wishart_p(I,m),B\sim Wishart_p(I,n)$ are two random matrices. Here $C,\mu_p,\sigma_p$ depends ONLY on $(\gamma, \phi)$ defined as canonical angle parameters $$\sin^2(\frac{\gamma}{2})=\frac{\min(p,n)-.5}{m+n-1}$$ $$\sin^2(\frac{\phi}{2})=\frac{\max(p,n)-.5}{m+n-1}$$ and also on $s_0$ if $s_0 < 0$. $F_1$ is the Tracy-Widom distribution.

To see more clearly how this will work in a regression case, we can actually see from definition that $\theta_{1,p}$ defined above is the solution to the generalized eigenvalue problem for random matrices $A,B$ such that $\det[x(A + B) − A] = 0$. A special case is $\det[P_{X}x − I] = 0$ where $P_X$ is chosen to be a projection matrix(OR the estimated covariance matrix) in regression setting. But as expected as we saw in the result above, as dimension of the model $p$ increase, such a convergence will slow down at the rate of $O(p^{-2/3})$, which allows more freedom for "extremely large/small eigen values".

I plotted an illustration as below (NOT Exact simulation!)

So for most Gaussian data, the largest eigenvalues will exponentially converge to a certain distribution, which means that the decaying rate of the largest eigenvalues is determined by the family of Tracy-Widom distribution indexed by (sample size) $m(p),n(p)$. Moreover we can even study the speed of convergence in a more precise way using the "central limit theorem"(circular law) for random matrices such as [Jiang].

[Bandeira&van Handel] now gives a sharp bound on spectral norm.

2.2 Behavior of eigenvalues of product uniform random matrices

Results of similar spirit concerning non-Gaussian data could be found in [Absil et.al]. In that paper the authors concern about the uniform distribution over a Grassmanian and their Theorem 1 (not copied here) also showed that the distribution of canonical angles between the spaces spanned by columns of two uniformly distributed random matrices $A,B$ are also dominated by Tracy-Widom distribution (multiplied by a gamma-polynomial term in $p,n(p)$). The cosines of canonical angles are exactly the singular values of $A'B$.[link3, Absil et.al]. An earlier related result is obtained in [Diaconis&Shahshahani] which considers only unitary random matrices.

Therefore the relative magnitude $\frac{\lambda_1}{\lambda_2}$ instead of $|\lambda_1-\lambda_2|$ are used as criterion when we choose which $k$ principle components are to remain in the model. This also explains a well-known drawback of PCA, which is that PCA is scale-dependent.

2.3 Conclusion

The asymptotic distribution of eigenvalues of a (OR asymptotically w.r.t. $m,n$) Gaussian random matrices is dominated by the Tracy-Widom distribution and its convergence behavior is controlled by an exponential rate associated with deviation $s$ and a polynomial rate associated with dimension $p$. Therefore whether the eigenvalues of a given random matrix concentrates together depends on the $m(p),n(p),p$ and eventually on the dimension $p$.

$\bullet$ The eigenvalues of a given Gaussian random matrix will not decay rapidly if we make the growth of the sample size satisfy $m(p),n(p)\rightarrow\infty$ in such a way that $$\lim_{p\rightarrow\infty}\frac{\min(n,p)}{(m+n)}>0,\lim_{p\rightarrow\infty}\frac{p}{m}<1$$. If it decays rapidly, we have some existing methods to penalized it as mentioned in part 1.

$\bullet$ The eigenvalues of product uniform distributions (on Grassmanians of dim $p$) will not decay rapidly under certain constraints.

Therefore it is safe to say that in most applications, as long as $m(p),n(p)$(usually sample sizes) grows with $p$(usually model dimension) in a reasonable manner, the eigenvalues will neither concentrate severely NOR change their relative magnitude too dramatically.

On the other hand, if a random matrix is actually generated from a model with varying (increasing) sample size $N$, then it will most likely have a set of rapidly decaying eigenvalues under mild constraints as discussed above.

In case where the growth conditions are not satisfied, the word "sparse" come into play.

$\blacksquare$3. Return to OP's last question.

Therefore the focus of (statistics) research is not on which family of random models have rapidly decaying eigenvalues such that we can avoid using variable selection on these "morbid models"; the focus of (some) statistical methods are on how to wrap up these possible morbidity and make variable selection based on such eigenvalues still work.

Later I have a nice discussion with some expert in the field, who pointed out that the eigenvalue problem can be framed as the distribution of roots of a polynomial with random coefficients(coefficients as random variables) and simultaneously a kernel if the covariance matrix is p.d.

The distribution of eigenvalues of a random matrix is described in terms of the distribution of roots of a determinant equation, which turns out to be a random polynomial[Deift].

And by now I think your question that

Have random models with this kind of behavior been studied, and do their singular values behave as expected?

Reference [Anderson]Anderson, Theodore Wilbur. "Asymptotic theory for principal component analysis." The Annals of Mathematical Statistics 34.1 (1963): 122-148.

[Hastie&Tibshirani]Trevor Hastie and Robert Tibshirani "Statistical Learning with Sparsity: The Lasso and Generalizations."

[Johnstone]Johnstone, Iain M. "Multivariate analysis and Jacobi ensembles: Largest eigenvalue, Tracy–Widom limits and rates of convergence." Annals of statistics 36.6 (2008): 2638.(Slides: http://www.stat.harvard.edu/NESS10/Johnstone-slides.pdf)

[Tracy&Widom]Tracy, Craig A., and Harold Widom. "On orthogonal and symplectic matrix ensembles." Communications in Mathematical Physics 177.3 (1996): 727-754.

[Absil et.al]Absil, P-A., Alan Edelman, and Plamen Koev. "On the largest principal angle between random subspaces." Linear Algebra and its applications 414.1 (2006): 288-294.