Sampling Matrix-Beta Variables From Wishart Factors

From Wishart Factors to Matrix-Beta Samples

Matrix-beta distributions describe random positive-definite matrices, but they are often easiest to sample indirectly. The useful construction starts with two independent Wishart matrices and applies a matrix transformation that preserves symmetry and enforces the appropriate support. The transformation must be a congruence, not entrywise division or multiplication.

Throughout, $A\sim W_p(a,I_p)$ means that $A$ is a $p\times p$ Wishart matrix with identity scale and degrees of freedom $a$. It has density proportional to $|A|^{(a-p-1)/2}\exp\{-\operatorname{tr}(A)/2\}$ on the positive-definite matrices, and zero elsewhere. We require $a>p-1$ for a nonsingular Wishart distribution. The same conventions apply to $B\sim W_p(b,I_p)$.

Definition: If $A\sim W_p(a,I_p)$ and $B\sim W_p(b,I_p)$ are independent, the matrix beta type I construction is $$ X=(A+B)^{-1/2}A(A+B)^{-1/2}, $$ and the matrix beta type II construction is $$ F=B^{-1/2}AB^{-1/2}. $$ Here $M^{1/2}$ denotes the unique symmetric positive-definite square root of a positive-definite matrix $M$. The respective parameter pairs are $(a/2,b/2)$.

For computation, a Cholesky factor is often more convenient than a symmetric square root. If $S=LL^T$, where $L$ is lower triangular with positive diagonal, then $L^{-1}AL^{-T}$ is the Cholesky-whitened version of $A$. It is symmetric even though $L^{-1}$ need not be. Both the symmetric-root and Cholesky versions produce the corresponding matrix-beta law.

Generating the Wishart Factors

There are two standard ways to generate a nonsingular identity-scale Wishart factor.

Gaussian construction: for integer $m\ge p$, generate an $m\times p$ matrix $G$ with independent $N(0,1)$ entries and set $W=G^TG$. Then $W\sim W_p(m,I_p)$, and $W$ is positive definite with probability one.
Bartlett construction: for any real $m>p-1$, generate a lower-triangular matrix $T$. Its diagonal entries are independent, with $T_{ii}^2\sim\chi^2_{m-i+1}$; its entries below the diagonal are independent $N(0,1)$, and all these entries are mutually independent. Set $W=TT^T$. Then $W\sim W_p(m,I_p)$.

In the Bartlett construction, the final diagonal degree of freedom is $m-p+1>0$, which is why $m>p-1$ is the nonsingularity condition. To obtain independent factors $A$ and $B$, use independent random draws for their two Gaussian matrices or two Bartlett factors.

Worked Example: A Bartlett Draw for Beta Type I

Take $p=2$, $a=5$, and $b=4$. For this illustrative realization, suppose the two independent Bartlett factors are $$ T_A=\begin{pmatrix}2&0\\0&\sqrt{3}\end{pmatrix}, \qquad T_B=\begin{pmatrix}\sqrt{2}&0\\0&\sqrt{3}\end{pmatrix}. $$ These are valid Bartlett outcomes: the squared diagonal entries for $T_A$ have the required chi-square degrees of freedom $5,4$, and those for $T_B$ have degrees of freedom $4,3$. The below-diagonal entries are zero in this particular illustrative draw.

Form the factors: $$ A=T_AT_A^T=\begin{pmatrix}4&0\\0&3\end{pmatrix}, \qquad B=T_BT_B^T=\begin{pmatrix}2&0\\0&3\end{pmatrix}. $$ Thus $S=A+B=6I_2$, so $S^{-1/2}=I_2/\sqrt{6}$. The beta type I sample is $$ X=S^{-1/2}AS^{-1/2} =\begin{pmatrix}2/3&0\\0&1/2\end{pmatrix}. $$ Its diagonal entries are strictly between zero and one. In particular, both $X$ and $I_2-X=\operatorname{diag}(1/3,1/2)$ are positive definite, as required.

Why the Type I Transformation Has the Beta Law

The transformation works because it separates the total matrix $S=A+B$ from the relative contribution of $A$. The following derivation also explains why the Cholesky implementation is valid.

