[proofplan] Represent the Wishart matrix as $W=X^\top X$, where the columns of $X$ are independent standard Gaussian vectors in $\mathbb{R}^n$. Apply Gram-Schmidt to these columns and write $X=QR$, with $Q$ having orthonormal columns and $R$ upper triangular with positive diagonal; then $W=R^\top R$, so $A=R^\top$ is the desired Cholesky factor. The distributional content follows from rotational invariance of the standard Gaussian law: at the $j$th step, the Gaussian coordinates in the already constructed orthonormal directions are independent $\mathcal{N}(0,1)$ variables, while the squared norm of the orthogonal residual is chi-squared with $n-j+1$ degrees of freedom and independent of those coordinates. [/proofplan]

[step:Construct the triangular factor by Gram-Schmidt] Let \begin{align*} X_1,\dots,X_p:\Omega \to \mathbb{R}^n \end{align*} denote the column random vectors of $X$. Let $\mathcal{N}_n(0,I_n)$ denote the standard Gaussian probability law on $\mathbb{R}^n$, meaning the law of a random vector whose $n$ coordinates in the standard basis are independent $\mathcal{N}(0,1)$ random variables. Since the entries of $X$ are independent $\mathcal{N}(0,1)$ random variables, the vectors $X_1,\dots,X_p$ are independent and each has distribution $\mathcal{N}_n(0,I_n)$. We recursively define orthonormal vectors $q_1,\dots,q_p$ and upper triangular entries $r_{ij}$ for $1 \leq i \leq j \leq p$. Set \begin{align*} r_{11} := |X_1|, \qquad q_1 := \frac{X_1}{r_{11}}. \end{align*} For $j \geq 2$, assuming $q_1,\dots,q_{j-1}$ have been defined, define \begin{align*} r_{ij} &:= q_i \cdot X_j, && 1 \leq i < j, \\ Y_j &:= X_j - \sum_{i=1}^{j-1} r_{ij} q_i, \\ r_{jj} &:= |Y_j|, \\ q_j &:= \frac{Y_j}{r_{jj}}. \end{align*} Because $X_1,\dots,X_j$ have a joint density on $\mathbb{R}^{nj}$ and $j \leq p \leq n$, the event that $X_1,\dots,X_j$ are linearly dependent has probability zero. Hence $r_{jj}>0$ for every $j$ almost surely, so the construction is valid almost surely. Let \begin{align*} Q:\Omega \to \mathbb{R}^{n \times p} \end{align*} be the random matrix with columns $q_1,\dots,q_p$, and let \begin{align*} R:\Omega \to \mathbb{R}^{p \times p} \end{align*} be the upper triangular random matrix with entries $R_{ij}=r_{ij}$ for $i \leq j$ and $R_{ij}=0$ for $i>j$. The recursion gives \begin{align*} X_j = \sum_{i=1}^j r_{ij} q_i \end{align*} for each $1 \leq j \leq p$, so \begin{align*} X = QR. \end{align*} Since the columns of $Q$ are orthonormal, $Q^\top Q=I_p$. Therefore \begin{align*} W = X^\top X = R^\top Q^\top Q R = R^\top R. \end{align*} Let $E \subset \Omega$ be the full-probability event on which the Gram-Schmidt construction above is valid for every $1 \leq j \leq p$. On $\Omega \setminus E$, define each $q_j$ to be the deterministic standard basis vector $u_j \in \mathbb{R}^n$ and define each $r_{ij}$ with $1 \leq i \leq j \leq p$ to be $0$; this gives everywhere-defined random variables without changing any almost-sure identity or distributional conclusion. Define \begin{align*} A:\Omega \to \mathbb{R}^{p \times p} \end{align*} by \begin{align*} A(\omega) := \begin{cases} R(\omega)^\top, & \omega \in E, \\ I_p, & \omega \in \Omega \setminus E. \end{cases} \end{align*} Then $A$ is a lower triangular random matrix on all of $\Omega$, its diagonal entries are positive on $E$, and \begin{align*} W=AA^\top \end{align*} on $E$, hence almost surely. [/step]

