Let's say I have a $p$-dimensional multivariate Gaussian distribution. And I take $n$ observations (each of them a $p$-vector) from this distribution and calculate the sample covariance matrix $S$. In this paper, the authors state that the sample covariance matrix calculated with $p > n$ is singular.
Some facts about matrix ranks, offered without proof (but proofs of all or almost all of them should be either given in standard linear algebra texts, or in some cases set as exercises after giving enough information to be able to do so):
If $A$ and $B$ are two conformable matrices, then:
(i) column rank of $A$ = row rank of $A$
(ii) $\text{rank}(A) = \text{rank}(A^T) = \text{rank}(A^TA) = \text{rank}(AA^T)$
(iii) $\text{rank}(AB)\leq \min(\text{rank}(A),\text{rank}(B))$
(iv) $\text{rank}(A+B) \leq \text{rank}(A) + \text{rank}(B)$
(v) if $B$ is square matrix of full rank, then $\text{rank}(AB) = \text{rank}(A)$
Consider the $n\times p$ matrix of sample data, $y$. From the above, the rank of $y$ is at most $\min(n,p)$.
Further, from the above clearly the rank of $S$ won't be larger than the rank of $y$ (by considering the computation of $S$ in matrix form, with perhaps some simplification).
If $n<p$ then $\text{rank}(y)<p$ in which case $\text{rank}(S)<p$.
The short answer to your question is that rank$(S) \le n - 1$. So if $p > n$, then $S$ is singular.
For a more detailed answer, recall that the (unbiased) sample covariance matrix can be written as
$$ S = \frac{1}{n-1}\sum_{i=1}^n (x_i - \bar{x})(x_i - \bar{x})^T. $$
Effectively, we are summing $n$ matrices, each having a rank of 1. Assuming the observations are linearly independent, in some sense each observation $x_i$ contributes 1 to rank$(S)$, and a 1 is subtracted from the rank (if $p > n$) because we center each observation by $\bar{x}$. However, if multicollinearity is present in the observations, then rank$(S)$ may be reduced, which explains why the rank might be less than $n - 1$.
A large amount of work has gone into studying this problem. For instance, a colleague of mine and I wrote a paper on this same topic, where we were interested in determining how to proceed if $S$ is singular when applied to linear discriminant analysis in the $p \gg n$ setting.
When you look at the situation the right way, the conclusion is intuitively obvious and immediate.
This post offers two demonstrations. The first, immediately below, is in words. It is equivalent to a simple drawing, appearing at the very end. In between is an explanation of what the words and the drawing mean.
The covariance matrix for $n$ $p$-variate observations is a $p\times p$ matrix computed by left-multiplying a matrix $\mathbb{X}_{np}$ (the recentered data) by its transpose $\mathbb{X}_{pn}^\prime$. This product of matrices sends vectors through a pipeline of vector spaces in which the dimensions are $p$ and $n$. Consequently the covariance matrix, qua linear transformation, will send $\mathbb{R}^n$ into a subspace whose dimension is at most $\min(p,n)$. It is immediate that the rank of the covariance matrix is no greater than $\min(p,n)$. Consequently, if $p\gt n$ then the rank is at most $n$, which--being strictly less than $p$--means the covariance matrix is singular.
All this terminology is fully explained in the remainder of this post.
(As Amoeba kindly pointed out in a now-deleted comment, and shows in an answer to a related question, the image of $\mathbb X$ actually lies in a codimension-one subspace of $\mathbb{R}^n$ (consisting of vectors whose components sum to zero) because its columns have all been recentered at zero. Therefore the rank of the sample covariance matrix $\frac{1}{n-1}\mathbb{X}^\prime \mathbb{X}$ cannot exceed $n-1$.)
Linear algebra is all about tracking dimensions of vector spaces. You only need to appreciate a few fundamental concepts to have a deep intuition for assertions about rank and singularity:
Matrix multiplication represents linear transformations of vectors. An $m\times n$ matrix $\mathbb{M}$ represents a linear transformation from an $n$-dimensional space $V^n$ to an $m$-dimensional space $V^m$. Specifically, it sends any $x\in V^n$ to $\mathbb{M}x = y \in V^m$. That this is a linear transformation follows immediately from the definition of linear transformation and basic arithmetical properties of matrix multiplication.
Linear transformations can never increase dimensions. This means that the image of the entire vector space $V^n$ under the transformation $\mathbb M$ (which is a sub-vector space of $V^m$) can have a dimension no greater than $n$. This is an (easy) theorem that follows from the definition of dimension.
The dimension of any sub-vector space cannot exceed that of the space in which it lies. This is a theorem, but again it is obvious and easy to prove.
The rank of a linear transformation is the dimension of its image. The rank of a matrix is the rank of the linear transformation it represents. These are definitions.
A singular matrix $\mathbb{M}_{mn}$ has rank strictly less than $n$ (the dimension of its domain). In other words, its image has a smaller dimension. This is a definition.
To develop intuition, it helps to see the dimensions. I will therefore write the dimensions of all vectors and matrices immediately after them, as in $\mathbb{M}_{mn}$ and $x_n$. Thus the generic formula
$$y_m = \mathbb{M}_{mn} x_n$$
is intended to mean that the $m\times n$ matrix $\mathbb M$, when applied to the $n$-vector $x$, produces an $m$-vector $y$.
Products of matrices can be thought of as a "pipeline" of linear transformations. Generically, suppose $y_a$ is an $a$-dimensional vector resulting from the successive applications of the linear transformations $\mathbb{M}_{mn}, \mathbb{L}_{lm}, \ldots, \mathbb{B}_{bc},$ and $\mathbb{A}_{ab}$ to the $n$-vector $x_n$ coming from the space $V^n$. This takes the vector $x_n$ successively through a set of vector spaces of dimensions $m, l, \ldots, c, b,$ and finally $a$.
Look for the bottleneck: because dimensions cannot increase (point 2) and subspaces cannot have dimensions larger than the spaces in which they lie (point 3), it follows that the dimension of the image of $V^n$ cannot exceed the smallest dimension $\min(a,b,c,\ldots,l,m,n)$ encountered in the pipeline.
This diagram of the pipeline, then, fully proves the result when it is applied to the product $\mathbb{X}^\prime \mathbb{X}$:
