I am trying to intuitively understand the Cholesky decomposition in gaussian process function sampling. I understand it as as the square root of the covariance matrix being the multivariate generalization of taking the sqrt of the variance and then transforming a standard normal variable. Let $z \sim \mathcal{N}(0, 1)$ and $\mathbf{LL^*} = Cholesky(\mathbf{\Sigma})$.

$$ y = z\sigma + \mu \\ \mathbf{y} = \mathbf{L}\mathbf{z} + \mathbf{\mu} $$

But what I am having a hard time justifying is this... $\mathbf{L}$ is a lower triangular matrix, so when I sample a GP with a vector of standard normal variables (like this https://katbailey.github.io/post/gaussian-processes-for-dummies/), the first output will only be a combination of a single column of $\mathbf{L}$ since the first row of $\mathbf{L}$ has a single entry. For example...

$$ \mathbf{Lz} = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 3 & 0 \\ 4 & 5 & 6 \end{bmatrix} \begin{bmatrix} z \\ z \\ z \end{bmatrix} $$

It seems kind of arbitrary and unnatural that an unbalanced matrix like this would work out to be an even sample when the last entry in $\mathbf{z}$ is a combination of all columns and the first entry in $\mathbf{z}$ is only a combination of one column.

Is there any intuition for this? I hope my explanation is clear.

  • $\begingroup$ The lower triangular is not the only possible square-root of a symmetric positive definite matrix, perhaps that's what making it harder for your understanding. Taking an unintuitive detour, have you considered actually computing the covariance of $\mathbb y$? $\endgroup$
    – Firebug
    Apr 2, 2021 at 13:00
  • $\begingroup$ In the vector case you must use a vector $\mathbf{z}$ with mean zero and its covariance set to the identity matrix. So in your formula for $\mathbf{L}\mathbf{z}$ you need three random variables $z_1$, $z_2$ and $z_3$ (which are independent). $\endgroup$
    – Yves
    Apr 2, 2021 at 13:48

1 Answer 1


$\newcommand{\0}{\mathbf 0}$Let $z \sim \mathcal N_n(\0, I)$ and $LL^T=\Sigma$. I'll ignore $\mu$ as I don't think it adds anything here.

Here's how we can think of $Lz$ as forming our multivariate Gaussian $X \sim \mathcal N_n(\0, \Sigma)$.

My main intuition here is that each new row of $L$ adds in a new independent variable that hasn't appeared in any of the previous terms, so depending on how we weight them we can marginally get a Gaussian with the right variance, since sums of independent Gaussians are Gaussian, and we can construct the correlation that we want based on how much weight we give to the new independent $z_i$ that we're adding in relative to the other $z_j$ already used in previous $X_j$.

More explicitly, let's use $\ell_i$ as the $i$th row of $L$ (I'll still treat them as column vectors). In $\ell_1$ we're just going to have a single non-zero element, namely $L_{11}$, and all this needs to do is get $\ell_1^Tz \stackrel{\text d}= X_1$. Taking $L_{11} = \sqrt{\Sigma_{11}}$ will do the job here as we are just scaling $z_1$ into $X_1$.

Next we'll have $\ell_2 = (L_{21}, L_{22}, 0, \dots, 0)^T$ so we're looking to get $X_2$ out of $L_{21}z_1 + L_{22} z_2$. Having $L_{22}$ non-zero is what allows us to prevent $X_1$ and $X_2$ from being perfectly correlated; analogously, if any diagonal elements of $L$ are zero then $LL^T$ is low rank and we have a singular covariance matrix which confirms the perfect correlation. On the other hand if we had $L_{21} = 0$ then we'd have $X_1$ and $X_2$ uncorrelated (and therefore independent since they're jointly Gaussian). The amount of $z_1$ and $z_2$ that we blend together is what allows us to go between these two extremes, and the overall magnitude of the weights $\|\ell_2\|^2$ is what controls the marginal variance of $X_2$.

