Distribution of the largest root of a matrix for Roy’s test in
multivariate analysis of variance

arXiv:1401.3987v4 [] 2 Apr 2017

Marco Chiania
a DEI, University of Bologna
V.le Risorgimento 2, 40136 Bologna, ITALY
marco.chiani@unibo.it

Abstract
Let X, Y denote two independent real Gaussian p × m and p × n matrices with
m, n ≥ p, each constituted by zero mean independent, identically distributed columns
with common covariance. The Roy’s largest root criterion, used in multivariate analysis of variance (MANOVA), is based on the statistic of the largest eigenvalue, Θ1 ,
−1
of (A + B) B, where A = XXT and B = YYT are independent central Wishart
matrices. We derive a new expression and efficient recursive formulas for the exact
distribution of Θ1 . The expression can be easily calculated even for large parameters,
eliminating the need of pre-calculated tables for the application of the Roy’s test.
Keywords: Roy’s test, Random Matrices, multivariate analysis of variance
(MANOVA), characteristic roots, largest eigenvalue, Tracy-Widom distribution,
Wishart matrices.
1. Introduction
The joint distribution of s non-null eigenvalues of a multivariate real beta matrix in the
null case can be written in the form [14, page 112], [2, page 331],
f (θ1 , . . . , θs ) = C(s, m, n)

s
Y

θim (1

i=1

n

− θi ) ·

s
Y

(θi − θj )

(1)

i<j

where 1 > θ1 ≥ θ2 · · · ≥ θs > 0, and C(s, m, n) is a normalizing constant given by

s
Y
Γ i+2m+2n+s+2
s/2
2


.
C = C(s, m, n) = π
(2)
Γ 2i Γ i+2m+1
Γ i+2n+1
2
2
i=1
This distribution arises in multivariate analysis of variance (MANOVA) and, with the
notation introduced above, is the distribution of the eigenvalues of (A + B)−1 B with
parameters
s = p,

m = (n − p − 1)/2,

n = (m − p − 1)/2 .

(3)

The marginal distribution of the largest eigenvalue, Θ1 , is of basic importance in testing hypotheses and constructing confidence regions in multivariate analysis of variance
Preprint submitted to JMVA

February 2014

(MANOVA) according to the Roy’s largest root criterion (see e.g. [2, page 333] and references therein), and is generally considered difficult to compute. For this reason, extensive studies have produced tables of upper percentage points for few specific (small)
values of s and some combinations of m, n (see e.g. [8, 15], [2, Table B.4]). The most
efficient numerical algorithm to compute the cumulative distribution function (CDF) of
Θ1 is provided in [3] based on [7], but for nonintegers m or n it requires infinite series
expansions which can results in computational time of several hours. Approximations
based on the Tracy-Widom distribution are discussed in [12].
In this paper we derive a simple expression for the exact CDF of Θ1 for arbitrary
s, n, m, and an iterative algorithm for its fast evaluation. The algorithm needs only
the incomplete beta function, and does not rely on numerical integration or series expansion. For instance, all results in [2, Table B.4] or in [12, Table 1] can be easily
computed instantaneously1 ; even for the most challenging cases analyzed in literature
[3, Table 1], which with previous methods required hours of computational time, no
more than a second is needed. Finally, we discuss some approximations based on the
Tracy-Widom distribution and its approximation [12, 4].
We remark that we study the case of real matrices: the complex analogous of our
problem, i.e., the case where X, Y are independent complex Gaussian, is much easier
and has been solved in [13].
Throughout the paper we indicate
R x with Γ(.) the gamma function, with B(a, b) the beta
function, with B (x; a, b) = 0 ta−1 (1 − t)b−1 dt the incomplete (lower) beta function
[1, Ch. 6], and with | · | the determinant.
2. Exact distribution of the largest eigenvalue for multivariate beta matrices in
the null case
The following is the main result of the paper.
Theorem 1. The CDF of the largest eigenvalue Θ1 for a multivariate beta matrix in
the null case is:
p
(4)
FΘ1 (θ1 ) = Pr {Θ1 ≤ θ1 } = C |A(θ1 )| .
When s is even, the elements of the s × s skew-symmetric matrix A(θ1 ) are:
ai,j (θ1 ) = E (θ1 ; m + j, m + i) − E (θ1 ; m + i, m + j)
where
E(x; a, b) ,

Z

