Multivariate Normal Distribution

Learning Objectives

The multivariate normal (MVN) distribution is the most important continuous distribution in statistics and machine learning. By the end of this section, you will master the MVN and understand why it appears everywhere in modern AI systems.

1. Define the MVN mathematically and understand every component of its PDF: the mean vector $\boldsymbol{\mu}$, the covariance matrix $\boldsymbol{\Sigma}$, and the normalizing constant
2. Visualize and interpret the elliptical contours of the MVN, understanding how eigenvalues and eigenvectors determine the shape and orientation
3. Compute marginal and conditional distributions of MVN, recognizing that they remain Gaussian with closed-form parameters
4. Understand the affine transformation property: linear transformations of MVN are still MVN, which is fundamental for neural network analysis
5. Appreciate high-dimensional behavior: the concentration of measure phenomenon and its implications for machine learning
6. Apply MVN concepts in practical contexts: Gaussian processes, VAEs, Bayesian deep learning, and probabilistic PCA
7. Implement MVN sampling and density evaluation efficiently using Cholesky decomposition and understand the Mahalanobis distance

The Big Picture: Why MVN Matters

"The multivariate normal is the hydrogen atom of statistics—simple enough to analyze completely, yet rich enough to capture the essential phenomena."

The multivariate normal extends the familiar bell curve to $d$ dimensions. Just as the univariate normal $N(\mu, \sigma^2)$ is parameterized by a mean and variance, the MVN $N(\boldsymbol{\mu}, \boldsymbol{\Sigma})$ is parameterized by:

• Mean vector $\boldsymbol{\mu} \in \mathbb{R}^d$: The center of the distribution, the expected value $E[\mathbf{X}]$
• Covariance matrix $\boldsymbol{\Sigma} \in \mathbb{R}^{d \times d}$: A symmetric positive definite matrix that encodes variances and correlations

The covariance matrix is the key to understanding MVN geometry. Its entries tell us:

• Diagonal entries $\Sigma_{ii}$: The variance of each component
• Off-diagonal entries $\Sigma_{ij}$: Covariances between pairs of components

Historical Origins: From Gauss to Modern AI

Carl Friedrich Gauss (1777-1855)

The normal distribution is often called the "Gaussian distribution" after Carl Friedrich Gauss, who derived it in the context of astronomical measurement errors around 1809. Gauss showed that if errors arise from the sum of many small independent effects, their distribution approaches the normal—the first statement of what we now call the Central Limit Theorem.

Francis Galton (1822-1911) and Karl Pearson (1857-1936)

Galton extended univariate ideas to the bivariate case, discovering regression to the mean and developing the concept of correlation. His student Karl Pearson formalized the bivariate normal and pioneered multivariate analysis, developing Principal Component Analysis (PCA).

R.A. Fisher and Modern Statistics

R.A. Fisher (1890-1962) developed the theoretical foundations for maximum likelihood estimation and multivariate analysis, showing how MVN structure enables powerful inference procedures. The Fisher information matrix—central to modern neural network optimization—derives its geometric meaning from the MVN.

Mathematical Definition

The Probability Density Function

A random vector $\mathbf{X} = (X_1, \ldots, X_d)^T$ follows a d-dimensional multivariate normal distribution, written$\mathbf{X} \sim N(\boldsymbol{\mu}, \boldsymbol{\Sigma})$, if its probability density function is:

Let's understand each component:

The Mahalanobis Distance

The exponent contains the squared Mahalanobis distance:

This is the "natural" distance for MVN—it measures how many standard deviations a point is from the mean, accounting for correlations. Points on the same contour ellipse have the same Mahalanobis distance.

Interactive PDF Explorer

Explore how the MVN parameters affect the distribution. Adjust the mean vector, variances, and correlation to see how the contour ellipses change. Toggle the principal axes to see the eigenvector directions, and enable marginal distributions to see the projections.

Geometry: Ellipses and Eigenvalues

The shape of MVN contours is determined by the eigendecomposition of the covariance matrix:

where $\mathbf{Q}$ contains the eigenvectors as columns, and$\boldsymbol{\Lambda}$ is diagonal with eigenvalues$\lambda_1, \ldots, \lambda_d$.

• Eigenvectors give the directions of the principal axes
• Eigenvalues give the variances along each principal axis
• √λ gives the standard deviation along each axis

Key Formulas for Bivariate Case

For the 2D case with covariance matrix:

The eigenvalues are:

And the principal axis makes angle θ with the x-axis:

Key Properties of MVN

1. Marginal Distributions

Every marginal distribution of an MVN is also normal:

More generally, any subset of components has a multivariate normal distribution with the corresponding subvector of means and submatrix of covariances.

2. Affine Transformation Property

