0
$\begingroup$

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?

Improve this question
$\endgroup$
12
  • 1
    $\begingroup$ If $x \sim \mathcal{N}(0,1^2)$ and $y$ is perfectly correlated with $x$, what value do you suppose that $y_1$ takes when you know $x_1=-0.896$? How can you use this information to construct a sampler for $y$? $\endgroup$ – Sycorax♦ Feb 6 '20 at 16:15
  • $\begingroup$ both the values will be same. the sampler of y shall be same as x. But thats why I am confused, according to formulation the sampler has to perform cholesky decomposition of the matrix with rank 1, which shall throw an error. Reference $\endgroup$ – GENIVI-LEARNER Feb 6 '20 at 16:21
  • 1
    $\begingroup$ 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. $\endgroup$ – Sycorax♦ Feb 6 '20 at 16:26
  • 1
    $\begingroup$ Correcting my previous comment. Your comment is correct -- we can have correlated vectors up to a certain degree, but if the vectors are so correlated that $\Sigma$ is singular, then we don't have a random vector in the technical sense. $\endgroup$ – Sycorax♦ Feb 6 '20 at 17:23
  • 3
    $\begingroup$ Does this answer your question? Generate normally distributed random numbers with non positive-definite covariance matrix (I probably should have searched for duplicates before writing my answer but here we are...) $\endgroup$ – Sycorax♦ Feb 6 '20 at 17:45
2
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ The software would have to perform lots of checks in such as where a scenario like in Gaussian Process where say taking 3000 samples and where say $x_{200}$ is perfectly correlated with $x_{1500}$ and various permutations of the concept! Since you clarified some doubts, I shall be looking into how is this dealt with in GP $\endgroup$ – GENIVI-LEARNER Feb 6 '20 at 17:55
  • 1
    $\begingroup$ I'm not sure that how hard it is to check this. It probably checks if the matrix is full rank, and if it's not, falls back to using an alternative strategy. Factorizing a very large matrix is expensive, but... that's the nature of working with large matrices. $\endgroup$ – Sycorax♦ Feb 6 '20 at 17:57
  • $\begingroup$ It is not that hard, see my answer. $\endgroup$ – marnix Apr 29 '20 at 21:44
  • $\begingroup$ @marnix cool! Looks like you implemented what my comment describes. $\endgroup$ – Sycorax♦ Apr 29 '20 at 21:50
0
$\begingroup$

+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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.