Can any covariance factorization $LL^\top$ be used for sampling? I thought that any factorization of the for $LL^\top$ of a covariance matrix could be used for correlating random noise according to the covariance. I tried doing this with the following code and three different $LL^\top$ factorizations which each gave separate answers.
Assuming that the covariance is $XX^\top$, I have tried the following, with $z \sim N(0, 1)$:
$$
\begin{aligned}
&X^\top z \\
&V^{\frac{1}{2}}L^\top z \\
&\text{chol}(XX^\top)z
\end{aligned}
$$
where $V$ and $L$ are the eigenvalues and eigenvectors respectively. Am I wrong in thinking that any $LL^\top$ can be used to sample according to the covariance, or have I gone wrong somewhere?
blue = N(0, 1) samples
orange = X^T @ z
green sqrt(V)Lz
red = chol(XX^T)z


from typing import Any, Dict

import pandas as pd  # type: ignore
import seaborn as sns  # type: ignore
import torch
from matplotlib import pyplot as plt  # type: ignore

T = torch.Tensor


def main() -> None:
    for i in range(10):
        cov = torch.randn(2, 2)

        L_full = torch.clone(cov)
        V_eig, L_eig = torch.linalg.eig(cov.mm(cov.T))
        V_eig, L_eig = V_eig.float().sqrt(), L_eig.float()
        L_chol = torch.linalg.cholesky(cov.mm(cov.T))

        samples = torch.randn(2, 1000)

        def make(t: str, samples: T) -> Dict[str, Any]:
            return {"type": [t] * 1000, "x": samples[0, :].tolist(), "y": samples[1, :].tolist()}

        factor_sample = L_full.T @ samples
        eig_sample = torch.diag_embed(V_eig) @ L_eig  @ samples
        chol_sample = L_chol @ samples

        out: Dict[str, Any] = {"type": [], "x": [], "y": []}
        for d in map(make, ("samples", "factor", "eig", "chol"), (samples, factor_sample, eig_sample, chol_sample)):  # type: ignore
            for k in d:
                out[k] += d[k]

        sns.kdeplot(data=pd.DataFrame(out), x="x", y="y", hue="type", ax=plt.gca())
        plt.show()


if __name__ == "__main__":
    main()


EDIT
I can see now that the covariance factors are invariant to rotations since $L Q^\top Q L^\top = LL^\top$ so it must be that the other non-cholesky factors have some random rotation $Q$ added to them.
My question now becomes:

*

*Why do the non-cholesky factors have a rotation added to them?

*Why is the cholesky factor the "correct" rotation. (I am assuming it is correct because it is what is used in libraries.)


Looks like some simple mixups with matrix transposes. The following works for me:
import torch
cov = torch.randn(2, 2)

L_full = torch.clone(cov)
V_eig, L_eig = torch.eig(cov.mm(cov.T),eigenvectors=True)
V_eig, L_eig = V_eig[:,0].float().sqrt(), L_eig.float().T
L_chol = torch.linalg.cholesky(cov.mm(cov.T))

samples = torch.randn(2, 5000)

def make(t,samples):
    return {"type": [t] * 5000, "x": samples[0, :].tolist(), "y": samples[1, :].tolist()}

A=L_full.T
B=torch.diag_embed(V_eig) @ L_eig  
C=L_chol.T
factor_sample = A.T @ samples
eig_sample = B .T  @ samples
chol_sample = C.T @ samples

print(np.cov(factor_sample))
print(np.cov(eig_sample))
print(np.cov(chol_sample))

In general,if $A^TA=B^TB$, then $A=QB$ for some orthogonal matrix $Q$. To see this, note that $A$ and $B$ must have the same singular values, so we have $A=U_1\Sigma V_1$ and $B=U_2\Sigma V_2$ where $\Sigma$ is the diagonal matrix of singular values. By the condition, $V_1 \Sigma^2 V_1^T=V_2\Sigma^2 V_2^T$, and by uniqueness of the SVD, we see that $V_1=V_2$ (if there are repeated singular values, then this becomes a little more delicate). So we can take $Q=U_1U_2^t$. This explains why the different samples are all related by rotations. As far as which is the ``correct" rotation, it is truly immaterial. This is because if $z\sim N(0,I_n)$ and $A^TA=B^TA$, then $Az$ and $Bz$ have the same distribution, namely $N(0,A^TA)=N(0,B^TB)$.
