Any covariance matrix is positive semi-definite: A proof

Let $\mathbf{X}=(X_1,X_2,\cdots ,X_n{)}^{\mathrm{\top }}\in {\mathbb{R}}^n$ be a vector of random variables. The covariance matrix $\mathrm{\Sigma }$ of $\mathbf{X}$ is a square ($n\times n$) matrix whose elements are covariances between the components of $\mathbf{X}$. That is,

$${\mathrm{\Sigma }}_{ij}=\mathrm{Cov}(X_i,X_j)$$

where $\mathrm{Cov}(X_i,X_j)$ is the covariance between $X_i$ and $X_j$, $i,j=1,2,\cdots ,n$: 

$$\mathrm{Cov}(X_i,X_j)=\mathbb{E}[(X_i-\mathbb{E}[X_i])(X_j-\mathbb{E}[X_j])].$$

Here, $\mathbb{E}[\cdot ]$ indicates the expectation value of a random variable. Any covariance matrix has the following properties:

1. Symmetric. That is, $$\mathrm{\Sigma }={\mathrm{\Sigma }}^{\mathrm{\top }}.$$
2. Positive semi-definite. That is,$$\mathrm{\forall }\mathbf{v}\in {\mathbb{R}}^n,{\mathbf{v}}^{\mathrm{\top }}\mathrm{\Sigma }\mathbf{v}\ge 0.$$

See also: Positive definite matrix (Wolfram MathWorld)

The symmetry is obvious from the definition of the covariance matrix. 

Now, let us prove that the covariance matrix is positive semi-definite. First, note that the covariance matrix can be expressed as

$$\mathrm{\Sigma }={\int }_{{\mathbb{R}}^n}\rho (\mathbf{x})(\mathbf{x}-\mathit{\mu })(\mathbf{x}-\mathit{\mu }{)}^{\mathrm{\top }}d\mathbf{x}$$

where $\rho (\mathbf{x})$ is the density of $\mathbf{X}$, and

$$\mathit{\mu }=\mathbb{E}[\mathbf{X}]={\int }_{{\mathbb{R}}^n}\rho (\mathbf{x})\mathbf{x}d\mathbf{x}\in {\mathbb{R}}^n$$

is the expectation value of $\mathbf{X}$. 

Let $\mathbf{v}\in {\mathbb{R}}^n$ be an arbitrary vector. We have

$$\begin{array}{rcl}{\mathbf{v}}^{\mathrm{\top }}\mathrm{\Sigma }\mathbf{v}& =& {\int }_{{\mathbb{R}}^n}\rho (\mathbf{x}){\mathbf{v}}^{\mathrm{\top }}(\mathbf{x}-\mathit{\mu })(\mathbf{x}-\mathit{\mu }{)}^{\mathrm{\top }}\mathbf{v}d\mathbf{x}\\ & =& {\int }_{{\mathbb{R}}^n}\rho (\mathbf{x})[(\mathbf{x}-\mathit{\mu }{)}^{\mathrm{\top }}\mathbf{v}{]}^{\mathrm{\top }}[(\mathbf{x}-\mathit{\mu }{)}^{\mathrm{\top }}\mathbf{v}]d\mathbf{x}\\ & =& {\int }_{{\mathbb{R}}^n}\rho (\mathbf{x})[(\mathbf{x}-\mathit{\mu }{)}^{\mathrm{\top }}\mathbf{v}{]}^{2}d\mathbf{x}\ge 0\end{array}$$

because the integrand is non-negative.

Remark (2022-10-14): In the initial version of this post, I used the eigenvalue decomposition of the covariance matrix to prove the positive semi-definiteness. But that is not only unnecessary but overly complicated. The above proof is simpler and more direct.