Generating Wishart Samples With Bartlett Factors

Why Generate a Wishart Matrix Through Factors?

A Wishart matrix represents a random scatter matrix, so a direct simulation based on its entries must preserve symmetry and positive definiteness. The Bartlett construction achieves both by generating a lower-triangular random factor and multiplying it by its transpose. It also makes the degrees of freedom and scale parameter explicit.

The construction is useful when a simulation needs many Wishart draws, or when a model specifies a nonidentity scale matrix. The essential choices are the dimension \(p\), the degrees of freedom \(\nu\), and a Cholesky factor of the scale matrix. The resulting matrix has dimension \(p \times p\), regardless of the number of random factors used.

Definition: For an integer \(\nu \ge p\) and a positive-definite \(p \times p\) matrix \(\Sigma\), write \(W \sim W_p(\nu,\Sigma)\) for a Wishart matrix with \(\nu\) degrees of freedom and scale matrix \(\Sigma\). Under the Gaussian scatter convention, if \(x_1,\ldots,x_\nu\) are independent \(N_p(0,\Sigma)\) vectors, then \(W=\sum_{j=1}^{\nu}x_jx_j^T\) has this distribution. The scale matrix is \(\Sigma\), not its inverse.

The Bartlett Factor and Its Sampling Rule

The Bartlett factor is a lower-triangular \(p \times p\) random matrix \(T\). Its diagonal entries are positive square roots of chi-square variables; its entries below the diagonal are standard normal variables. All these variables are mutually independent.

Theorem: Let \(\nu \ge p\) be an integer. Construct a lower-triangular matrix \(T\) as follows: \(T_{ii}^2 \sim \chi^2_{\nu-i+1}\) for \(i=1,\ldots,p\); \(T_{ij} \sim N(0,1)\) for \(i>j\); and \(T_{ij}=0\) for \(i<j\). Take all specified random variables independently, and take \(T_{ii}\) to be the positive square root of \(T_{ii}^2\). Then \(TT^T \sim W_p(\nu,I_p)\). If \(\Sigma=LL^T\) is a Cholesky factorization, with \(L\) lower triangular and \(\Sigma\) positive definite, then \(W=LTT^TL^T\) has distribution \(W_p(\nu,\Sigma)\).

There are \(p\) diagonal chi-square draws, with degrees of freedom \(\nu,\nu-1,\ldots,\nu-p+1\), and \(p(p-1)/2\) standard normal draws. The last diagonal chi-square has one degree of freedom when \(\nu=p\), so it is still positive almost surely. Consequently, \(T\) is nonsingular almost surely, and \(TT^T\), as well as \(LTT^TL^T\), is positive definite almost surely.

Factor the scale.
Find the lower-triangular Cholesky factor \(L\) satisfying \(\Sigma=LL^T\).

Draw the Bartlett entries.
Draw the diagonal squares independently from \(\chi^2_{\nu-i+1}\), and draw each entry below the diagonal independently from \(N(0,1)\).

Form the output.
Set \(A=LT\), then return \(W=AA^T=LTT^TL^T\). The product has \(p\) rows and \(p\) columns.

Why the Factorization Gives a Wishart Matrix

The identity-scale result follows by applying Gram–Schmidt to independent Gaussian row vectors. The successive residual lengths produce the diagonal chi-square factors, while the projection coordinates produce the independent standard normal entries.

Proof of the identity-scale result: Let \(Y\) be a \(p \times \nu\) matrix of independent \(N(0,1)\) entries, and denote its rows by \(y_1,\ldots,y_p\), viewed as vectors in \(\mathbb R^\nu\). By the Gaussian scatter definition, \(YY^T \sim W_p(\nu,I_p)\). With probability one, these rows are linearly independent because \(\nu\ge p\) and a Gaussian matrix has full row rank almost surely.

