  Introduced by Kiefer and Wolfowitz \cite{KW56}, the nonparametric maximum
likelihood estimator (NPMLE) is a widely used methodology for learning mixture
odels and empirical Bayes estimation. Sidestepping the non-convexity in mixture
likelihood, the NPMLE estimates the mixing distribution by maximizing the total
likelihood over the space of probability measures, which can be viewed as an
extreme form of overparameterization.
  In this paper we discover a surprising property of the NPMLE solution.
Consider, for example, a Gaussian mixture model on the real line with a
subgaussian mixing distribution. Leveraging complex-analytic techniques, we
show that with high probability the NPMLE based on a sample of size $n$ has
$O(\log n)$ atoms (mass points), significantly improving the deterministic
upper bound of $n$ due to Lindsay \cite{lindsay1983geometry1}. Notably, any
such Gaussian mixture is statistically indistinguishable from a finite one with
$O(\log n)$ components (and this is tight for certain mixtures). Thus, absent
any explicit form of model selection, NPMLE automatically chooses the right
model complexity, a property we term \emph{self-regularization}. Extensions to
other exponential families are given. As a statistical application, we show
that this structural property can be harnessed to bootstrap existing Hellinger
risk bound of the (parametric) MLE for finite Gaussian mixtures to the NPMLE
for general Gaussian mixtures, recovering a result of Zhang
\cite{zhang2009generalized}.
