I came across the method for generating multivariate normal samples on wikipedia: https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Drawing_values_from_the_distribution
A widely used method for drawing (sampling) a random vector x from the N-dimensional multivariate normal distribution with mean vector μ and covariance matrix Σ works as follows:
Find any real matrix A such that A AT = Σ. When Σ is positive-definite, the Cholesky decomposition is typically used, and the extended form of this decomposition can always be used (as the covariance matrix may be only positive semi-definite) in both cases a suitable matrix A is obtained. An alternative is to use the matrix A = UΛ½ obtained from a spectral decomposition Σ = UΛUT of Σ. The former approach is more computationally straightforward but the matrices A change for different orderings of the elements of the random vector, while the latter approach gives matrices that are related by simple re-orderings. In theory both approaches give equally good ways of determining a suitable matrix A, but there are differences in computation time.
Let z = (z1, …, zN)T be a vector whose components are N independent standard normal variates (which can be generated, for example, by using the Box–Muller transform).
Let x be μ + Az. This has the desired distribution due to the affine transformation property.
Why does the cholesky decomposition matrix 'A' multiplied by the vector of samples chosen from the standard normal distribution 'z' plus 'μ' give us our result (ie x = μ + Az)?
Why does this work? What is the proof?