This is perhaps the most powerful property. If $\mathbf{X} \sim N(\boldsymbol{\mu}, \boldsymbol{\Sigma})$and $\mathbf{A}$ is a matrix with $\mathbf{b}$ a vector, then:

3. Sum of Independent Gaussians

If $\mathbf{X} \sim N(\boldsymbol{\mu}_X, \boldsymbol{\Sigma}_X)$ and$\mathbf{Y} \sim N(\boldsymbol{\mu}_Y, \boldsymbol{\Sigma}_Y)$ are independent:

Conditional Distributions

One of the most beautiful properties of MVN: conditionals are also Gaussian, with closed-form parameters. Partition $\mathbf{X}$ as:

Then the conditional distribution is:

with:

High Dimensions: The Concentration of Measure

In machine learning, we often work with very high-dimensional data (images with millions of pixels, language models with billions of parameters). The MVN exhibits surprising behavior in high dimensions that is crucial to understand.

Key High-Dimensional Phenomena

1. Shell Concentration: For standard MVN in $d$ dimensions, almost all samples lie on a thin shell at distance $\sqrt{d}$ from the origin. The origin—the mode of the distribution—has near-zero probability!
2. Orthogonality: Random pairs of samples become nearly orthogonal as $d \to \infty$. The expected cosine between two random unit vectors goes to 0.
3. Curse of Dimensionality: The volume of the unit ball shrinks exponentially. Most of the volume concentrates in the "corners" of the hypercube.

AI/ML Applications

1. Variational Autoencoders (VAEs)

VAEs use MVN as both prior and approximate posterior:

The reparameterization trick uses the affine property: sample$\mathbf{z} = \boldsymbol{\mu} + \boldsymbol{\sigma} \odot \boldsymbol{\epsilon}$where $\boldsymbol{\epsilon} \sim N(\mathbf{0}, \mathbf{I})$.

2. Gaussian Processes

A GP defines a prior over functions where any finite collection of function values is MVN:

The conditional distribution formula gives us posterior predictions with uncertainty!

3. Bayesian Neural Networks

Weight distributions are often modeled as MVN:

The Laplace approximation approximates the posterior as Gaussian using the Hessian.

4. Probabilistic PCA

PPCA assumes latent factors are MVN:

Marginalizing over $\mathbf{z}$, we get $\mathbf{x} \sim N(\boldsymbol{\mu}, \mathbf{W}\mathbf{W}^T + \sigma^2\mathbf{I})$.

5. Fisher Information and Natural Gradient

For MVN, the Fisher information has special structure enabling efficient natural gradient descent. This underpins K-FAC and other second-order optimization methods for neural networks.

Python Implementation

Basic MVN Operations

1import numpy as np
2from scipy.stats import multivariate_normal
3import matplotlib.pyplot as plt
4
5# Define the multivariate normal distribution
6# Mean vector (2D)
7mu = np.array([1.0, 2.0])
8
9# Covariance matrix (2x2, symmetric positive definite)
10# Variances on diagonal, covariances off-diagonal
11sigma_x, sigma_y = 1.5, 1.0
12rho = 0.7  # Correlation coefficient
13cov = rho * sigma_x * sigma_y
14
15Sigma = np.array([
16    [sigma_x**2, cov],
17    [cov, sigma_y**2]
18])
19
20# Create the multivariate normal distribution
21mvn = multivariate_normal(mean=mu, cov=Sigma)
22
23# Sample from the distribution
24np.random.seed(42)
25samples = mvn.rvs(size=1000)
26
27# Evaluate PDF at a specific point
28point = np.array([1.5, 2.5])
29pdf_value = mvn.pdf(point)
30print(f"PDF at {point}: {pdf_value:.4f}")
31
32# Compute eigenvalues and eigenvectors of Sigma
33eigenvalues, eigenvectors = np.linalg.eig(Sigma)
34print(f"\nEigenvalues: {eigenvalues}")
35print(f"Eigenvectors:\n{eigenvectors}")
36
37# Principal axis angle
38angle = np.arctan2(eigenvectors[1, 0], eigenvectors[0, 0])
39print(f"Principal axis angle: {np.degrees(angle):.1f}°")
40
41# Verify: marginal distributions
42print(f"\nMarginal X: N({mu[0]}, {Sigma[0,0]:.3f})")
43print(f"Marginal Y: N({mu[1]}, {Sigma[1,1]:.3f})")
44
45# Conditional distribution Y | X = x_obs
46x_obs = 2.0
47cond_mean = mu[1] + rho * (sigma_y / sigma_x) * (x_obs - mu[0])
48cond_var = sigma_y**2 * (1 - rho**2)
49print(f"\nConditional Y|X={x_obs}: N({cond_mean:.3f}, {cond_var:.3f})")

Sampling via Cholesky Decomposition

