  The Hilbert scheme of points $\mathrm{Hilb}^n(S)$ of a smooth surface $S$ is
a well-studied parameter space, lying at the interface of algebraic geometry,
commutative algebra, representation theory, combinatorics, and mathematical
physics. The foundational result is a classical theorem of Fogarty, stating
that $\mathrm{Hilb}^n(S)$ is a smooth variety of dimension $2n$. In recent
years there has been growing interest in a natural generalization of
$\mathrm{Hilb}^n(S)$, the nested Hilbert scheme $\mathrm{Hilb}^{(n_1,n_2)}(S)$,
which parametrizes nested pairs of zero-dimensional subschemes $Z_1 \supseteq
Z_2$ of $S$ with $\mathrm{deg} (Z_i)=n_i$. In contrast to Fogarty's theorem,
$\mathrm{Hilb}^{(n_1,n_2)}(S)$ is almost always singular, and very little is
known about its singularities. In this paper we aim to advance the knowledge of
the geometry of these nested Hilbert schemes. Work by Fogarty in the 70's shows
that $\mathrm{Hilb}^{(n,1)}(S)$ is a normal Cohen-Macaulay variety, and Song
more recently proved that it has rational singularities. In our main result, we
prove that the nested Hilbert scheme $\mathrm{Hilb}^{(n,2)}(S)$ has rational
singularities. We employ an array of tools from commutative algebra to prove
this theorem. Using Gr\"obner bases, we establish a connection between
$\mathrm{Hilb}^{(n,2)}(S)$ and a certain variety of matrices with an action of
the general linear group. This variety of matrices plays a central role in our
work, and we analyze it by various algebraic techniques, including square-free
Gr\"obner degenerations, the Stanley-Reisner correspondence, and the
Kempf-Lascoux-Weyman technique of calculating syzygies. Along the way, we also
obtain results on classes of irreducible and reducible nested Hilbert schemes,
dimension of singular loci, and $F$-singularities in positive characteristic.
