Matrix-Beta Normalizing Constants From Gamma Integrals

Why the Normalizing Constant Matters

A matrix-variate density is not a probability model until its integral over the full support is one. For the matrix-beta type I distribution, the density has a simple determinant kernel on the set of symmetric matrices between zero and the identity. Its normalizing constant is less immediate: it must account for all the matrix coordinates and their dependence. The matrix gamma integral provides a direct derivation, rather than requiring the constant to be guessed or borrowed from the scalar beta distribution.

The key idea is to start with two independent positive-definite matrix integrals, then change variables to their sum and a matrix-valued proportion. The resulting integral separates into a gamma integral and a matrix-beta integral. Keeping track of the Jacobian is essential: it is what produces the determinant power needed for the matrix gamma function to normalize the beta density.

The Matrix Gamma Integral

Let \(\mathcal S_p^+\) denote the set of real \(p\times p\) symmetric positive-definite matrices. For a symmetric matrix \(S\), \(|S|\) denotes its determinant and \(\operatorname{tr}(S)\) its trace. Integration over symmetric matrices uses Lebesgue measure on the \(p(p+1)/2\) distinct entries. For example, \(dS=\prod_{i\leq j}dS_{ij}\).

The condition \(a>(p-1)/2\) is needed for convergence at the boundary of the positive-definite cone. The integral can be derived explicitly by writing \(S=LL^{\mathsf T}\), where \(L\) is lower triangular with positive diagonal entries. The Cholesky change of variables has Jacobian \(dS=2^p\prod_{i=1}^p l_{ii}^{p+1-i}\prod_{i\geq j}dl_{ij}\). This Jacobian follows by induction on \(p\): after separating the last row of \(L\), the last off-diagonal block contributes \(|L_{p-1}|\) and the last diagonal coordinate contributes \(2l_{pp}\). Combining these with the dimension-\((p-1)\) Jacobian gives the displayed product.

Here is the integral calculation. Since \(|S|=\prod_i l_{ii}^2\) and \(\operatorname{tr}(S)=\sum_{i\geq j}l_{ij}^2\), each strict lower-triangular coordinate contributes \(\int_{-\infty}^{\infty}e^{-x^2}\,dx=\sqrt{\pi}\). There are \(p(p-1)/2\) such coordinates. For diagonal coordinate \(l_{ii}>0\), the determinant power and Jacobian together give exponent \(2a-i\), so its integral, including the factor 2 from the Jacobian, is \(2\int_0^\infty e^{-x^2}x^{2a-i}\,dx=\Gamma(a-(i-1)/2)\). The latter integral is finite precisely when \(2a-i>-1\); requiring this for every \(i\) gives \(a>(p-1)/2\). The nonnegative integrand permits the product separation by Tonelli's theorem. Multiplying the factors yields \(\Gamma_p(a)\), as claimed.

This derivation also explains why \(\Gamma_p\) is not simply an ordinary gamma function with a changed argument. Each diagonal coordinate has a different power, so the product contains gamma functions at successively shifted arguments.

Deriving the Matrix-Beta Integral

Define the matrix interval \(0<U<I\) to consist of symmetric matrices \(U\) for which both \(U\) and \(I-U\) are positive definite. For real \(a,b>(p-1)/2\), consider the integral of \(|U|^{a-(p+1)/2}|I-U|^{b-(p+1)/2}\) over this interval. The following derivation determines its value and establishes its finiteness under the stated conditions.

To prove the result, take the product of two matrix gamma integrals. For \(X,Y\in\mathcal S_p^+\), its integrand is \(|X|^{a-(p+1)/2}|Y|^{b-(p+1)/2}\exp\{-\operatorname{tr}(X+Y)\}\). The double integral equals \(\Gamma_p(a)\Gamma_p(b)\), by Tonelli's theorem.

Set \(S=X+Y\), and write \(S=LL^{\mathsf T}\) using its Cholesky factor. For each fixed \(S\), express \(X=LUL^{\mathsf T}\); then \(Y=L(I-U)L^{\mathsf T}\), with \(0<U<I\). These transformations give a one-to-one correspondence between \(X,Y\in\mathcal S_p^+\) and \(S\in\mathcal S_p^+\), \(0<U<I\). The Jacobian can be checked in two parts. First, \((X,Y)\mapsto(S,X)\), with \(S=X+Y\), has absolute Jacobian one. Second, at fixed \(S\), the congruence map \(U\mapsto LUL^{\mathsf T}\) on symmetric matrices has absolute Jacobian \(|L|^{p+1}=|S|^{(p+1)/2}\). To see the congruence formula, orthogonal changes of coordinates have absolute Jacobian one, and for a diagonal matrix with positive entries \(d_i\), the diagonal coordinates scale by \(d_i^2\) and the off-diagonal coordinates by \(d_id_j\). The product of these factors is \(\prod_i d_i^{p+1}\). A singular-value decomposition reduces general invertible \(L\) to this case. The full change of variables is block triangular in \((S,U)\), so its Jacobian is the same fixed-\(S\) factor.

Now \(|X|=|S||U|\), \(|Y|=|S||I-U|\), and \(\operatorname{tr}(X+Y)=\operatorname{tr}(S)\). Including the Jacobian, the power of \(|S|\) becomes \((a-(p+1)/2)+(b-(p+1)/2)+(p+1)/2=a+b-(p+1)/2\). Thus the transformed double integral factors as \[ \left(\int_{\mathcal S_p^+}|S|^{a+b-(p+1)/2}e^{-\operatorname{tr}(S)}\,dS\right) \left(\int_{0<U<I}|U|^{a-(p+1)/2}|I-U|^{b-(p+1)/2}\,dU\right). \] The first factor is \(\Gamma_p(a+b)\), since \(a+b>(p-1)/2\). Equating this product to \(\Gamma_p(a)\Gamma_p(b)\) proves the formula. Dividing the kernel by this integral gives a density that integrates to one. \(\square\)

Worked Constants and Density Checks

Why the Formula Is Useful—and a Common Pitfall

The multivariate gamma function converts an integral over the positive-definite cone into a closed-form constant. The beta identity then uses that same function to normalize a density over the bounded matrix interval. This is useful whenever a matrix-beta kernel appears in a model or derivation: the shape conditions and normalizing constant can be checked from one integral identity, while the dimension \(p\) remains explicit.

A common mistake is to substitute the scalar beta constant \(\Gamma(a)\Gamma(b)/\Gamma(a+b)\) directly, or to use only \(\Gamma_p(a)\Gamma_p(b)\) and omit the denominator. Neither gives the matrix-beta integral. The correct constant is a ratio of multivariate gamma functions, and the denominator comes from the integral over the sum matrix \(S\). A separate arithmetic pitfall is to simplify that ratio without checking the factors: in the \(p=2,a=2,b=3\) example, the numerator is \(3\pi^2/4\), and division by \(315\pi/2\) gives \(\pi/210\), not \(\pi/105\). The density constant must therefore be \(210/\pi\); using twice that kernel integral would make the proposed density integrate to \(1/2\).