Cholesky decomposition lower triangular in Gaussian process sampling 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.
 A: $\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.