i, j = 1, . . . , s

(5)

x

ta−1 (1 − t)n B (t; b, n + 1) dt.

(6)

0

When s is odd, the elements of the (s + 1) × (s + 1) skew-symmetric matrix A(x1 ) are
1 On

a current desktop computer in less than 0.1 seconds.

2

as in (5), with the additional elements
B (θ1 ; m + i, n + 1)

ai,s+1 (θ1 )

=

as+1,j (θ1 )

=

−aj,s+1 (θ1 )

as+1,s+1 (θ1 )

=

0

i = 1, . . . , s
j = 1, . . . , s

(7)
(8)
(9)

Note that ai,j (θ1 ) = −aj,i (θ1 ) and ai,i (θ1 ) = 0.
Moreover, the elements ai,j (θ1 ) can be computed iteratively, starting from the beta
function, without numerical integration or series expansion.
Proof. The proof is based on the approach introduced in [4] for Wishart and GOE
matrices.

Denoting ξ(x) = xm (1 − x)n , x = [x1 , x2 , . . . , xs ], and with V(x) = xi−1
the
j
Vandermonde matrix, we have for the eigenvalues in ascending order
f (xs , . . . , x1 ) = C

Y

(xj − xi )

i<j

s
Y

ξ(xi ) = C |V(x)|

i=1

where now 0 < x1 ≤ · · · ≤ xs < 1.
The CDF of the largest eigenvalue is then
Z
Z
FΘ1 (θ1 ) =
...

s
Y

ξ (xi )

(10)

i=1

f (xs , . . . , x1 )dx

(11)

0≤x1 <...<xs ≤θ1

Z
=

C

Z
|V(x)|

...
0≤x1 <...<xs ≤θ1

s
Y

ξ (xi ) dx .

(12)

i=1

To evaluate this integral we recall that for a generic s × s matrix Φ(x) with elements
{Φi (xj )} the following identity holds [6]
Z
Z
|Φ(x)| dx = Pf (A)
(13)
...
a≤x1 <...<xs ≤b

p

where Pf (A) = |A| is the Pfaffian, and the skew-symmetric matrix A is s × s for s
even, and (s + 1) × (s + 1) for s odd, with
Z bZ b
ai,j =
sgn(y − x)Φi (x)Φj (y)dxdy
i, j = 1, . . . , s.
(14)
a

a

Rb
For s odd the additional elements are ai,s+1 = −as+1,i = a Φi (x)dx, i = 1, . . . , s,
and as+1,s+1 = 0.
Using (13) in (12) with a = 0, b = θ1 , Φi (x) = xi−1 ξ(x) = xi−1 xm (1 − x)n with
some simple manipulations gives Theorem 1.
To avoid the numerical integration in (6), we first observe that the incomplete beta
functions can be computed iteratively by the relation
B(x; a + 1, b) =

a
xa (1 − x)b
B(x; a, b) −
.
a+b
a+b
3

(15)

R
This
relation is obtained for example with integration by parts of ta (1 − t)b−1 dt =
R
f (t)g 0 (t)dt with f (t) = ta − ta+1 , g 0 (t) = (1 − t)b−2 .
Then, from (6) and (15) the following identities can be easily verified:
E(x; a, a)

=

E(x; a, b + 1)

=

E(x; b, a)

=

1
2
B (x; a, n + 1)
2
E(x; a, b)
B(x; a + b, 2n + 2)
b
−
b+n+1
b+n+1
B (x; a, n + 1) B (x; b, n + 1) − E(x; a, b) .

(16)
(17)
(18)