Apply Gram–Schmidt to the rows in order. At step \(i\), let \(q_1,\ldots,q_{i-1}\) be the orthonormal directions obtained from the preceding rows. Decompose \(y_i=\sum_{j=1}^{i-1} t_{ij}q_j+r_i\), where \(r_i\) is perpendicular to those directions. Conditional on the preceding rows, \(y_i\) is a standard Gaussian vector independent of them. In any fixed orthonormal basis, its coordinates are independent \(N(0,1)\). Thus, conditional on the preceding rows, the projection coefficients \(t_{i1},\ldots,t_{i,i-1}\) are independent \(N(0,1)\); the residual has \( \nu-i+1\) independent standard normal coordinates in the perpendicular subspace. Therefore \(\lVert r_i\rVert^2\sim\chi^2_{\nu-i+1}\), independently of those projection coefficients. These conditional distributions do not depend on the preceding rows, so the new coefficients and residual length are also independent of the preceding steps. Induction gives mutual independence of all the coefficients and residual lengths used in the factor.

Set \(t_{ii}=\lVert r_i\rVert\). The entries \(t_{ij}\) for \(j>i\) are zero, producing a lower-triangular \(T\) with exactly the stated Bartlett distributions. The row inner products satisfy \(y_i y_k^T=\sum_{j=1}^{\min(i,k)}t_{ij}t_{kj}\), because the Gram–Schmidt directions are orthonormal. This is the \((i,k)\) entry of \(TT^T\). Hence \(YY^T=TT^T\), proving \(TT^T\sim W_p(\nu,I_p)\). \(\square\)

To obtain the scale-matrix result, use the same \(Y\) and set \(X=LY\). Each column of \(X\) is a linear transformation by \(L\) of an independent \(N_p(0,I_p)\) vector, so the columns are independent \(N_p(0,LL^T)=N_p(0,\Sigma)\) vectors. Their scatter is \(XX^T=LYY^TL^T\). Since \(YY^T\) has the same distribution as \(TT^T\), it follows that \(XX^T\) has the same distribution as \(LTT^TL^T\). By the Gaussian scatter definition, \(XX^T\sim W_p(\nu,\Sigma)\). This proves the stated scale transformation.

Worked Constructions

Worked Example: A Two-Dimensional Identity-Scale Draw

Take \(p=2\) and \(\nu=3\). The diagonal degrees of freedom are \(3\) and \(2\), not \(3\) and \(3\). Consider Bartlett factor values with \(T_{11}^2=3\), \(T_{22}^2=2\), and \(T_{21}=1\). These illustrate the calculation for a draw; they are not claims that a continuous random variable assigns positive probability to any exact value.

The factor and its product are

T= \begin{pmatrix} \sqrt{3} & 0\\ 1 & \sqrt{2} \end{pmatrix}, \qquad TT^T= \begin{pmatrix} 3 & \sqrt{3}\\ \sqrt{3} & 3 \end{pmatrix}.

The factor is \(2 \times 2\), and the output is also \(2 \times 2\). The off-diagonal entry is \(\sqrt{3}\), because it is the inner product of the two rows of \(T\). The diagonal entry in position \((2,2)\) is \(1^2+(\sqrt{2})^2=3\). Since \(\det(T)=\sqrt{6}\ne0\), this particular output is positive definite.

Worked Example: Incorporating a Nonidentity Scale Matrix

Let \(p=2\), \(\nu=4\), and \(\Sigma=\begin{pmatrix}4&2\\2&2\end{pmatrix}\). A Cholesky factor is \(L=\begin{pmatrix}2&0\\1&1\end{pmatrix}\), since \(LL^T=\begin{pmatrix}4&2\\2&2\end{pmatrix}\). The Bartlett diagonal degrees of freedom are \(4\) and \(3\). Use illustrative factor values \(T_{11}^2=5\), \(T_{22}^2=3\), and \(T_{21}=-1\).

First form the Bartlett matrix, then multiply on the left by the scale factor:

T= \begin{pmatrix} \sqrt{5}&0\\ -1&\sqrt{3} \end{pmatrix}, \qquad A=LT= \begin{pmatrix} 2\sqrt{5}&0\\ \sqrt{5}-1&\sqrt{3} \end{pmatrix}.

