Sampling Multivariate Normals via Spectral Decomposition
You need to generate samples from a multivariate normal distribution $X \sim N(0, \Sigma)$ where $\Sigma$ is an $N \times N$ positive semi-definite covariance matrix. The standard approach is Cholesky decomposition, but suppose Cholesky is off the table -- maybe $\Sigma$ is singular or nearly singular, or you just want to understand the alternative.
- Explain how to use the spectral (eigenvalue) decomposition of $\Sigma$ to generate these samples. Give the step-by-step procedure.
- Prove that your method produces the correct covariance structure, i.e., that $\text{Cov}(X) = \Sigma$.
- When would you prefer spectral decomposition over Cholesky in practice?
Hints
- Any factorization $\Sigma = A A^\top$ lets you transform independent standard normals into correlated ones with the right covariance. What factorization does eigendecomposition give you?
- Write $\Sigma = V \Lambda V^\top$ and think about what $V \Lambda^{1/2}$ does geometrically -- it scales along the principal axes and then rotates.
- Set $X = V \Lambda^{1/2} Z$ where $Z \sim N(0, I)$, and verify by computing $E[X X^\top]$. For the comparison with Cholesky, think about what happens when $\Sigma$ has zero eigenvalues.
Worked Solution
How to Think About It: The core idea is simple: if you can write $\Sigma = A A^\top$ for any matrix $A$, then $X = A Z$ where $Z \sim N(0, I)$ gives you $X \sim N(0, \Sigma)$. Cholesky gives you one such $A$ (lower triangular), but spectral decomposition gives you another one that is sometimes better. The spectral version factors $\Sigma$ through its eigenvectors and eigenvalues, which has a nice geometric interpretation: you are generating independent noise along the principal axes of the ellipsoid defined by $\Sigma$, then rotating into the original coordinate system.
Key Insight: Any symmetric PSD matrix $\Sigma$ can be written as $\Sigma = V \Lambda V^\top$ where $V$ is orthogonal and $\Lambda$ is diagonal with non-negative entries. Taking elementwise square roots of $\Lambda$ gives you a "square root" of $\Sigma$: set $A = V \Lambda^{1/2}$, and $A A^\top = V \Lambda^{1/2} (V \Lambda^{1/2})^\top = V \Lambda V^\top = \Sigma$.
The Method:
- Compute the eigendecomposition $\Sigma = V \Lambda V^\top$, where $V = [v_1 \mid \cdots \mid v_N]$ is the matrix of eigenvectors and $\Lambda = \text{diag}(\lambda_1, \ldots, \lambda_N)$.
- Form the matrix square root: $\Lambda^{1/2} = \text{diag}(\sqrt{\lambda_1}, \ldots, \sqrt{\lambda_N})$. If any $\lambda_i$ is slightly negative due to numerical noise, clip it to zero.
- Generate $Z \sim N(0, I_N)$ -- a vector of $N$ independent standard normals.
- Compute $X = V \Lambda^{1/2} Z$.
Proof of Correctness:
$$\text{Cov}(X) = E[X X^\top] = V \Lambda^{1/2} E[Z Z^\top] \Lambda^{1/2} V^\top = V \Lambda^{1/2} I \Lambda^{1/2} V^\top = V \Lambda V^\top = \Sigma$$
Since $X$ is a linear transformation of a Gaussian vector $Z$, $X$ is also Gaussian, so matching the mean (zero) and covariance ($\Sigma$) fully determines the distribution.
Implementation:
```python import numpy as np
def sample_mvn_spectral(Sigma, n_samples=1): eigenvalues, V = np.linalg.eigh(Sigma) eigenvalues = np.maximum(eigenvalues, 0) # clip numerical noise A = V * np.sqrt(eigenvalues) # broadcasting: V @ diag(sqrt(eig)) Z = np.random.randn(len(Sigma), n_samples) return (A @ Z).T # shape: (n_samples, N) ```
Practical Considerations:
- Singular covariance: Cholesky requires $\Sigma$ to be strictly positive definite. Spectral decomposition handles singular $\Sigma$ gracefully -- zero eigenvalues just mean you don't generate noise along those directions.
- Near-singular matrices: Eigendecomposition is more numerically stable when $\Sigma$ has a large condition number.
- Speed: Cholesky is faster, $O(N^3/3)$ vs $O(N^3)$ for full eigendecomposition. If $\Sigma$ is well-conditioned and you just need samples, Cholesky wins on speed.
- Dimensionality reduction: If only $k \ll N$ eigenvalues are meaningfully large, you can truncate: use only the top $k$ eigenvectors and generate $k$-dimensional $Z$. This gives you approximate sampling in $O(kN)$ per sample, which is useful in high-dimensional settings like PCA-based risk models.
Answer: Generate $Z \sim N(0, I_N)$ and compute $X = V \Lambda^{1/2} Z$ where $\Sigma = V \Lambda V^\top$ is the eigendecomposition. This works because $\text{Cov}(X) = V \Lambda V^\top = \Sigma$. Prefer spectral over Cholesky when $\Sigma$ is singular, nearly singular, or when you want to exploit low-rank structure for dimensionality reduction.
Intuition
The fundamental trick behind all multivariate normal sampling is the same: find any matrix $A$ such that $A A^\top = \Sigma$, then multiply independent standard normals by $A$. Cholesky and spectral decomposition are just two different ways to factor $\Sigma$. The spectral version is geometrically transparent -- the eigenvectors define the axes of the covariance ellipsoid, the eigenvalues define how stretched each axis is, and sampling amounts to generating noise along each axis independently and rotating into the original frame.
This shows up constantly in quant work. Risk models (like Barra or Axioma) often deliver covariance matrices that are singular or nearly singular because the number of factors is smaller than the number of assets. Cholesky will choke on these; spectral decomposition handles them naturally and also lets you do truncated sampling -- keep only the top $k$ eigenvectors and you get fast, approximate simulation that captures most of the variance. This is the backbone of scenario generation and Monte Carlo in production risk systems.