Covariance

Covariance Matrix

Earlier, we introduced the variance as a measure of how a single random variable deviates from its mean. However, real-world phenomena rarely involve just one variable. In machine learning, a data point is typically a vector of features, and understanding how these features relate to one another is essential for tasks ranging from dimensionality reduction to generative modeling. This motivates the transition from univariate to multivariate statistics, which begins with a systematic way to measure how pairs of random variables vary together.

Here, \((X - \mathbb{E}[X])\) and \((Y - \mathbb{E}[Y])\) are mean deviations, representing how far \(X\) and \(Y\) deviate from their respective means. The covariance measures how these deviations vary together (joint variation).

To work with two random variables together, we need a way to describe their joint behavior. The individual distributions of \(X\) and \(Y\), each described by its own pmf or pdf, are not enough. Two random variables may have identical individual distributions while being either independent or perfectly correlated. The required object is the joint distribution, which packages how \(X\) and \(Y\) take values together.

The one-variable statement is the Expectation of a Function of a Random Variable established earlier. In the discrete case the two-variable form follows by the same regrouping argument, collecting the pairs \((x, y)\) on which \(g\) takes each value and summing the joint mass over each group. The continuous case rests on the change-of-variables rule for joint densities, which we take for granted here and which belongs to the treatment of product measures.

Before deriving the standard computational identity for covariance, we prove the two-variable case of linearity of expectation. The proof given earlier covered a single random variable scaled and shifted by constants.

With bivariate linearity in hand, we obtain a computational identity for covariance that is often more convenient than the definition.

Independence of random variables has a powerful consequence for expectations of products. The expectation factors across the two variables.

The converse is not true in general. Zero covariance indicates only the absence of linear dependence, not full independence. A standard counter-example is \(X \sim \mathcal{N}(0, 1)\) and \(Y = X^2\). The standard normal density \(\varphi(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2}\) is an even function (\(\varphi(-x) = \varphi(x)\)), and \(x^3\) is an odd function, so \(\mathbb{E}[X^3] = \int_{\mathbb{R}} x^3 \varphi(x)\, dx = 0\) by odd-times-even symmetry on \(\mathbb{R}\). Together with \(\mathbb{E}[X] = 0\), this gives \(\operatorname{Cov}[X, Y] = \mathbb{E}[X^3] - \mathbb{E}[X]\,\mathbb{E}[X^2] = 0\). Yet \(Y\) is fully determined by \(X\), so \(X\) and \(Y\) are far from independent.

A positive (resp. negative) covariance indicates a positive (resp. negative) linear association between \(X\) and \(Y\). The absence of such linear association is what zero covariance certifies.

In practice, we extend this idea to \(n\) random variables. Consider a vector \(\boldsymbol{x} \in \mathbb{R}^n\) whose entries represent random variables \(X_1, \ldots, X_n\).

Each diagonal entry \(\Sigma_{ii}\) is the variance of \(X_i\), since \[ \operatorname{Cov}[X_i, X_i] = \mathbb{E}[(X_i - \mathbb{E}[X_i])^2] = \sigma_i^2. \] The total variance is defined as the trace of \(\Sigma\): \[ \operatorname{tr}(\Sigma) = \sum_{i=1}^n \operatorname{Var}[X_i]. \]

The covariance matrix has several important structural properties that connect probability theory to linear algebra.

First, \(\Sigma\) is symmetric because covariance itself is symmetric: \[ \Sigma_{ij} = \operatorname{Cov}[X_i, X_j] = \operatorname{Cov}[X_j, X_i] = \Sigma_{ji}. \] The same fact can be read at the matrix level. For every realization of the random vector, the outer product \((\boldsymbol{x} - \boldsymbol{\mu})(\boldsymbol{x} - \boldsymbol{\mu})^\top\) is a symmetric matrix, and the expectation of a symmetric-matrix-valued random variable is again symmetric (since \(\mathbb{E}\) is taken entrywise and commutes with transposition).

Since \(\Sigma\) is symmetric, by the spectral theorem, it is orthogonally diagonalizable: \[ \Sigma = P D P^\top \] where \(P\) is an orthogonal matrix whose columns are unit eigenvectors of \(\Sigma\), and \(D\) is a diagonal matrix of the corresponding eigenvalues.

