While the element-wise definition ($cov(x_i, x_j)$) is useful for understanding the relationship between two specific variables, we rarely calculate the matrix element-by-element in practice. Instead, we use vector algebra to compute the entire matrix at once.
Given a set of $m$ data points (samples) $\{\mathbf{x}^{(1)}, \mathbf{x}^{(2)}, \dots, \mathbf{x}^{(m)}\}$, where each $\mathbf{x}^{(i)} \in \mathbb{R}^n$, and the empirical mean vector:$$\mathbf{\bar{x}} = \frac{1}{m} \sum_{i=1}^{m} \mathbf{x}^{(i)}$$The Empirical Covariance Matrix $\mathbf{C}$ is defined as the average of the "spread" of the data around the mean. Writing the centered sample as $\mathbf{\tilde x}^{(i)} = \mathbf{x}^{(i)} - \mathbf{\bar{x}}$, this is:$$\mathbf{C} = \frac{1}{m-1} \sum_{i=1}^{m} \mathbf{\tilde x}^{(i)} (\mathbf{\tilde x}^{(i)})^T$$
To see why this holds, look at the $(j, k)$-th entry of the matrix resulting from the sum above. Recall the centered sample $\mathbf{\tilde x}^{(i)} = \mathbf{x}^{(i)} - \mathbf{\bar{x}}$.
The product of a column vector and a row vector—the outer product—results in a matrix:$$\mathbf{\tilde x}^{(i)} (\mathbf{\tilde x}^{(i)})^T = \begin{pmatrix} \tilde x_1^{(i)} \\ \vdots \\ \tilde x_n^{(i)} \end{pmatrix} \begin{pmatrix} \tilde x_1^{(i)} & \dots & \tilde x_n^{(i)} \end{pmatrix} = \begin{pmatrix} \tilde x_1^{(i)}\tilde x_1^{(i)} & \dots & \tilde x_1^{(i)}\tilde x_n^{(i)} \\ \vdots & \ddots & \vdots \\ \tilde x_n^{(i)}\tilde x_1^{(i)} & \dots & \tilde x_n^{(i)}\tilde x_n^{(i)} \end{pmatrix}$$
Comparing this to your previous definition:
The outer product $\mathbf{\tilde x}^{(i)}(\mathbf{\tilde x}^{(i)})^T$ represents the contribution to the variance by a single sample.
If we organize our $m$ measurements (centered by subtracting the mean $\bar{\mathbf{x}}$) as rows in an $m \times n$ matrix $\tilde{\mathbf{X}}$:
$$\tilde{\mathbf{X}} = \begin{pmatrix} — & (\mathbf{\tilde x}^{(1)})^T & — \\ — & (\mathbf{\tilde x}^{(2)})^T & — \\ & \vdots & \\ — & (\mathbf{\tilde x}^{(m)})^T & — \end{pmatrix}$$
The covariance matrix $\mathbf{C}$ (which must be $n \times n$) is then:
$$\mathbf{C} = \frac{1}{m-1} \tilde{\mathbf{X}}^T \tilde{\mathbf{X}}$$
Why this inner product reproduces the covariance. Writing the product out makes it explicit. The transpose $\tilde{\mathbf{X}}^T$ is an $n \times m$ matrix whose rows are the components and whose columns are the samples, while $\tilde{\mathbf{X}}$ is $m \times n$ with the samples as rows:
$$ \mathbf{C} = \frac{1}{m-1} \tilde{\mathbf{X}}^T \tilde{\mathbf{X}} = \frac{1}{m-1} \begin{pmatrix} \tilde x_1^{(1)} & \tilde x_1^{(2)} & \tilde x_1^{(3)} & \dots \\ \tilde x_2^{(1)} & \tilde x_2^{(2)} & \tilde x_2^{(3)} & \dots \\ \vdots & \vdots & \vdots & \\ \end{pmatrix} \begin{pmatrix} \tilde x_1^{(1)} & \tilde x_2^{(1)} & \dots \\ \tilde x_1^{(2)} & \tilde x_2^{(2)} & \dots \\ \tilde x_1^{(3)} & \tilde x_2^{(3)} & \dots \\ \vdots & \vdots & \\ \end{pmatrix} = \frac{1}{m-1} \begin{pmatrix} \sum_i \tilde x_1^{(i)} \tilde x_1^{(i)} & \sum_i \tilde x_1^{(i)} \tilde x_2^{(i)} & \dots \\ \sum_i \tilde x_2^{(i)} \tilde x_1^{(i)} & \sum_i \tilde x_2^{(i)} \tilde x_2^{(i)} & \dots \\ \vdots & \vdots & \ddots \\ \end{pmatrix} $$
The row-times-column rule makes the $(j,k)$-th entry a sum over the $m$ samples,
$$C_{jk} = \frac{1}{m-1} \sum_{i=1}^{m} \tilde x_j^{(i)} \tilde x_k^{(i)},$$
which is exactly the empirical covariance $cov(x_j, x_k)$ from the outer-product form above. So the single matrix product $\tilde{\mathbf{X}}^T \tilde{\mathbf{X}}$ carries out all $n \times n$ of those sample-sums at once — that is why this form is the workhorse in code.
Alternatively, we can organize in a centered data matrix $\hat{\mathbf{X}}$ each data vector as a column, i.e. $\hat{\mathbf{X}}$ is an $n \times m$-matrix:
$$\mathbf{C} = \frac{1}{m-1} \hat{\mathbf{X}} \hat{\mathbf{X}}^T$$
The empirical covariance can be expressed at different levels of abstraction. Each form is mathematically equivalent but serves a different purpose in implementation and theory.
| Level | Formula | Best used for... |
|---|---|---|
| Scalar | $cov(x_1, x_2) = \frac{1}{m-1} \sum_{i=1}^{m} (x_1^{(i)} - \bar{x}_1)(x_2^{(i)} - \bar{x}_2)$ | Understanding the statistical relationship between two specific variables. |
| Outer Product | $\mathbf{C} = \frac{1}{m-1} \sum_{i=1}^{m} \mathbf{\tilde x}^{(i)}(\mathbf{\tilde x}^{(i)})^T$ | Theoretical derivations, Eigen-analysis, and Rank-1/Rank-$\mu$ updates in CMA-ES. |
| Data Matrix | $\mathbf{C} = \frac{1}{m-1} \tilde{\mathbf{X}}^T\tilde{\mathbf{X}}$ | High-performance computation and software implementations (assuming samples as rows). |
Note on Notation: In the Data Matrix form above, $\tilde{\mathbf{X}}$ is an $m \times n$ matrix where each row is a centered observation vector $\mathbf{\tilde x}^{(i)} = \mathbf{x}^{(i)} - \bar{\mathbf{x}}$.