Theorem: Suppose $A\sim W_p(a,I_p)$ and $B\sim W_p(b,I_p)$ are independent, with $a,b>p-1$. Set $S=A+B$, let $S=LL^T$ be its Cholesky factorization, and define $X=L^{-1}AL^{-T}$. Then $S$ and $X$ are independent, and $X$ has density proportional to $$ |X|^{(a-p-1)/2}|I_p-X|^{(b-p-1)/2} $$ on the set where $X$ and $I_p-X$ are positive definite.

Proof. The joint density of $A,B$ is proportional to $$ |A|^{(a-p-1)/2}|B|^{(b-p-1)/2} \exp\{-\operatorname{tr}(A+B)/2\}. $$ Use the change of variables $S=A+B$, $A=LXL^T$, with $L$ the Cholesky factor of $S$. Then $B=L(I_p-X)L^T$. For fixed $S$, the Jacobian of the congruence map $X\mapsto LXL^T$ on symmetric matrices is $|L|^{p+1}=|S|^{(p+1)/2}$. The full change of variables has this same Jacobian factor because the $S$ coordinate is retained, making its derivative block triangular.

The determinants and trace become $$ |A|=|S||X|,\qquad |B|=|S||I_p-X|,\qquad \operatorname{tr}(A+B)=\operatorname{tr}(S). $$ After including the Jacobian, the joint density factors into a term in $S$, $$ |S|^{(a+b-p-1)/2}\exp\{-\operatorname{tr}(S)/2\}, $$ and a term in $X$, $$ |X|^{(a-p-1)/2}|I_p-X|^{(b-p-1)/2}. $$ The first term is integrable because $a+b>p-1$; the second is supported on $X>0$ and $I_p-X>0$. The factorization proves independence and the stated density. The principal symmetric-root version has the same law: for each fixed $S$, its whitened matrix and the Cholesky-whitened matrix differ by orthogonal conjugation, and the displayed density is invariant under such conjugation. $\square$

Worked Example: A Non-diagonal Type I Sample and Its Support

Consider $p=2$, $a=b=3$, and the following possible positive-definite factor matrices: $$ A=\begin{pmatrix}3&1\\1&2\end{pmatrix}, \qquad B=\begin{pmatrix}3&-1\\-1&4\end{pmatrix}. $$ Both can occur as Wishart-factor outputs. For example, $A=G_A^TG_A$ and $B=G_B^TG_B$, where $$ G_A=\begin{pmatrix}\sqrt{3}&1/\sqrt{3}\\0&\sqrt{5/3}\end{pmatrix}, \qquad G_B=\begin{pmatrix}\sqrt{3}&-1/\sqrt{3}\\0&\sqrt{11/3}\end{pmatrix}. $$ Indeed, multiplying each matrix by its transpose in the order $G^TG$ gives the stated $A$ and $B$. Their determinants are $5$ and $11$, respectively, and their first diagonal entries are positive, so each is positive definite.

Their sum is $S=6I_2$, hence $$ X=S^{-1/2}AS^{-1/2} =\frac{1}{6}\begin{pmatrix}3&1\\1&2\end{pmatrix} =\begin{pmatrix}1/2&1/6\\1/6&1/3\end{pmatrix}. $$ For a symmetric $2\times2$ matrix, positive definiteness follows if its upper-left entry and determinant are positive. Here $X_{11}=1/2>0$ and $$ |X|=\frac{1}{2}\frac{1}{3}-\left(\frac{1}{6}\right)^2 =\frac{5}{36}>0. $$ Also $$ I_2-X=\begin{pmatrix}1/2&-1/6\\-1/6&2/3\end{pmatrix}, \qquad |I_2-X|=\frac{1}{2}\frac{2}{3}-\left(-\frac{1}{6}\right)^2 =\frac{11}{36}>0. $$ Thus $0<X<I_2$, establishing the strict beta type I support check directly.

Type II Sampling and Its Density

For beta type II, whiten $A$ by the denominator factor $B$, rather than by the sum. Given independent $A\sim W_p(a,I_p)$ and $B\sim W_p(b,I_p)$, compute $F=B^{-1/2}AB^{-1/2}$. Alternatively, if $B=LL^T$, use $F=L^{-1}AL^{-T}$. The support is the full cone of positive-definite matrices.

Result: Under the Wishart assumptions above, the type II sample $F=B^{-1/2}AB^{-1/2}$ has density proportional to $$ |F|^{(a-p-1)/2}|I_p+F|^{-(a+b)/2},\qquad F>0. $$

