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?
is addressed.
Reference
[Anderson]Anderson, Theodore Wilbur. "Asymptotic theory for principal component analysis." The Annals of Mathematical Statistics 34.1 (1963): 122-148.
[link1]https://math.stackexchange.com/questions/1467581/what-does-it-mean-for-a-probability-distribution-to-have-a-finite-fourth-moment
[Vershynin]http://www-personal.umich.edu/~Vershyninv/papers/non-asymptotic-rmt-plain.pdf
[Hastie&Tibshirani]Trevor Hastie and Robert Tibshirani "Statistical Learning with Sparsity: The Lasso and Generalizations."
[link2]https://stats.stackexchange.com/questions/263879/is-an-admissible-minimax-rule-unique
[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.
[link3]angle between subspaces
[Jiang]Jiang, Tiefeng. "Approximation of Haar distributed matrices and limiting distributions of eigenvalues of Jacobi ensembles." Probability theory and related fields 144.1-2 (2009): 221-246.
[Diaconis&Shahshahani]Diaconis, Persi, and Mehrdad Shahshahani. "On the eigenvalues of random matrices." Journal of Applied Probability (1994): 49-62.http://statweb.stanford.edu/~cgates/PERSI/papers/random_matrices.pdf
[Pastur et.al]Pastur, Leonid Andreevich, Mariya Shcherbina, and Mariya Shcherbina. Eigenvalue distribution of large random matrices. Vol. 171. Providence, RI: American Mathematical Society, 2011.
[Deift]Deift, Percy. Orthogonal polynomials and random matrices: a
Riemann-Hilbert approach. Vol. 3. American Mathematical Soc., 2000.
[Bandeira&van Handel]Bandeira, Afonso S., and Ramon van Handel. "Sharp nonasymptotic bounds on the norm of random matrices with independent entries." The Annals of Probability 44.4 (2016): 2479-2506.
The conjecture is true. Moreover we show that
if $W(N,N,1)$ is replaced by any unitary $N \times N$ matrix
then $\max(2k-N,0)$ is still a lower bound on the multiplicity,
and is usually sharp (though for $k<N$ there are easy exceptions
such as the $N \times N$ identity matrix).
First observe that if $T$ is any "sub-unitary" matrix
(i.e. $|Tv| \leq |v|$ for all $v$)
then the vectors $v$ such that $|Tv| = |v|$ form a vector space
whose dimension is the multiplicity of $1$ as a singular value of $T$.
Indeed that multiplicity is the dimension of the $1$-eigenspace of $T^* T$;
but if $T^* T v = v$ then $|Tv|^2 = (Tv,Tv) = (v, T^* T v) = |v|^2$,
and conversely if $|v| = |Tv|$ then
$|v|^2 = (v, T^* T v) \leq \left|v\right| \left|T^*T v\right| \leq |v|^2$
(since $\|T^* T\| \leq 1$), with equality only if $T^* T v = v$.
Now suppose $T$ is obtained as the principal $k \times k$ minor of
an $N \times N$ unitary matrix $U$. Then for any $v \in {\bf C}^k$
we find $Tv$ by padding $v$ by $N-k$ zeros to a vector $\bar v \in {\bf C}^n$,
applying $U$, and discarding the last $N-k$ entries of the result $U \bar v$.
Therefore $|Tv| \leq |v|$ with equality if and only if
$U\bar v$ is supported on the first $k$ coordinates.
This identifies the space $\{ v : |Tv| = |v| \}$ with
the kernel of the $k \times (N-k)$ submatrix of $U$
formed by intersecting the first $k$ rows with the last $N-k$ columns.
Since this submatrix has rank at most $\min(N-k,k)$,
its kernel has dimension at least $\max(2k-N,0)$.
In the special case that $U$ is the matrix of the discrete Fourier transform,
any square submatrix supported on the first $m$ rows or the first $m$ columns
is a Vandermonde matrix or its transpose, and is thus nonsingular.
Therefore our $k \times (N-k)$ submatrix has a nonsingular square submatrix
of order $\min(N-k,k)$. Hence it attains the upper bound $\min(N-k,k)$
on the rank of a $k \times (N-k)$ matrix, and the singular value $1$
has the lowest possible multiplicity $\max(2k-N,0)$, QED.
To show that this holds for "typical" unitary $U$, we need only observe that
$m \times n$ matrices of rank $\min(m,n)$ form a nonempty open subset of
the space of $m \times n$ matrices. Hence unitary matrices $U$ whose
relevant $k \times (N-k)$ submatrix satisfies this condition form
the preimage of an open set under a continuous map, and therefore
constitute an open subset of the group ${\rm U}_N$ of unitary matrices.
This open subset is nonempty because it contains the discrete
Fourier transform, and is thus dense because ${\rm U}_N$ is connected.
Best Answer
Your numerical findings have a simple explanation: Let me denote the upper left block of the $N\times N$ unitary matrix by $R$ and the upper right block by $T$; the matrix $R$ has dimension $fN\times fN$, the upper right block has dimension $fN\times(1-f)N$; unitarity requires that $$RR^\dagger+TT^\dagger=I,$$ with $\dagger$ indicating the conjugate transpose and $I$ the identity matrix. For $f\geq 1/2$ the matrix $TT^\dagger$ has rank $\leq (1-f)N$, so it has at least $fN-(1-f)N=(2f-1)N$ eigenvalues equal to zero.
Hence for $f\geq 1/2$ the matrix $RR^\dagger$ has at least $(2f-1)N$ eigenvalues equal to 1, or equivalently, $R$ has $\geq (2f-1)N$ singular values equal to unity. For generic matrices the $\geq$ will be an equality.