Why Multivariate Distribution with Covariance Matrix with Rank 1 Works?

Question

So I want to sample from multivariate normal distribution and have this code where mean is 0 and I added covariance matrix with all entries to 1 implying all random variables are equally correlated.

import numpy as np
X = [0,1,2]
samples = np.random.multivariate_normal([0,0,0], [[1,1,1],[1,1,1],[1,1,1]]) 
print(samples)

>> samples [-0.89635305 -0.89635305 -0.89635305]

The question is for computing the trinormal distribution the cholesky decomposition of the covariance matrix has to be done, but here the rank of the matrix is 1, so why the code works and doesnt throw error?

It only gives warning if the covariance matrix is following:

 [[1,0,1],[0,1,0],[1,1,1]] 

Any explanation for this?

Answer 1

The recipe you're using assumes that your covariance matrix is full rank. Cholesky decomposition is only unique if $\Sigma$ is positive definite. And your comment about inverting a rank-deficient matrix should make it clear why this is a degenerate case. I'm not sure exactly what numpy is doing, but clearly it's smart enough to detect that $\Sigma$ is not full rank and use an alternative strategy. You can read the source to find out how numpy works.

Technically, when $\Sigma$ is singular, you don't define a distribution over random vectors. Perfect correlation implies the $n$-vectors lie in a $n$-plane, which has volume 0, i.e. its $\mathcal{L}^n$ measure is 0, so it's not a distribution on $\mathbb{R}^n$. So your intuition about it "not working" because $\Sigma$ is singular is correct. The numpy implementation is playing fast-and-loose with the definition of "random vector" and inferring that, loosely, you want 3 values that are perfectly correlated. This is well-defined in 1 dimension, so it's just replicating the single value 3 times.

Putting this all together, if the vectors are so correlated that $\Sigma$ is singular, then we don't have a random vector in the technical sense, but it's possible to write software to "side-step" the technical issue so that it "works" in a certain sense.

Answer 2

+1 for the accepted answer. For reference, here is some code that implements drawing samples from a multivariate normal distribution with a possibly rank-deficient covariance matrix:

import numpy as np

def eigh_sample(mean, cov, small=1e-9): 
    s, v = np.linalg.eigh(cov)
    s[abs(s) < small] = 0
    return mean + v * np.sqrt(s) @ np.random.standard_normal(len(mean))

>>> eigh_sample([0,0,0], [[1,1,1],[1,1,1],[1,1,1]])
array([-0.57289804, -0.57289804, -0.57289804])

This is a simple modification of this snippet.