To see the form of this density, hold $B$ fixed and write $A=B^{1/2}FB^{1/2}$. The Jacobian from $F$ to $A$ is $|B|^{(p+1)/2}$. In the product Wishart density, this gives an $F$-factor $|F|^{(a-p-1)/2}$, while the remaining $B$-dependent determinant power is $|B|^{(a+b-p-1)/2}$. The exponential term is $\exp\{-\operatorname{tr}(B(I_p+F))/2\}$. Integrating over $B>0$, change variables to $H=(I_p+F)^{1/2}B(I_p+F)^{1/2}$. The determinant and Jacobian contribute the factor $|I_p+F|^{-(a+b)/2}$; the remaining integral over $H>0$ is independent of $F$. This establishes the stated density form without needing its normalizing constant.

Worked Example: A Cholesky Type II Sample

Let $p=2$, $a=3$, $b=4$, and suppose a pair of Wishart draws is $$ A=\begin{pmatrix}2&1\\1&2\end{pmatrix}, \qquad B=I_2. $$ Both factors are positive definite: $|A|=3>0$, and $B$ is the identity. The Cholesky factor of $B$ is $L=I_2$, so $$ F=L^{-1}AL^{-T}=A=\begin{pmatrix}2&1\\1&2\end{pmatrix}. $$ The leading diagonal entry is $2>0$, and $|F|=3>0$, so $F$ is positive definite. Its eigenvalues are $1$ and $3$, also confirming that it lies in the type II support. The choice $B=I_2$ makes the arithmetic transparent; in actual simulation, $B$ is generated as an independent Wishart draw and need not equal the identity.

A Reliable Sampling and Checking Procedure

Check dimensions and degrees of freedom.
Choose $p$ and require $a,b>p-1$ for nonsingular factors. For Gaussian construction, use integer degrees of freedom at least $p$.

Generate independent factors.
Use separate Gaussian matrices or Bartlett decompositions to produce $A$ and $B$.

Apply the matching transform.
For type I, factor $S=A+B=LL^T$ and set $X=L^{-1}AL^{-T}$. For type II, factor $B=LL^T$ and set $F=L^{-1}AL^{-T}$.

Check symmetry and support.
Numerically symmetrize only to remove rounding error, for example replacing $X$ by $(X+X^T)/2$. Check positive definiteness of $X$ and $I_p-X$ for type I, or of $F$ for type II.

The support checks have a structural explanation. For type I, $X=L^{-1}AL^{-T}$ is positive definite because it is a nonsingular congruence of $A$, and $I_p-X=L^{-1}BL^{-T}$ is positive definite for the same reason. For type II, $F=L^{-1}AL^{-T}>0$. These checks apply to exact arithmetic; with floating-point computations, a nearly singular factor can produce small numerical eigenvalues, so use a scale-aware tolerance rather than interpreting tiny rounding errors as a mathematical failure.

Common Mistake

A tempting but incorrect shortcut is to divide the entries of $A$ by the entries of $A+B$ to make a type I sample. Entrywise division does not implement matrix whitening: it need not preserve symmetry, and it does not guarantee that both $X$ and $I_p-X$ are positive definite. The correct operation is the congruence $L^{-1}AL^{-T}$, where $LL^T=A+B$. Likewise, type II uses $L^{-1}AL^{-T}$ with $LL^T=B$, not $B^{-1}A$, which is generally nonsymmetric. The congruence form makes the support guarantees explicit and is numerically preferable to computing a matrix inverse directly.

Key takeaway: Generate two independent nonsingular Wishart factors, then whiten the numerator by the sum for beta type I or by the denominator for beta type II. Cholesky congruences preserve symmetry and make the positive-definite support checks transparent.

Check Your Understanding

Answer these questions by referring to the transformations and support conditions.

For a $p\times p$ Bartlett factor with degrees of freedom $m$, what is the smallest diagonal chi-square degree of freedom, and what condition ensures it is positive?
If $S=A+B=LL^T$, write the Cholesky-based type I sample and the corresponding expression for $I_p-X$.
Why does a nonsingular congruence of a positive-definite matrix remain positive definite?
State the support check required for a type II sample.
Explain why $B^{-1}A$ is not the recommended type II sampling transform.