In particular, the total variance equals the sum of the eigenvalues. Write \(P_{ik}\) for the entries of \(P\) and \(\lambda_1, \ldots, \lambda_n\) for the diagonal entries of \(D\). The \(i\)-th diagonal entry of \(\Sigma = P D P^\top\) is \[ \begin{align*} \Sigma_{ii} &= \sum_{k=1}^n P_{ik} \lambda_k (P^\top)_{ki} \\\\ &= \sum_{k=1}^n \lambda_k P_{ik}^2, \end{align*} \] so \[ \begin{align*} \operatorname{tr}(\Sigma) &= \sum_{i=1}^n \operatorname{Var}[X_i] \\\\ &= \sum_{i=1}^n \sum_{k=1}^n \lambda_k P_{ik}^2 \\\\ &= \sum_{k=1}^n \lambda_k \sum_{i=1}^n P_{ik}^2 \\\\ &= \sum_{k=1}^n \lambda_k, \end{align*} \] where the last step uses that each column of \(P\) is a unit vector. Note that the individual eigenvalues need not coincide with the individual variances \(\operatorname{Var}[X_i]\). In general only the totals agree.

As a consequence, the eigenvalues of \(\Sigma\) are non-negative. Permuting the columns of \(P\) and the diagonal entries of \(D\) together, we may order them so that \[ \lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_n \geq 0. \] The eigenvalue ordering is where Principal Component Analysis begins.

The matrix \(\mathbb{E}[\boldsymbol{x}\boldsymbol{x}^\top]\) is called the autocorrelation matrix, denoted \(R_{\boldsymbol{x}\boldsymbol{x}}\): \[ \begin{align*} R_{\boldsymbol{x}\boldsymbol{x}} &= \mathbb{E}[\boldsymbol{x}\boldsymbol{x}^\top] \\\\ &= \begin{bmatrix} \mathbb{E}[X_1^2] & \mathbb{E}[X_1 X_2] & \cdots & \mathbb{E}[X_1 X_n] \\\\ \mathbb{E}[X_2 X_1] & \mathbb{E}[X_2^2] & \cdots & \mathbb{E}[X_2 X_n] \\\\ \vdots & \vdots & \ddots & \vdots \\\\ \mathbb{E}[X_n X_1] & \mathbb{E}[X_n X_2] & \cdots & \mathbb{E}[X_n^2] \end{bmatrix}. \end{align*} \]

Principal Component Analysis (PCA)

The previous section established that every covariance matrix admits an orthogonal eigendecomposition with non-negative eigenvalues. This spectral structure has a direct and powerful interpretation. The eigenvectors identify the directions along which the data varies the most, and the eigenvalues quantify how much variance lies along each direction. Exploiting this insight leads to Principal Component Analysis (PCA), one of the most widely used techniques in data science for dimensionality reduction.

Before describing PCA, we establish the precise sense in which "eigenvalue equals variance along its eigenvector." This is what makes the spectral decomposition of \(\Sigma\) statistically meaningful in the first place.

This proposition has an immediate optimization consequence. Expanding any unit vector \(\boldsymbol{v}\) in the orthonormal eigenbasis \(\{\boldsymbol{p}_1, \ldots, \boldsymbol{p}_n\}\) as \(\boldsymbol{v} = \sum_i c_i \boldsymbol{p}_i\) with \(\sum_i c_i^2 = 1\), a direct computation using \(\Sigma \boldsymbol{p}_i = \lambda_i \boldsymbol{p}_i\) and orthonormality yields \(\boldsymbol{v}^\top \Sigma\, \boldsymbol{v} = \sum_i c_i^2 \lambda_i\), a convex combination of the eigenvalues. With the eigenvalues ordered as above, \(\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_n \geq 0\), this combination is maximized by placing all weight on the largest eigenvalue: \[ \begin{align*} \max_{\|\boldsymbol{v}\|=1} \operatorname{Var}[\boldsymbol{v}^\top \boldsymbol{x}] &= \max_{\|\boldsymbol{v}\|=1} \boldsymbol{v}^\top \Sigma\, \boldsymbol{v} \\\\ &= \lambda_1, \end{align*} \] attained at \(\boldsymbol{v} = \boldsymbol{p}_1\). The variance-maximizing directions are exactly the unit eigenvectors for the largest eigenvalue.

The orthogonal diagonalization of the covariance matrix \(\Sigma\) is important in statistics and machine learning for analyzing data and reducing dimensionality: \[ \Sigma = P D P^\top. \] Here, the column vectors of \(P\) (unit eigenvectors of \(\Sigma\)) are called the principal components (PCs) of the data, and the diagonal entries of \(D\) (eigenvalues of \(\Sigma\)) represent the variances along these principal components. Each eigenvalue indicates "how much" of the total variance of \(\Sigma\) is captured by its corresponding principal component.

The first principal component (PC1) corresponds to the largest eigenvalue and represents the direction in the data space where the distribution varies the most (the direction maximizing the variance of the data). The second principal component (PC2) captures the next largest variance and is orthogonal to PC1. Similarly, each subsequent PC captures no more variance than the one before it and maintains orthogonality to all previous PCs.