The output \(W=AA^T\) is

W= \begin{pmatrix} 20&10-2\sqrt{5}\\ 10-2\sqrt{5}&9-2\sqrt{5} \end{pmatrix}.

For example, \(W_{12}=(2\sqrt{5})(\sqrt{5}-1)=10-2\sqrt{5}\), and \(W_{22}=(\sqrt{5}-1)^2+(\sqrt{3})^2=9-2\sqrt{5}\). The multiplication is \(LTT^TL^T=(LT)(LT)^T\); it preserves symmetry and, because \(L\) and \(T\) are nonsingular, gives a positive-definite output.

Worked Example: Checking a Three-Dimensional Bartlett Draw

Take \(p=3\) and \(\nu=4\). The diagonal degrees of freedom must be \(4,3,2\), in that order. Use \(T_{11}=2\), \(T_{22}=\sqrt{3}\), \(T_{33}=\sqrt{2}\), and below-diagonal values \(T_{21}=1\), \(T_{31}=-1\), \(T_{32}=2\). Then

T= \begin{pmatrix} 2&0&0\\ 1&\sqrt{3}&0\\ -1&2&\sqrt{2} \end{pmatrix}, \qquad TT^T= \begin{pmatrix} 4&2&-2\\ 2&4&-1+2\sqrt{3}\\ -2&-1+2\sqrt{3}&7 \end{pmatrix}.

The \((2,3)\) entry is the row inner product \(1(-1)+(\sqrt{3})(2)+0(\sqrt{2})=-1+2\sqrt{3}\). The \((3,3)\) entry is \((-1)^2+2^2+(\sqrt{2})^2=7\). Both the factor and output have three rows and columns. The determinant of \(T\) is \(2\sqrt{6}\ne0\), so the displayed \(TT^T\) is positive definite.

Common Mistakes and Practical Checks

A factor construction is easy to implement, but a small indexing or multiplication error changes the distribution. In particular, the diagonal degrees of freedom decrease as the row index increases, and the Cholesky factor must be applied on both sides of the Bartlett product.

Using the same chi-square degrees of freedom on every diagonal. The \(i\)th diagonal square has \(\nu-i+1\) degrees of freedom. For example, when \(p=3,\nu=4\), the sequence is \(4,3,2\), not \(4,4,4\).
Multiplying by the scale factor on only one side. \(LT T^T\) is not the required symmetric output in general. Form \(A=LT\), then compute \(AA^T=LTT^TL^T\).
Confusing scale with inverse scale. If \(\Sigma=LL^T\), the construction uses \(L\) and yields scale \(\Sigma\). Substituting a factor of \(\Sigma^{-1}\) would generate a different distribution.
Ignoring the dimension restriction for a full-rank factor. The stated construction assumes integer \(\nu\ge p\). If \(\nu<p\) in the Gaussian scatter construction, the scatter has rank at most \(\nu\) and is singular; it cannot be represented by a nonsingular \(p \times p\) Bartlett factor of the form given here.

For a numerical implementation, check that \(L\) is lower triangular and \(LL^T=\Sigma\), that \(T\) is lower triangular with the specified diagonal degrees of freedom, and that the output is computed as \(AA^T\). These checks catch many dimension and parameterization errors before a large simulation is run.

Check Your Understanding

Use the construction and proofs above to answer the following questions.

For \(p=4\) and \(\nu=6\), what are the four chi-square degrees of freedom on the diagonal of \(T\), and how many standard normal entries are below the diagonal?
Why does \(\nu\ge p\) ensure that the final diagonal entry of \(T\) is positive almost surely?
If \(\Sigma=LL^T\), explain why \(LTT^TL^T\) has scale \(\Sigma\) rather than \(\Sigma^{-1}\).
For a \(3 \times 3\) Bartlett factor, which entries contribute to \((TT^T)_{23}\)?
What changes about the rank of the Gaussian scatter if the integer degrees of freedom satisfy \(\nu<p\)?