Therefore, the elements of the matrix in (4) can be evaluated iteratively without any
numerical integration.
Theorem 1 can thus be used for an efficient computation of the exact CDF of the largest
eigenvalue for Roy’s test. In fact, using (16), (17), and (18), the CDF in (4) can be
simply evaluated, without numerical integrations or series expansion, by Algorithm 1
reported below.
Algorithm 1 Algorithm: CDF of the largest eigenvalue for Roy’s test
Input: s, m, n, θ1
Needed:
Incomplete beta function B (x; a, b)
Output: FΘ1 (θ1 ) = Pr {Θ1 ≤ θ1 }
A=0
for i = 1 → s do
bi = B(θ1 ; m + i, n + 1)2 /2
for j = i → s − 1 do
m+j
B(θ1 ; 2m + i + j, 2n + 2)
bj+1 =
bj −
m+j+n+1
m+j+n+1
ai,j+1 = B (θ1 ; m + i, n + 1) B (θ1 ; m + j + 1, n + 1) − 2bj+1
end for
end for
if s is odd then
Append one column to A with ai,s+1 = B (θ1 ; m + i, n + 1) , i = 1, . . . , s
Append one zero row to A
end if
A = A − AT
p
return FΘ1 (θ1 ) = C(s, m, n) |A|
Implementing directly the algorithm in Mathematica on a personal computer, for each
value x1 we obtain the exact CDF in (4) in less than 0.1 seconds for all tables in [15],
[2, Table B.4] and [12, Table 1]. For s = 54, m = −1/2, n = 45/2, which corresponds to the most challenging case analyzed in literature [3, Table 1, last row], the
exact CDF is computed with the new expression in less than a second. For comparison,
the approach based on series expansions of the hypergeometric function of a matrix argument requires hours. For larger parameters for which no other methods are available
4

in literature, like for example s = 200, m = −1/2, n = 299/2, the computational time
is less than fifteen seconds on a common personal computer. Therefore, the upper percentage points for Roy’s test can be the evaluated exactly and almost instantaneously
for parameters of interest in applied statistics.
3. Approximations
In [12] it is shown that the logit of Θ1 approaches the Tracy-Widom law for large s.
More precisely, it is shown that, for m ≥ −1/2, n ≥ 0 and large s, the following
approximation holds
log (Θ1 /(1 − Θ1 )) − µ D
≈ TW1
(19)
σ
where TW1 denotes a random variable with Tracy-Widom distribution of order 1 [16,
17, 9, 10, 18], and the values of µ, σ are given, using the parameters in (3), by [11, 12]


γ+φ
µ = 2 log tan
(20)
2
1
16
(21)
σ3 =
(m + n − 1)2 sin2 (γ + φ) sin γ sin φ


m + n − 2p
(22)
γ = arccos
m+n−1


m−n
φ = arccos
.
(23)
m+n−1
Moreover, in [4] it is shown that the TW1 can be closely approximated by a shifted
gamma distribution, so that


x+α
F1 (x) , Pr(TW1 ≤ x) ' P k,
x > −α
(24)
δ
x = F1−1 (y) ' δ P −1 (k, y) − α

(25)

where P (a, x) is the regularized lower incomplete gamma function, P −1 (a, y) is its
inverse, and the constants k = 46.446, δ = 0.186054, α = 9.84801 have been chosen
to match the moments of the approximation to that of the Tracy-Widom.
Thus, using (24) in (19) we obtain for the CDF of the Roy’s statistic in the null case


log (θ1 /(1 − θ1 )) − µ + σα
(26)
FΘ1 (θ1 ) ' P k,
δσ
and for its inverse, useful for evaluating the percentiles,

exp σ[δP −1 (k, y) − α] + µ
−1
θ1 = FΘ1 (y) '
.
1 + exp {σ[δP −1 (k, y) − α] + µ}

(27)

An example comparing exact (by Algorithm 1) with approximate (by (26)) distributions
is reported in Fig. 1.
5

CDFHΘ1 L
1

0.8

0.6

0.4

0.2

0
0

0.2

0.4

0.6

0.8

1

Θ1

Figure 1: CDF of the largest eigenvalue, real beta matrix, m = −1/2, n = 100. From left to right:
s = 5, 15, 100. Comparison between the exact distribution (4) (solid line) and the approximation in (26)
(dotted). Note that for s = 100 the two curves are almost indistinguishable.

While Theorem 1 with Algorithm 1 gives the exact distribution, the two approximations
above, which are in simple closed forms, can be used for a rapid, approximate design of
the Roy’s test. For instance, with s = 5, m = −1/2, n = 1000 the 80th-percentile obtained by root-finding with the exact CDF of Algorithm 1 is θ1 = 0.008501, while (27)
gives θ1 = 0.008609. For larger matrices and larger values of the CDF the approximation generally improves. For example, with s = 200, m = −1/2, n = 299/2 the 99thpercentile obtained by Algorithm 1 is θ1 = 0.827760, while (27) gives θ1 = 0.827761.
4. Remarks for the complex multivariate beta
For completeness we recall that, when X, Y are two independent complex Gaussian,
the analogous of (1) is the complex multivariate beta, where the joint distribution of
the eigenvalues is [13]
f (θ1 , . . . , θs ) = C 0 (s, m, n)