For example, suppose we have an observed data matrix \(X \in \mathbb{R}^{m \times 10}\), which contains \(m\) data points, each represented as a vector \(\boldsymbol{x}_i \in \mathbb{R}^{10}\). We compute the sample covariance matrix and diagonalize it: \[ S = \frac{1}{m-1}(X-\bar{X})^\top (X-\bar{X}) = PDP^\top, \] where \(\bar{X} \in \mathbb{R}^{m \times 10}\) is the matrix whose every row equals the column-mean vector \(\bar{\boldsymbol{x}} = \frac{1}{m}\sum_{i=1}^m \boldsymbol{x}_i\). We then examine the eigenvalues to determine how much variance each principal component captures, comparing with the total variance of \(S\).

Suppose PC1 captures 40% of the total variance, PC2 30%, and PC3 25%. We project the data onto the subspace spanned by the first 3 most significant principal components. Even though we discard the remaining 7 dimensions (which carry only a small fraction of the total variance, typically interpreted as noise or fine-scale variation), we still retain the most significant patterns (trends) in the data. \[ \boldsymbol{x}_i \in \mathbb{R}^{10} \to \boldsymbol{z}_i \in \mathbb{R}^3. \] The vector \(\boldsymbol{z}_i\) is known as the latent vector.

This dimensionality reduction process is called principal component analysis (PCA). By reducing dimensions, PCA enables efficient analysis of large datasets while retaining their most meaningful structure. PCA is widely used in machine learning, image processing, and exploratory data analysis for dimensionality reduction and noise filtering.

PCA with Singular Value Decomposition (SVD)

In the derivation above, PCA was formulated through the eigendecomposition of the sample covariance matrix \(S = PDP^\top\). While conceptually clean, this approach requires explicitly forming \(X^\top X\), which can be expensive for high-dimensional data and may suffer from numerical instability. The singular value decomposition (SVD) provides an alternative that avoids these issues while yielding the same principal components.

Given a data matrix in mean-deviation form \(X \in \mathbb{R}^{m \times n}\) (so that \(\bar{X} = 0\), where \(m\) is the number of data points and \(n\) is the number of features), the sample covariance matrix simplifies to \[ S = \frac{1}{m-1} X^\top X. \] Note: in PCA, it is common to work with \(X^\top X \in \mathbb{R}^{n \times n}\) rather than \(XX^\top \in \mathbb{R}^{m \times m}\) because we are typically more interested in understanding the relationships among the \(n\) features than among the \(m\) data points.