[step:Identify the distribution of one Gram-Schmidt step]Fix $j \in \{1,\dots,p\}$. Let $\mathcal{F}_{j-1}$ denote the $\sigma$-algebra generated by $X_1,\dots,X_{j-1}$; for $j=1$, let $\mathcal{F}_0 := \{\varnothing,\Omega\}$. On the full-probability event where $q_1,\dots,q_{j-1}$ are defined, extend this orthonormal family to an orthonormal basis \begin{align*} q_1,\dots,q_{j-1},e_j,\dots,e_n \end{align*} of $\mathbb{R}^n$, where each $e_k$ is chosen as an $\mathcal{F}_{j-1}$-measurable vector. To construct the completion measurably, process the fixed deterministic basis $u_1,\dots,u_n$ of $\mathbb{R}^n$ in order, project each $u_m$ onto the orthogonal complement of the span of the vectors already selected, and keep the first vector whose projected norm is nonzero. The ordering gives a fixed tie-breaking rule on the full-probability rank event; on its complement, define the remaining $e_k$ arbitrarily, for instance by the same deterministic basis. Each projection and norm is a Borel function of $q_1,\dots,q_{j-1}$, so the selected vectors are $\mathcal{F}_{j-1}$-measurable. Define the $\mathcal{F}_{j-1}$-measurable orthogonal matrix \begin{align*} O_j:\Omega \to \mathbb{R}^{n \times n} \end{align*} whose rows are $q_1,\dots,q_{j-1},e_j,\dots,e_n$ in that order. Since $X_j$ is independent of $\mathcal{F}_{j-1}$ and has distribution $\mathcal{N}_n(0,I_n)$, conditioning on $\mathcal{F}_{j-1}$ fixes the matrix $O_j$. For every deterministic orthogonal matrix $O \in \mathbb{R}^{n \times n}$, the transformed vector $OX_j$ again has law $\mathcal{N}_n(0,I_n)$ by rotational invariance of the standard Gaussian law. Therefore, conditional on $\mathcal{F}_{j-1}$, the coordinates of $X_j$ in this orthonormal basis are independent standard normal random variables. Define these coordinate random variables by \begin{align*} Z_i &:= q_i \cdot X_j, && 1 \leq i < j, \\ Z_k &:= e_k \cdot X_j, && j \leq k \leq n. \end{align*} Then, conditional on $\mathcal{F}_{j-1}$, the family \begin{align*} Z_1,\dots,Z_{j-1},Z_j,\dots,Z_n \end{align*} consists of independent $\mathcal{N}(0,1)$ random variables. By construction, for $1 \leq i<j$, \begin{align*} r_{ij}=q_i \cdot X_j=Z_i. \end{align*} The residual vector is the [orthogonal projection](/theorems/437) of $X_j$ onto the orthogonal complement of $\operatorname{span}\{q_1,\dots,q_{j-1}\}$: \begin{align*} Y_j = \sum_{k=j}^n Z_k e_k. \end{align*} Hence \begin{align*} r_{jj}^2 = |Y_j|^2 = \sum_{k=j}^n Z_k^2. \end{align*} Therefore, conditional on $\mathcal{F}_{j-1}$, \begin{align*} r_{1j},\dots,r_{j-1,j} \end{align*} are independent $\mathcal{N}(0,1)$ random variables, $r_{jj}^2$ has distribution $\chi^2_{n-j+1}$, and $r_{jj}$ is independent of $r_{1j},\dots,r_{j-1,j}$.[/step]

