Moments Of Matrix-Variate Beta Distributions

Why These Moments Have a Simple Form

A matrix-valued random variable contains several entries, so its mean might appear to require a separate calculation for every entry. For a standard matrix beta type I variable, symmetry reduces that task: its mean is a scalar multiple of the identity matrix. The scalar is determined by the expected trace, which can be derived directly from the density.

The conditions on the parameters matter. They first ensure that the matrix beta density is a probability density; once it is normalized, the bounded support ensures that its first moments, and indeed all polynomial moments, are finite. We will separate these two facts rather than infer existence from a formal algebraic expression for the mean.

Definition and Parameter Conditions

Let $X$ be a $p\times p$ real symmetric random matrix. Write $0<X<I_p$ to mean that both $X$ and $I_p-X$ are positive definite, where $I_p$ is the $p\times p$ identity matrix.

Definition: $X$ has a standard matrix-variate beta type I distribution with parameters $a,b$, written $X\sim \operatorname{Beta}_p^{\mathrm I}(a,b)$, if it has support $0<X<I_p$ and density, with respect to Lebesgue measure on the independent entries of a symmetric matrix, $$ f(X)=\frac{1}{B_p(a,b)} |X|^{a-(p+1)/2}|I_p-X|^{b-(p+1)/2},\qquad a,b>\frac{p-1}{2}. $$ Here $B_p(a,b)$ is the normalizing constant, and $|\cdot|$ denotes determinant.

The restrictions $a,b>(p-1)/2$ are the standard conditions for this density to be integrable and define a probability distribution. They are not optional qualifications on the moment formula: if either restriction fails, the standard matrix beta type I probability law is not defined by this density. For $p=1$, the density reduces to the ordinary beta density on $(0,1)$, with parameters $a$ and $b$.

Theorem: If $X\sim \operatorname{Beta}_p^{\mathrm I}(a,b)$ with $a,b>(p-1)/2$, then $$ E[X]=\frac{a}{a+b}I_p \qquad\text{and}\qquad E[\operatorname{tr}(X)]=\frac{pa}{a+b}. $$ Every polynomial moment of the entries of $X$ is finite under these conditions.

Why the Mean Must Be a Multiple of the Identity

The density is unchanged by orthogonal conjugation: for every orthogonal matrix $Q$, replacing $X$ by $QXQ^{T}$ preserves the support, determinant factors, and Lebesgue measure. Therefore $QXQ^{T}$ has the same distribution as $X$. The expectation exists because $X$ is bounded on its support, so $E[X]=QE[X]Q^{T}$ for every orthogonal $Q$.

A symmetric matrix unchanged by conjugation by every orthogonal matrix must be a scalar multiple of $I_p$. To see this, conjugation by diagonal sign-change matrices forces each off-diagonal entry to be zero, and conjugation by permutation matrices forces all diagonal entries to be equal. Thus $E[X]=cI_p$ for some scalar $c$. Taking traces gives $\operatorname{tr}(E[X])=pc$. It remains to calculate the expected trace.

Deriving the Expected Trace

Set $\alpha=a-(p+1)/2$ and $\beta=b-(p+1)/2$. Apart from its normalizing constant, the density is $g(X)=|X|^\alpha|I_p-X|^\beta$. Consider the symmetric matrix-valued function $V(X)=X(I_p-X)=X-X^2$. Its eigenvalue in an eigendirection of $X$, with eigenvalue $\lambda$, is $\lambda(1-\lambda)$; it vanishes at both ends of the interval $0<\lambda<1$. This makes it suitable for integration by parts on the support.

Derivation: On the vector space of symmetric matrices, use coordinates $x_{ii}$ and $x_{ij}$ for $i<j$, and define divergence as the sum of the derivatives of each coordinate component with respect to its matching coordinate. For $V=X-X^2$, $$ \operatorname{div}V=\frac{p(p+1)}{2}-(p+1)\operatorname{tr}(X). $$ Also, the directional derivative of $\log g(X)$ in direction $V(X)$ is $$ D_V\log g(X)=\alpha\operatorname{tr}\!\left(X^{-1}V\right) -\beta\operatorname{tr}\!\left((I_p-X)^{-1}V\right) =\alpha\bigl(p-\operatorname{tr}(X)\bigr)-\beta\operatorname{tr}(X). $$ Consequently, $$ \operatorname{div}\!\bigl(g(X)V(X)\bigr) =g(X)\bigl(ap-(a+b)\operatorname{tr}(X)\bigr). $$

For completeness, the divergence expression follows by differentiating the components. For a diagonal coordinate, $\partial V_{ii}/\partial x_{ii}=1-2x_{ii}$. For an off-diagonal coordinate $i<j$, $\partial V_{ij}/\partial x_{ij}=1-x_{ii}-x_{jj}$. Summing the diagonal terms gives $p-2\operatorname{tr}(X)$; summing the off-diagonal terms gives $p(p-1)/2-(p-1)\operatorname{tr}(X)$. Their sum is the stated divergence.

For the directional derivative, the determinant differential gives $D_V\log|X|=\operatorname{tr}(X^{-1}V)=p-\operatorname{tr}(X)$. Likewise, $D_V\log|I_p-X|=-\operatorname{tr}((I_p-X)^{-1}V)=-\operatorname{tr}(X)$. Combining these expressions with the divergence yields the displayed formula: the constant coefficient is $p(p+1)/2+p\alpha=pa$, and the coefficient of $\operatorname{tr}(X)$ is $-(p+1)-\alpha-\beta=-(a+b)$.