Now apply the SVD to \(X\). We write \(X = U \Sigma V^\top\), where \(U \in \mathbb{R}^{m \times m}\) and \(V \in \mathbb{R}^{n \times n}\) are orthogonal and \(\Sigma \in \mathbb{R}^{m \times n}\) is the rectangular diagonal matrix of singular values \(\sigma_1 \geq \sigma_2 \geq \cdots \geq 0\). (Here \(\Sigma\) denotes the SVD's singular-value matrix in the standard linear-algebra convention. The population covariance matrix, written with the same letter in the earlier sections, is meant only where it is named as such.) Substituting, \[ \begin{align*} X^\top X &= (U\Sigma V^\top)^\top (U\Sigma V^\top) \\\\ &= V\Sigma^\top U^\top U\Sigma V^\top \\\\ &= V\,(\Sigma^\top \Sigma)\,V^\top, \end{align*} \] using \(U^\top U = I_m\) (orthogonality of \(U\)).

The product \(\Sigma^\top \Sigma\) is the \(n \times n\) square diagonal matrix \(\operatorname{diag}(\sigma_1^2, \ldots, \sigma_n^2)\), where \(\sigma_i = 0\) for \(i \gt \min(m, n)\). Dividing by \(m-1\) yields the sample covariance matrix in eigendecomposition form: \[ \begin{align*} S &= \frac{1}{m-1} X^\top X \\\\ &= V \cdot \operatorname{diag}\!\left(\frac{\sigma_1^2}{m-1}, \ldots, \frac{\sigma_n^2}{m-1}\right) \cdot V^\top. \end{align*} \]

Comparing with the eigendecomposition \(S = PDP^\top\) from the previous section, and writing \(\boldsymbol{v}_i\) for the \(i\)-th column of \(V\) (the \(i\)-th right singular vector), we may take \[ \begin{align*} P &= V, \\\\ \boldsymbol{p}_i &= \boldsymbol{v}_i, \\\\ \lambda_i &= \frac{\sigma_i^2}{m-1}. \end{align*} \] The right singular vectors can therefore be taken as the principal components, and the eigenvalues of \(S\) are the squared singular values divided by \(m-1\).

To interpret these eigenvalues as variances, recall that the Variance Along a Direction proposition was stated for a random vector \(\boldsymbol{x}\) with population covariance \(\Sigma\). The proposition was proved for an arbitrary random vector, so it applies in particular to the empirical distribution that places mass \(1/m\) on each observed data point. That distribution has mean \(\mathbf{0}\) and covariance matrix \(X^\top X / m\), so \(\boldsymbol{v}^\top (X^\top X / m)\, \boldsymbol{v}\) is the variance of the projections onto \(\boldsymbol{v}\). The sample covariance matrix rescales this by \(m/(m-1)\), so \(\boldsymbol{v}^\top S \boldsymbol{v}\) is the Bessel-corrected estimate of that variance, the unbiased sample variance.

Thus the eigenvalues of \(S\) are precisely the sample variances along the principal components: \[ \begin{align*} \widehat{\operatorname{Var}}\!\left[\boldsymbol{v}_i^\top \boldsymbol{x}\right] &= \frac{\sigma_i^2}{m-1}, \\\\ \text{sample standard deviation along PC } i &= \frac{\sigma_i}{\sqrt{m-1}}. \end{align*} \]

SVD offers several numerical advantages for PCA. For high-dimensional data, computing \(X^\top X\) explicitly would be expensive or memory-intensive, but SVD can be computed directly from \(X\). This gives access to the principal components and variances without ever forming \(X^\top X\).

                        import numpy as np

                        # Random data matrix ( m data points with n features)
                        def generate_data(m, n):
                            data = np.random.randn(m, n) 

                            # Make some correlations 
                            data[:, 2] = 0.7 * data[:, 0] + 0.5 * data[:, 2]
                            data[:, 3] = 0.2 * data[:, 0] + 0.5 * data[:, 1] 
                            data[:, 4] = -0.3 * data[:, 1] + 0.2 * data[:, 2]
                            data[:, 5] = 0.4 * data[:, 0] + 0.1 * data[:, 1]
                            data[:, 6] = 0.8 * data[:, 3] + -0.3 * data[:, 2]

                            # Mean-deviation form 
                            return (data - np.mean(data, axis=0) )

                        # PCA with covariance matrix (As reference)
                        def pca_via_covariance(data):
                            
                            # "Sample" covariance matrix 
                            # Note: Dividing by (m-1) provides "unbiased" estimate of the population covariance. 
                            cov_matrix = np.dot(data.T, data) / (data.shape[0] - 1)

                            # Eigendecomposition 
                            # np.linalg.eigh() is for symmetric matrix (real, diagonaizable), better than np.linalg.eig() 
                            eigvals, eigvecs = np.linalg.eigh(cov_matrix) 

                            #  Make eigenvalues & eigenvectors in descending order
                            idx = np.argsort(eigvals)[::-1] 
                            eigvals = eigvals[idx]
                            eigvecs = eigvecs[:, idx] 

                            # Each variance vs total variance
                            ratio = eigvals / np.sum(eigvals)

                            return eigvals, eigvecs, ratio

                        # PCA with SVD (we use this function)
                        def pca_with_svd(data):
                            
                            # Singular Value Decomposition (we don't need the matrix U: use "_" )
                            _, S, vt = np.linalg.svd(data, full_matrices=False)
                            
                            # Get eigenvalues via singular values: lambda_i = (S_i)^2  / m - 1
                            eigvals = (S ** 2) / (data.shape[0] - 1)
                            
                            # Each variance vs total variance
                            ratio = eigvals / np.sum(eigvals)
                            
                            return eigvals, vt.T, ratio 

                        # Threshold t for the number of PCs
                        def threshold(var_ratio, t):
                            
                            # Compute cumulative variance
                            cum_variance = np.cumsum(var_ratio)
                            
                            # Find the number of components for retaining a fraction t of the variance
                            num_pcs = np.argmax(cum_variance >= t) + 1
                            return num_pcs

                        # Dimension of data 
                        m = 1000000 
                        n = 10 

                        # Run PCA
                        eigvals, eigvecs, ratio = pca_with_svd(generate_data(m, n))

                        # Threshold for variance retention
                        t = 0.95
                        num_pcs = threshold(ratio, t)

                        # Print results
                        print("Eigenvalues:")
                        for i, val in enumerate(eigvals):
                            print(f"  Lambda {i + 1}: {val:.6f}")

                        print("\nExplained Variance Ratio (%):")
                        for i, var in enumerate(ratio):
                            print(f"  PC{i + 1}: {var * 100:.2f}%")

                        print(f"\nTo retain {t * 100:.0f}% of variance, use the first {num_pcs} PCs.")
                            

While covariance captures the magnitude and direction of linear association, its values depend on the scales of the variables involved. In the next page, we introduce the correlation coefficient and correlation matrix, which standardize covariance to produce "scale-free" measures of linear relationship.