[guided]Fix the column index $j$. The information from the previous columns is encoded in the $\sigma$-algebra $\mathcal{F}_{j-1}$ generated by $X_1,\dots,X_{j-1}$. Once this information is fixed, the vectors $q_1,\dots,q_{j-1}$ are deterministic orthonormal vectors. We complete $q_1,\dots,q_{j-1}$ to an orthonormal basis \begin{align*} q_1,\dots,q_{j-1},e_j,\dots,e_n \end{align*} of $\mathbb{R}^n$. The vectors $e_j,\dots,e_n$ are chosen measurably from the previous data, so conditioning on $\mathcal{F}_{j-1}$ makes the entire basis fixed. The crucial input is rotational invariance of the standard Gaussian vector. Define the $\mathcal{F}_{j-1}$-measurable orthogonal matrix \begin{align*} O_j:\Omega \to \mathbb{R}^{n \times n} \end{align*} whose rows are $q_1,\dots,q_{j-1},e_j,\dots,e_n$. Since $X_j$ is independent of $\mathcal{F}_{j-1}$ and $X_j\sim \mathcal{N}_n(0,I_n)$, conditioning on $\mathcal{F}_{j-1}$ fixes $O_j$ while leaving $X_j$ with standard Gaussian law. For each deterministic orthogonal matrix $O \in \mathbb{R}^{n \times n}$, rotational invariance gives $OX_j\sim \mathcal{N}_n(0,I_n)$, so its coordinates are independent $\mathcal{N}(0,1)$ random variables. Thus, conditional on $\mathcal{F}_{j-1}$, the coordinate variables \begin{align*} Z_i &:= q_i \cdot X_j, && 1 \leq i < j, \\ Z_k &:= e_k \cdot X_j, && j \leq k \leq n \end{align*} are independent standard normal random variables. The Gram-Schmidt coefficients before the diagonal are exactly the coordinates of $X_j$ in the already constructed directions: \begin{align*} r_{ij}=q_i\cdot X_j=Z_i,\qquad 1\leq i<j. \end{align*} The residual removes precisely those coordinates, so it lies in the orthogonal complement: \begin{align*} Y_j = X_j-\sum_{i=1}^{j-1}(q_i\cdot X_j)q_i = \sum_{k=j}^n Z_k e_k. \end{align*} Taking the Euclidean norm and using orthonormality gives \begin{align*} r_{jj}^2=|Y_j|^2=\sum_{k=j}^n Z_k^2. \end{align*} A sum of squares of $n-j+1$ independent standard normal variables has distribution $\chi^2_{n-j+1}$. Because the variables $Z_1,\dots,Z_{j-1}$ are independent of $Z_j,\dots,Z_n$, the off-diagonal coefficients $r_{1j},\dots,r_{j-1,j}$ are independent of the diagonal length $r_{jj}$.[/guided]

[step:Promote the conditional distributions to joint independence] For each $j$, the preceding step gives the conditional distribution of the new block \begin{align*} B_j := (r_{1j},\dots,r_{j-1,j},r_{jj}) \end{align*} given $\mathcal{F}_{j-1}$. This conditional distribution does not depend on the realised values of $X_1,\dots,X_{j-1}$. Moreover, within $B_j$, the variables $r_{1j},\dots,r_{j-1,j}$ and $r_{jj}$ are independent, with \begin{align*} r_{ij} &\sim \mathcal{N}(0,1), && 1 \leq i<j, \\ r_{jj}^2 &\sim \chi^2_{n-j+1}. \end{align*} Since the conditional law of $B_j$ given $\mathcal{F}_{j-1}$ is fixed and does not depend on $\mathcal{F}_{j-1}$, the block $B_j$ is independent of $\mathcal{F}_{j-1}$. Because $B_1,\dots,B_{j-1}$ are functions of $X_1,\dots,X_{j-1}$, the block $B_j$ is independent of $B_1,\dots,B_{j-1}$. Induction over $j=1,\dots,p$ shows that all entries \begin{align*} \{r_{ij}:1\leq i\leq j\leq p\} \end{align*} are mutually independent, with the distributions just stated. [/step]

[step:Translate the upper triangular QR factor into the lower triangular Bartlett factor] Recall that $A=R^\top$. Hence, for $1 \leq j \leq p$, \begin{align*} a_{jj}=r_{jj}, \end{align*} and for $1 \leq j<i\leq p$, \begin{align*} a_{ij}=r_{ji}. \end{align*} The independence of the triangular entries of $R$ therefore gives the independence of the lower triangular entries of $A$. The diagonal distributions become \begin{align*} a_{jj}^2=r_{jj}^2\sim \chi^2_{n-j+1}, \qquad 1\leq j\leq p, \end{align*} and the strictly lower triangular distributions become \begin{align*} a_{ij}=r_{ji}\sim \mathcal{N}(0,1), \qquad 1\leq j<i\leq p. \end{align*} Together with $W=AA^\top$ almost surely, this proves the Bartlett decomposition. [/step]