I am trying to find the time complexity of sampling from a multivariate normal distribution with an $n\times n$ covariance matrix $C$. I have a found an answer here and have a couple of questions.
The answer says if we decompose $C=LL^T$ once (whose complexity is $O(n^3)$), then $Lx$ has the desired covariance structure. Here $x$ is an $n-$vector of independent standard normal random variables. I don't understand why $Lx$ has the desired covariance structure?
Is it true that sampling from a standard normal distribution (or any continuous distribution) is $O(1)$ and hence the complexity of sampling $x$ is $O(n)$? Sampling complexity clearly depends on the algorithm used but I was wondering if say the state of the art algorithm in Python has the complexity of $O(1)$ for sampling from a continuous distribution?
If $X \sim \mathcal{N}(0, I)$ then $LX \sim \mathcal{N}(0, LL^T) = \mathcal{N}(0, C)$ because $L$ is a linear transformation. One way to prove this is with characteristic functions, see the wikipedia page.
One way to draw from standard normal distribution is to compute $\Phi^{-1}(U)$, where $\Phi$ is the CDF of the normal distribution and $U$ is a uniform distribution on $[0,1]$. You need to start with a uniform random variable because that is what the pseudorandom number generator gives you. Another option is the Box-Muller transform. I'm not sure what exactly numpy is doing, but you can look at the source code. Both of these approaches would be $O(1)$ assuming that a uniform random variable can be generated in constant time.