s
Y

θim (1 − θi )n ·

i=1

s
Y

2

(θi − θj )

(28)

i<j

where 1 > θ1 ≥ θ2 · · · ≥ θs > 0, and
C 0 (s, m, n) =

s
Y

Γ (m + n + s + i)
.
Γ
(i)
Γ (i + m) Γ (i + n)
i=1

In this case we can write
f (θ1 , . . . , θs ) = C 0 (s, m, n) |V(θ)|

2

s
Y
i=1

6

ξ(θi ) .

(29)

Therefore, by applying [5, Corollary 2] the CDF of Θ1 in the complex case is simply
given by [13]
FΘ1 (θ1 ) = Pr {Θ1 ≤ θ1 } = C 0 (s, m, n) |A(θ1 )|

(30)

where the elements of the s × s matrix A(θ1 ) are:
ai,j (θ1 ) = B(θ1 , m + i + j − 1, n + 1)

i, j = 1, . . . , s .

(31)

Moreover, expressions analogous to (24) and (27) can be derived, by using the TracyWidom limiting behavior for Θ1 in the complex case [11] and the approximation of the
Tracy-Widom of order 2 with a shifted gamma [4].
References
[1] M. Abramowitz, I. A. Stegun, Handbook of Mathematical Functions wih Formulas, Graphs, and Mathematical Tables, United States Department of Commerce,
Washington, D.C., 1970.
[2] T. W. Anderson, An Introduction to Multivariate Statistical Analysis, Wiley, New
York, 2003.
[3] R. W. Butler, R. L. Paige, Exact distributional computations for Roy’s statistic and
the largest eigenvalue of a Wishart distribution, Statistics and Computing 21 (2)
(2011) 147–157.
[4] M. Chiani, Distribution of the largest eigenvalue for real Wishart and Gaussian
random matrices and a simple approximation for the Tracy-Widom distribution,
Journal of Multivariate Analysis 129 (2014) 69 – 81, arXiv:1209.3394.
[5] M. Chiani, M. Z. Win, A. Zanella, On the Capacity of Spatially Correlated MIMO
Rayleigh Fading Channels 49 (10) (2003) 2363–2371.
[6] N. De Bruijn, On some multiple integrals involving determinants, J. Indian Math.
Soc 19 (1955) 133–151.
[7] R. D. Gupta, D. S. P. Richards, Hypergeometric functions of scalar matrix argument are expressible in terms of classical hypergeometric functions, SIAM J.
Math. Anal. 16 (1985) 852–858.
[8] D. Heck, Charts of some upper percentage points of the distribution of the largest
characteristic root, The Annals of Mathematical Statistics 31 (3) (1960) 625–642.
[9] K. Johansson, Shape Fluctuations and Random Matrices, Communications in
Mathematical Physics 209 (2000) 437–476.
[10] I. Johnstone, On the distribution of the largest eigenvalue in principal components
analysis, The Annals of Statistics 29 (2) (2001) 295–327.

7

[11] I. Johnstone, Multivariate analysis and Jacobi ensembles: Largest eigenvalue,
Tracy-Widom limits and rates of convergence, The Annals of Statistics 36 (8)
(2008) 2638–2716.
[12] I. M. Johnstone, Approximate null distribution of the largest root in multivariate
analysis, The annals of Applied Statistics 3 (4) (2009) 1616–1633.
[13] C. G. Khatri, Distribution of the largest or the smallest characteristic root under
null hypothesis concerning complex multivariate normal populations, Ann. Math.
Stat. 35 (1964) 1807–1810.
[14] R. J. Muirhead, Aspects of Multivariate Statistical Theory, Wiley, New York,
1982.
[15] K. S. Pillai, Upper percentage points of the largest root of a matrix in multivariate
analysis, Biometrika 54 (1-2) (1967) 189–194.
[16] C. Tracy, H. Widom, Level-spacing distributions and the Airy kernel, Communications in Mathematical Physics 159 (1) (1994) 151–174.
[17] C. Tracy, H. Widom, On orthogonal and symplectic matrix ensembles, Communications in Mathematical Physics 177 (1996) 727–754.
[18] C. Tracy, H. Widom, The distributions of random matrix theory and their applications, New Trends in Mathematical Physics (2009) 753–765.

8