1import numpy as np
2import matplotlib.pyplot as plt
3
4def sample_mvn_cholesky(mu, Sigma, n_samples):
5    """
6    Sample from MVN using Cholesky decomposition.
7    This is how most libraries implement MVN sampling.
8
9    X = μ + L @ Z, where Z ~ N(0, I) and L @ L.T = Σ
10    """
11    d = len(mu)
12
13    # Cholesky decomposition: Σ = L @ L.T
14    L = np.linalg.cholesky(Sigma)
15
16    # Sample standard normal vectors
17    Z = np.random.randn(n_samples, d)
18
19    # Transform: X = μ + L @ Z
20    X = mu + Z @ L.T
21
22    return X
23
24# Example: Bivariate normal
25mu = np.array([0, 0])
26Sigma = np.array([[1.0, 0.6],
27                  [0.6, 1.0]])
28
29# Generate samples
30np.random.seed(42)
31samples = sample_mvn_cholesky(mu, Sigma, 2000)
32
33# Verify statistics
34print("Sample mean:", samples.mean(axis=0))
35print("Sample covariance:\n", np.cov(samples.T))
36print("True covariance:\n", Sigma)
37
38# The Cholesky factor L
39L = np.linalg.cholesky(Sigma)
40print("\nCholesky factor L:")
41print(L)
42print("\nVerify L @ L.T = Σ:")
43print(L @ L.T)

Mahalanobis Distance and Confidence Ellipses

1import numpy as np
2from scipy.stats import chi2
3
4def mahalanobis_distance(x, mu, Sigma):
5    """
6    Compute Mahalanobis distance from x to μ.
7
8    D² = (x - μ)ᵀ Σ⁻¹ (x - μ)
9
10    This is the "natural" distance metric for MVN.
11    """
12    diff = x - mu
13    Sigma_inv = np.linalg.inv(Sigma)
14    return np.sqrt(diff @ Sigma_inv @ diff)
15
16# Example
17mu = np.array([0, 0])
18Sigma = np.array([[1.0, 0.8],
19                  [0.8, 1.0]])
20
21# Compare two points at same Euclidean distance
22point1 = np.array([1.5, 0.0])   # Along x-axis
23point2 = np.array([1.06, 1.06]) # Along principal axis
24
25# Euclidean distances
26eucl_1 = np.linalg.norm(point1 - mu)
27eucl_2 = np.linalg.norm(point2 - mu)
28
29# Mahalanobis distances
30maha_1 = mahalanobis_distance(point1, mu, Sigma)
31maha_2 = mahalanobis_distance(point2, mu, Sigma)
32
33print(f"Point 1: {point1}")
34print(f"  Euclidean: {eucl_1:.3f}, Mahalanobis: {maha_1:.3f}")
35print(f"\nPoint 2: {point2}")
36print(f"  Euclidean: {eucl_2:.3f}, Mahalanobis: {maha_2:.3f}")
37
38# For MVN: D² follows chi-squared distribution!
39# Probability that D² <= threshold
40d = 2  # dimensions
41threshold = 2.0
42prob = chi2.cdf(threshold**2, df=d)
43print(f"\nP(D² <= {threshold}²) = {prob:.3f}")
44print(f"This defines a {prob*100:.1f}% confidence ellipse")

Common Pitfalls

Test Your Understanding

Summary

The multivariate normal distribution is foundational to modern statistics and machine learning. Here are the key takeaways:

Mathematical Foundations

• PDF: $f(\mathbf{x}) = \frac{1}{(2\pi)^{d/2}|\Sigma|^{1/2}} \exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^T\Sigma^{-1}(\mathbf{x}-\boldsymbol{\mu})\right)$
• Parameters: Mean vector $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$
• Contours: Ellipsoids whose axes align with eigenvectors of $\boldsymbol{\Sigma}$

Key Properties

• Marginals are Gaussian: Any subset of components is MVN
• Conditionals are Gaussian: With linear mean and constant variance
• Affine transformations preserve Gaussianity: $\mathbf{AX} + \mathbf{b} \sim N(\mathbf{A}\boldsymbol{\mu}+\mathbf{b}, \mathbf{A}\boldsymbol{\Sigma}\mathbf{A}^T)$
• Uncorrelated = Independent: Only for MVN!

High-Dimensional Behavior

• Probability concentrates on a shell at distance $\sqrt{d}$ from origin
• Random vectors become nearly orthogonal
• This affects weight initialization, normalization, and optimization in neural networks

AI/ML Applications

• VAEs: MVN priors and posteriors with reparameterization
• Gaussian Processes: Prior over functions, exact posterior inference
• Bayesian Neural Networks: Gaussian weight posteriors
• Probabilistic PCA: Latent variable models
• Natural Gradient: Fisher information geometry