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.)

 A: 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)$.