We continue this process where for $X_{i+1}$ we take a linear combination of $z_1, \dots, z_i$ and "mix in" enough of $z_{i+1}$ to get the right correlation structure, and the overall magnitude of the weights can be scaled to get $\text{Var}(X_{i+1})$ right.

Related to this, we know that the eigenvalues of a triangular matrix are on the diagonal and $L$ (and therefore $\Sigma$) will be full rank if and only if all of these diagonal elements are positive. These diagonal elements are exactly the weights that we give to the previously-unseen fluctuation $z_{i+1}$ when making $X_{i+1}$, so this is one way to see that the closer $\Sigma$ is to being low rank, the more correlated elements of $X$ are.

Update: I think of this process as iteratively building $\Sigma$. When we start with $X_1 = \ell_1^Tz = L_{11}z_1$ we correctly get the upper $1\times 1$ block of $\Sigma$ that governs the distribution of $X_1$. Next we add in $X_2$ and $$ \text{Cov}(X_1, X_2) = \text{Cov}(X_2, X_1) = \text{Cov}(\ell_1^Tz, \ell_2^Tz) = \ell_1^T \text{Var}(z)\ell_2 = \ell_1^T\ell_2 $$ so this shows that $\Sigma_{12}$ and $\Sigma_{21}$ come from the overlap between $\ell_1$ and $\ell_2$ (which makes sense since $z_2$ is independent of $z_1$ so the only $z_i$ that contribute to the covariance are the ones that appear in both, and that's the elements of $z$ corresponding to $\ell_1$ and $\ell_2$ both being nonzero). We also get $\Sigma_{22}$ from $\|\ell_2\|^2$.

In general we'll have $\text{Cov}(X_i, X_j) = \ell_i^T\ell_j$ and this comes from $z_1$ through $z_{\min\{i,j\}}$ since that's as much overlap as there is between the $z_k$ going into $X_i$ and $X_j$. So as we advance $X_i$ to $X_{i+1}$ we fill in the covariances between $X_{i+1}$ and $X_1, \dots, X_i$ so we go from the upper $i\times i$ block of $\Sigma$ to the upper $(i+1)\times(i+1)$, and so on.

  • $\begingroup$ Thank you for the detailed answer, it was a very nice read and confirmed most of what I already knew, but I think it missed the main point of my question. The main point is that all the $\mathbf{x_i}$'s should be jointly Gaussian, but in this case the first row of $\mathbf{L}$ is only getting mixed in correlations from itself and none of the other off diagonal elements. How can it be that the last $z_i$ gets correlations from every other element in $\mathbf{\Sigma}$, but the first $z_i$ only gets one? $\endgroup$
    – Joff
    Apr 8, 2021 at 2:07
  • $\begingroup$ @deltaskelta I've just updated to hopefully address your comment $\endgroup$
    – jld
    Apr 8, 2021 at 15:20
  • $\begingroup$ maybe I am just really dense, but I still don't get how that answers the question. You say "As we move from $X_i$ to $X_{i+1}$ we fill in the covariances" so we still have the $z_1$ is just a standard normal sample with no covariances built into the equation, right? $\endgroup$
    – Joff
    Apr 9, 2021 at 1:14
  • $\begingroup$ @deltaskelta in your example we have $\ell_1 = (1, 0, 0)^T$ and $\ell_3 = (4, 5, 6)^T$ so $X_1 = \ell_1^Tz = z_1$ and $X_3 = \ell_3^Tz = 4z_1 + 5z_2 + 6z_3$. Then $$\text{Cov}(X_1, X_3) = \text{Cov}(z_1, 4z_1 + 5z_2 + 6z_3) = 4 = \ell_1^T\ell_3$$since $z_1 \perp z_2, z_3$, so the covariance between $X_1$ and $X_3$ is entirely from the common $z_1$ they share, since the terms that $X_3$ has that $X_1$ doesn't have are all independent of the terms in $X_1$. Does that help? $\endgroup$
    – jld
    Apr 9, 2021 at 15:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.