To integrate this identity, first restrict to matrices whose eigenvalues lie in $(\varepsilon,1-\varepsilon)$, and apply the divergence theorem there. The boundary flux vanishes as $\varepsilon$ decreases to zero. At a boundary eigenvalue $\lambda=\varepsilon$, the normal component of $V$ contributes a factor of order $\varepsilon$, while the density contributes $\varepsilon^\alpha$; the resulting boundary term has order $\varepsilon^{\alpha+1}$. Since $\alpha+1=a-(p-1)/2>0$, it tends to zero. At $\lambda=1-\varepsilon$, the corresponding factor is of order $\varepsilon^{\beta+1}$, which tends to zero because $b>(p-1)/2$. The other eigenvalue factors are integrable under the same parameter conditions. Thus the integral of the divergence over the full support is zero.

0=\int_{0<X<I_p}\operatorname{div}(gV)\,dX =\int_{0<X<I_p}g(X)\bigl(ap-(a+b)\operatorname{tr}(X)\bigr)\,dX.

Divide by the finite normalizing integral $\int g(X)\,dX$. This gives $ap-(a+b)E[\operatorname{tr}(X)]=0$, and hence $E[\operatorname{tr}(X)]=pa/(a+b)$. Since $E[X]=cI_p$, we have $pc=pa/(a+b)$, so $c=a/(a+b)$. Both claimed formulas follow. $\square$

Worked Calculations

Worked Example: A Two-Dimensional Mean Matrix

Let $X\sim\operatorname{Beta}_2^{\mathrm I}(3,5)$. First check the parameter conditions: $(p-1)/2=1/2$, and both 3 and 5 exceed $1/2$. The theorem therefore applies. The scalar multiplying the identity is $a/(a+b)=3/(3+5)=3/8$, so $$ E[X]=\frac{3}{8}I_2 =\begin{pmatrix}3/8&0\\0&3/8\end{pmatrix}. $$ The expected trace is $pa/(a+b)=2(3)/8=3/4$. This agrees with the trace of the mean matrix, $3/8+3/8=3/4$.

Worked Example: Expected Trace in Dimension Three

Suppose $Y\sim\operatorname{Beta}_3^{\mathrm I}(2,4)$. Here $(p-1)/2=1$, so $a=2>1$ and $b=4>1$; the density is valid. The expected mean matrix is $$ E[Y]=\frac{2}{2+4}I_3=\frac13 I_3. $$ The diagonal entries of the mean are each $1/3$, and its off-diagonal entries are zero. Therefore $E[\operatorname{tr}(Y)]=1/3+1/3+1/3=1$. The direct trace formula gives the same answer: $pa/(a+b)=3(2)/6=1$.

Worked Example: The Scalar Case

For $p=1$, take $Z\sim\operatorname{Beta}_1^{\mathrm I}(2,3)$. The parameter threshold is zero, so both parameters are allowed. The mean formula becomes $E[Z]=2/(2+3)=2/5$. Since a $1\times1$ matrix has trace equal to its sole entry, the expected trace is also $2/5$, matching $pa/(a+b)=1(2)/5$. This verifies that the matrix result reduces to the ordinary beta mean in dimension one.

Existence, Interpretation, and a Common Pitfall

Every eigenvalue of $X$ lies between zero and one. Consequently, the operator norm of $X$ is less than one, and each entry is bounded in absolute value by one. Any fixed polynomial in the entries is therefore bounded on the support and has a finite expectation under the probability density. The first-moment conditions do not impose additional restrictions beyond the conditions needed for the density to exist.

The expected trace also has a useful scale interpretation: it is the expected sum of the $p$ eigenvalues of $X$. Its value is $p$ times the mean eigenvalue in the sense of the average eigenvalue, $a/(a+b)$. This does not say that the eigenvalues are independent, nor does it determine their individual joint distribution. The mean matrix summarizes first moments, not the full dependence structure.

Pitfall: applying the formula outside the parameter domain. One might set $p=3$, $a=0.8$, $b=2$, and mechanically report $E[X]=0.8/2.8\,I_3$. The arithmetic is correct, but the conclusion is not: $(p-1)/2=1$, and $a=0.8$ fails the required condition. The standard matrix beta type I density is not normalizable for these parameters, so there is no such probability distribution whose mean is given by that expression.
Pitfall: treating the mean matrix as a covariance matrix. The mean is $a/(a+b)I_p$, not a covariance matrix or a statement that the entries of $X$ are independent. It describes the expected matrix itself; dependence between entries is not settled by this result.
Pitfall: forgetting the dimension in the trace. The expected trace is $pa/(a+b)$, not $a/(a+b)$. The latter is the multiplier of $I_p$; taking the trace adds that multiplier $p$ times.

Key takeaway: Under $a,b>(p-1)/2$, the standard matrix beta type I law is orthogonally invariant and supported on bounded matrices. Integration by parts gives $E[\operatorname{tr}(X)]=pa/(a+b)$; invariance then yields $E[X]=a/(a+b)I_p$. No stronger parameter conditions are needed for these moments.

Check Your Understanding

Use the support, parameter conditions, and moment identities to answer the following.

For $p=4$, what lower bound must both $a$ and $b$ strictly exceed for the standard density to be valid?
If $X\sim\operatorname{Beta}_2^{\mathrm I}(4,6)$, calculate $E[X]$ and $E[\operatorname{tr}(X)]$.
Why does orthogonal conjugation invariance imply that $E[X]$ is a scalar multiple of the identity?
Do the conditions for the density to exist impose additional restrictions for the first moments to be finite? Explain using the support.
For $p=3$, $a=0.9$, and $b=2$, can the standard matrix beta type I mean formula be used? Identify the condition that decides.