' $
Department of Mathematics
Recursive Formula for
E[
∏k
i=0
T r
{(W Σ
−1
)
m
i
}], Where
W
∼ W
p
(I, n) in Finite and Asymptotic
Regime
Jolanta Pielaszkiewicz, Dietrich von Rosen and Martin Singull
LiTH-MAT-R--2015/04--SE
Department of Mathematics Link¨oping University S-581 83 Link¨oping, Sweden.
Recursive formula for E
Q
ki=0Tr{(W Σ
−1)
mi},
where W ∼ W
p(Σ, n) in finite and asymptotic regime
Jolanta Pielaszkiewicz∗, Dietrich von Rosen,∗ and Martin Singull∗
∗Department of Mathematics,
Link¨oping University, SE–581 83 Link¨oping, Sweden. E-mail: jolanta.pielaszkiewicz@liu.se
E-mail: martin.singull@liu.se
Department of Energy and Technology,
Swedish University of Agricultural Sciences, SE–750 07 Uppsala, Sweden.
E-mail: Dietrich.von.Rosen@slu.se
Abstract
In this paper, we give a general recursive formula for E[Qk
i=0Tr{Wmi}], where
W ∼ Wp(I, n) denotes a real Wishart matrix. Formulas for fixed n, p are presented
as well as asymptotic versions when np n,p→∞→ c, i.e., when the so called Kolmogorov condition holds. Finally, we show application of the asymptotic moment relation when deriving moments for the Marchenko-Pastur distribution (free Poisson law). A numer-ical illustration using implementation of the main result is also performed.
Keywords: Wishart matrix; spectral distribution; eigenvalue distribution; moments; random matrices; free Poisson law; Marchenko-Pastur law.
1
Introduction
Multivariate analysis and random matrix theory are useful tools in mathematics, financial mathematics, statistics, engineering, physics as well as other disciplines too. The matrix dis-tribution which is nowadays known as the Wishart disdis-tribution was first derived by Wishart (1928). It is usually regarded as an extension of the chi-square distribution. Wishart ma-trices are, commonly used in statistics. For example, under a multivariate Gaussian model (normality), the distribution of the sample-covariance matrix is Wishart distributed. More-over, many test statistics considered in classical multivariate analysis are given as a function of one or several Wishart matrices. In this work the expectation of the product of traces of
Wishart matrices is studied. These quantities can be used when approximating densities. For example, Kawasaki and Seo (2012) considered tests for mean vectors with unequal co-variance matrices and then used moments of trace products of Wishart matrices. Moreover, when expanding the Stieltjes transform moments of products of traces appear.
Let X ∈ Rp×nfollow the central matrix normal distribution, denoted X ∼ Np,n(0, Σ, In),
where the dispersion matrix Σ in this work is assumed to be positive definite, denoted Σ > 0, and In is the identity matrix of size n × n. Alternatively, one can think about
a set of n independently distributed p-dimensional column vectors Xi, i = 1, . . . , n, each
distributed according to a multivariate normal distribution, i.e., Np(0, Σ). Then, W =
XX0 = Pn
i=1XiXi0, where X = (X1, . . . , Xn) : p × n and X0 denotes the transpose of X,
follows a central real Wishart distribution, W ∼ Wp(Σ, n). In fact, if there exist a normally
distributed X, with mean 0, and W = XX0 then W is said to be Wishart distributed. Another way of defining the Wishart distributions is to utilize the Laplace transform which provides a somewhat more general definition of the Wishart distribution.
Let for an arbitrary matrix A the matrix Akdenote
k times
z }| {
AA · · · A where usual matrix multi-plication is applied. Our main result provides a new recursive, non–combinatorial, formula for E[ k Y i=0 Tr{Wmi}], mi ≥ 0, i = 0, . . . , k, (1)
where W ∼ Wp(I, n), E[·] denotes expectation and the trace Tr{·} is defined as the sum of
the diagonal elements of a square matrix. The special case E[Tr{Wl}] coincides with the free moments of the empirical spectral distribution for W . Asymptotically (np → c > 0) delivers the moments of the free Poisson law, see Marchenko and Pastur (1967). Observe that since the free Poisson distribution is compactly supported the moments carry the full information about the distribution. Thus, it is interesting to put energy into the derivation of moment relations which can be interpreted and easily used in computations. Moreover, our expressions can be used to prove recursive formulas for the Catalan numbers. More about Catalan number can be found in Grimaldi (2011).
The new moment formula, given in this paper, is an extension of results provided by Fujikoshi, Ulyanov and Shimizu (2011), see Theorem 2.2.6, when W ∼ Wp(I, n). It is also
an extension of formulas given in Gupta and Nagar (2000). We will compare our results with formulas for the complex Wishart matrix presented in Haagerup and Thorbjornsen (2003) and Hanlon, Stanley and Stembridge (1992). The presented result provides also a formula for E[(Tr{W })k] as an alternative to the results given by Letac and Massam (2004), see further Section 2.
Furthermore, a relatively new idea when deriving moment relations is to apply so called umbra calculations. For example, see Di Nardo (2014) who studied relations similar to as well as more general than (1).
The paper is organized as follows. In Section 2 some previously obtained results for special cases of (1) are recalled, which indeed motivate this research. Section 3 comprises the derivation of the main result, i.e., a general recursive formula for (1). A number of corollaries will be presented as a consequences of the formula. Moreover, in Section 4
a theorem about asymptotic expressions is given. Then, in Section 5, we recall the free Poisson law and use our result to rediscover its moments. A recursive formula for Catalan numbers is also presented. At the end of the paper, in Section 6, a flowchart is given for computing moments explicitly together with a few numerical results.
2
A Survey of Known Results
In this section for both real and complex Wishart matrices a number of results that are of the same form as (1) will be considered.
2.1 Real Wishart matrix
Let W ∼ Wp(Σ, n). Then, one can easily see that Σ−1/2W Σ−1/2 ∼ Wp(I, n). Thus,
E[Tr{Σ−1W }] equals E[Tr{V }] if V ∼ Wp(I, n).
The expectation of the power of the trace of the Wishart matrix, W ∼ Wp(Σ, n),
has been derived by Nel (1971) for specific powers k. The general formula involve zonal polynomials and following Theorem 3.3.23 in Gupta and Nagar (2000) is given by
E[(Tr{W })k] = 2k X κ n 2 κ Zκ(Σ),
where Zκ(·) stands for zonal polynomial corresponding to κ, κ = (k1, . . . , km) is a partition
of k, such that ki ≥ 0 for all i = 1, . . . , m and Piki = k, [a]κ = Qmi=1[a − i + 1]ki with
[a]k = (a + k)!/a! for a ∈ C and k ∈ N0. The alternative version of the closed formula for
the expectation of the power of the trace of the Wishart matrix W can be found in Mathai (1980).
Using results for zonal polynomials, see e.g., Subrahmaniam (1976), E[(Tr{Σ−1W })k] = 2k np 2 k .
For any k ∈ N, a theorem stated in Letac and Massam (2004) gives an alternative formula E[(Tr{Σ−1W })k] =
X
(i)∈Ik
k!
i1! · · · ik!1i1· · · kik
(np)i1+...+ik2i2+2i3+...+(k−1)ik, (2)
where the set Ik consists of k-tuples (i) = (i1, . . . , ik) such that i1+ 2i2+ . . . + kik= k and
ij, j = 1, . . . , k are non-negative integers. Moreover, several other results and references
are mentioned in Gluckmuller (1998).
Above the expectation of the power of the trace of a real central Wishart matrix has been presented. Other formulations are about the expectation of the trace of the power, i.e., E[Tr{(Σ−1W )l}]. For a non-central Wishart distribution, V ∼ Wp(Σ, n, Θ), i.e., a Wishart
matrix which is defined via X ∼ Np,n(M, Σ, I) as V = XX0 with Θ = M M0, it is well
known that the following holds (see Gupta and Nagar (2000) for more details) E[Tr{V Σ−1}] = np + Tr{Θ},
Moreover, for W ∼ Wp(Σ, n) (see e.g., Fujikoshi, Ulyanov and Shimizu (2011))
E[Tr{W2}] = (n + n2)Tr{Σ2} + n(Tr{Σ})2. (4) Letac and Massam (2008) derived an expression for
E[
k
Y
i=0
Tr{(W Hi)}], (5)
where H1, . . . , Hk are arbitrary real symmetric matrices. Choosing Hi appropriately
mo-ments of all monomials are available. Using umbral calculations a generalization of (5) to a non-central Wishart matrix has been suggested by Di Nardo (2014).
Note that the notations in the papers of Letac and Massam (2008) and Di Nardo (2014) differ somewhat from ours since the Wishart distribution, e.g., in Letac and Massam is based on the Laplace transform with a different Σ.
2.2 Complex Wishart matrix
The complex p-dimensional column vector Xi follows the complex normal distribution Xi ∼
NC
p (µ, Σ) if the 2p-dimensional vector Yi= (<Xi0, =Xi0)0 ∼ N2p((<µ0, =µ0)0, ΣYi), where
ΣYi = 1 2 <Σ −=Σ =Σ <Σ ,
=Xi (resp. <Xi) denotes the imaginary (resp. real) part of the complex Xi and the matrix
=Σ is skew-symmetric.
An early reference to the complex Wishart distribution is Goodman (1963). Let again X = (X1, . . . , Xn) where Xi ∼ NpC(µ, Σ) which are supposed to be independently
dis-tributed. Then, a matrix WC ∼ WC
p(Σ, n) is complex Wishart distributed if WC = XX ∗,
where∗ denotes the conjugate transpose. An alternative definition of the complex Wishart distribution is based on the Laplace transform.
Let WC∼ WC
p(I, n). The matrix WC is considered in Hanlon, Stanley and Stembridge
(1992), where it was obtained that
E[Tr{WCk}] = 1 k k X j=1 (−1)j−1[n + k − j]k[p + k − j]k (k − j)!(j − 1)! , k ∈ N. (6) A corresponding recursive result was derived by Haagerup and Thorbjornsen (2003), i.e., for all k ∈ N the following holds
E[Tr{WC0}] = p, E[Tr{WC1}] = np, E[Tr{WCk+1}] = (2k + 1)(n + p) k + 2 E[Tr{W k C}] +(k − 1)(k 2− (n − p)2) k + 2 E[Tr{W k−1 C }]. (7)
The given recursive formula (7) was inspired by the Harer-Zagier recursion formula (see Harer and Zagier (1986)). The original paper Harer and Zagier (1986) gives a recursive relation for even moments of the spectral distribution of a complex self-adjoint random matrix Z = (Zij) formed of p2independent real standard normal distributed variables such
that Zij has mean 0 and variance 1, i.e., an expression for1pE[Tr{Z2k}]. It has been reproved and extended to the interesting case of Wishart matrices in Haagerup and Thorbjornsen (2003). The proof of (7) is mostly combinatorial while the proof given in Harer and Zagier (1986) combines advanced tools from analysis and theory of solving differential equations. Both the explicit result given in Hanlon, Stanley and Stembridge (1992) and the recursive result of Haagerup and Thorbjornsen (2003) derive E[Tr{WCk}], although equation (7) is said by Haagerup and Thorbjornsen (2003) to be more efficient than (6) for generation of moment tables.
2.3 Comparison of Moments for Complex and Real Wishart Matrices
The formulas regarding complex Wishart matrices WCand real Wishart matrices W differ for fixed n and p. In Table 1, E[Tr{W•k}] is presented for k ≤ 5, where W• denotes either
WC ∼ WC
p(I, n) or W ∼ Wp(I, n). There is a significant difference between the real and
complex cases. In particular, it means that (7) cannot be applied when WC is replaced by W .
Table 1: Comparison of E[Tr{Wk}] and E[Tr{Wk
C}], when k ≤ 5, WC ∼ WC,p(I, n) and
W ∼ Wp(I, n). The formulas have been derived using (7) and our result presented in
Theorem 3.1. k E[Tr{Wk}] E[Tr{WCk}] 1 np np 2 np(n + p + 1) np(n + p) 3 np(n2+ p2+ np + 3n + 3p + 4) np(n2+ p2+ np + 1) 4 np(n3+ p3+ 6n2p + 6np2+ 6n2+ 6p2+ 17np + 21n + 21p + 20) np(n3+p3+6n2p+6np2+5n+5p) 5 np(p4 + n4 + 10np3 + 10n3p + 20n2p2+ 61n2+ 61p2+ 163np + 10n3 + 10p3 + 55n2p + 55np2 + 148n + 148p + 132) np(n4 + p4 + 10n3p + 10np3 + 20n2p2+ 15n2+ 15p2+ 40np + 8) Using (7) E[Tr{WC2}] = 3(n + p) 3 E[Tr(WC)] = np(n + p),
which differs from the real case, i.e., E[Tr{W2}] = np(n + p + 1). We can also see the difference through a simple example.
be given as W = XX0 = X 2 11+ X122 X11X21+ X12X22 X11X21+ X12X22 X222 + X212 ! . Hence, since W ∼ W2(I2, 2), E[TrW ] = E[X112 + X122 + X222 + X212 ] = 4, E[Tr{W2}] = E[(X112 + X122 )2+ 2(X11X21+ X12X22)2+ (X222 + X212 )2] = 20,
which is consistent with the formulas (3) and (4). Similarly, let us consider the case with a 2 × 2 complex matrix X such that Xjk =
Yjk√+iZjk
2 , Yjk, Zjk ∼ N (0, 1) and Yjk, Zjk are
independent for all j, k = 1, 2. Then, Xjk has mean 0 and variance 1 and the matrix
WC= XX∗ ∼ WC 2(I2, 2). Furthermore, E[TrWC] = 1 2E[Y 2 11+ Z112 + Y122 + Z122 + Y212 + Z212 + Y222 + Z222 ] = 4, E[Tr{WC2}] = 1 4E (Y112 + Z112 + Y122 + Z122)2+ 2(Y12Y22+ Z12Z22+ Y11Y21+ Z11Z21)2 +2(Z12Y22− Y12Z22+ Z11Y21− Y11Z21)2+ (Y212 + Z 2 21+ Y 2 22+ Z 2 22) 2 = 16.
The difference between the results for the real and complex Wishart matrices of fixed size increases together with the power k. It can be observed that it is given by a polynomial of degree lower than k + 1.
3
A New Recursive Moment Formula
In this paper, the main goal is to generalize all the results for a real Wishart matrix which were presented in the previous by proposing a recursive formula for E[Qk
i=0Tr{(W Σ−1)mi}],
where W ∼ Wp(Σ, n) for any k ∈ N and mi ∈ N0, i = 1, . . . , k. Essential to the proof of
the forthcoming Theorem 3.1 is the use of the operator dXd , which can be considered as differentiation with respect to a symmetric matrix X. For Y ∈ Rq×r and X ∈ Rp×p
dY dX = X I ∂yij ∂xkl (gl⊗ gk)kl(ej⊗ di)0, kl= 1 : k = l, 1 2 : k 6= l, (8)
where I = {i, j, k, l : 1 ≤ i ≤ q, 1 ≤ j ≤ r, 1+ ≤ k ≤ p, 1 ≤ l ≤ n}, di, ej and gk are i-th,
j-th and k-th column of Iq, Ir and Ip, respectively, and ⊗ denotes the Kronecker product.
Rules of how to operate with (8) are given in A.
Moreover, the Wishart density function representing Wp(Σ, n) is used in the proof. Even
though the density exists only for n ≥ p, the result given below holds as well in the case p ≥ n, due to rotational invariance, i.e., Tr{Y Y0} = Tr{Y0Y }.
Theorem 3.1. Let W ∼ Wp(I, n). Then, the following recursive formula holds for all
k ∈ N and all m0, m1, . . . , mk such that m0= 0, mk∈ N, mi∈ N0, i = 1, . . . , k − 1
E k Y i=0 Tr{Wmi} = (n − p + mk− 1)E Tr{Wmk−1} k−1 Y i=0 Tr{Wmi} +2 k−1 X i=0 miE Tr{Wmk+mi−1} k−1 Y j=0 j6=i Tr{Wmj} + mk−1 X i=0 E Tr{Wi}Tr{Wmk−1−i} k−1 Y j=0 Tr{Wmj} . (9)
Remark 3.1. Note, that the presented formula is recursive with respect to the power mk.
Proof. Put mk = l + 1. Then,
Σ−1E k−1 Y i=0 Tr{(W Σ−1)mi}(W Σ−1)l =R Σ−1Qk−1 i=0Tr{(W Σ −1)mi}(W Σ−1)lf WdW, (10)
where fW denotes the density function for W , i.e., if Σ > 0 and n ≥ p
fW(W ) = 1 2pn2 Γp(n 2) |Σ|−n2|W | n−p−1 2 e− 1 2Tr{Σ −1W } , W > 0.
Note that as Tr{W Σ−1} is rotationaly invariant the result will also hold for p ≥ n.
On equation (10) we operate with dΣd−1 and then apply the trace Tr{·}. Note that LHS
(RHS) stands for left hand side (right hand side) of (10). Then, since E Qk−1
i=0 Tr{(W Σ−1)mi}(W Σ−1)l does not depend on Σ
d(LHS) dΣ−1 = dΣ−1 dΣ−1 E k−1 Y i=0 Tr{(W Σ−1)mi}(W Σ−1)l ⊗ Ip = 1 2(Ip2 + Kp,p) E k−1 Y i=0 Tr{(W Σ−1)mi}(W Σ−1)l ⊗ Ip ,
where Kp,p stands for the commutation matrix (see Appendix B or Kollo and von Rosen
(2005)). Moreover, Tr d(LHS) dΣ−1 = Tr 1 2(Ip2+ Kp,p) E k−1 Y i=0 Tr{(W Σ−1)mi}(W Σ−1)l ⊗ Ip = 1 2E k−1 Y i=0 Tr{(W Σ−1)mi}Tr{(W Σ−1)l}Tr{Ip} +1 2E k−1 Y i=0 Tr{(W Σ−1)mi}Tr{(W Σ−1)lIp}
= p + 1 2 E k−1 Y i=0 Tr{(W Σ−1)mi}Tr{(W Σ−1)l} . For the RHS the following holds
d(RHS) dΣ−1 = Z dfW dΣ−1vec 0 Σ−1 k−1 Y i=0 Tr{(W Σ−1)mi}(W Σ−1)ldW | {z } =B + Z d(Qk−1 i=0Tr{(W Σ−1)mi}Σ−1(W Σ−1)l) dΣ−1 fWdW | {z } =C . Since dfW dΣ−1 = 1 2(nvecΣ − vecW )fW, where (19) and (21) has been applied
B = E 1
2(nvecΣ − vecW )vec
0 Σ−1 k−1 Y i=0 Tr{(W Σ−1)mi}(W Σ−1)l and thus Tr{B} = n 2E k−1 Y i=0 Tr{(W Σ−1)mi}Tr{(W Σ−1)l} − 1 2E k−1 Y i=0 Tr{(W Σ−1)mi}Tr{(W Σ−1)l+1}.
Using the product rule (19) we have that C = C1+ C2, where
C1 = E d(Σ −1(W Σ−1)l) dΣ−1 k−1 Y i=0 Tr{(W Σ−1)mi} = E d(W −1(W Σ−1)l+1) dΣ−1 k−1 Y i=0 Tr{(W Σ−1)mi} (15) = E d((W Σ −1)l+1) dΣ−1 (I ⊗ W −1) k−1 Y i=0 Tr{(W Σ−1)mi} (20) = E d(W Σ−1) dΣ−1 X i+j=l (W Σ−1)i⊗ (W Σ−1)0j (I ⊗ W−1) k−1 Y i=0 Tr{(W Σ−1)mi} (15) = E dΣ −1 dΣ−1(I ⊗ W ) X i+j=l i,j≥0 (W Σ−1)i⊗ (W Σ−1)0j (I ⊗ W−1) k−1 Y i=0 Tr{(W Σ−1)mi} (12) = E Ip2 + Kp,p 2 X i+j=l i,j≥0 (W Σ−1)i⊗ W (W Σ−1)0jW−1 k−1 Y i=0 Tr{(W Σ−1)mi}
= 1 2E X i+j=l i,j≥0 (W Σ−1)i⊗ W (W Σ−1)0jW−1 k−1 Y i=0 Tr{(W Σ−1)mi} +1 2E Kp,p X i+j=l i,j≥0 (W Σ−1)i⊗ W (W Σ−1)0jW−1 k−1 Y i=0 Tr{(W Σ−1)mi} and C2 = E d(Qk−1 i=0 Tr{(W Σ−1)mi}) dΣ−1 vec 0 (Σ−1(W Σ−1)l) . Hence, Tr{C1} = 1 2 X i+j=l i,j≥0 E Tr{(W Σ−1)i}Tr{W (W Σ−1)0jW−1} k−1 Y i=0 Tr{(W Σ−1)mi} +1 2 X i+j=l i,j≥0 E Trn(W Σ−1)iW (W Σ−1)0jW−1o k−1 Y i=0 Tr{(W Σ−1)mi} = 1 2 X i+j=l i,j≥0 E Tr{(W Σ−1)i}Tr{(W Σ−1)0j} k−1 Y i=0 Tr{(W Σ−1)mi} +1 2 X i+j=l i,j≥0 E Tr{(W Σ−1)l} k−1 Y i=0 Tr{(W Σ−1)mi} = 1 2 X i+j=l i,j≥0 E Tr{(W Σ−1)i}Tr{(W Σ−1)0j} k−1 Y i=0 Tr{(W Σ−1)mi} +l + 1 2 E Tr{(W Σ−1)l} k−1 Y i=0 Tr{(W Σ−1)mi} .
Using again the product rule (19) we have
C2 = E d( Qk−2 i=0 Tr{(W Σ −1)mi}) dΣ−1 Tr{(W Σ −1)mk−1}vec0(Σ−1(W Σ−1)l) | {z } =D1 + E d(Tr{(W Σ −1)mk−1}) dΣ−1 k−2 Y i=0 Tr{(W Σ−1)mi}vec0(Σ−1(W Σ−1)l) | {z } =D2 .
The integral D2is computed using the chain rule (18) dΣd−1 = dW Σ −1
dΣ−1 dW Σd−1, and then (12),
(16) and (20), which gives
D2 = E d(Tr{(W Σ −1)mk−1}) dΣ−1 k−2 Y i=0 Tr{(W Σ−1)mi}vec0(Σ−1(W Σ−1)l) (18) = E d((W Σ −1)mk−1) dΣ−1 d(Tr{(W Σ−1)mk−1}) d((W Σ−1)mk−1) k−2 Y i=0 Tr{(W Σ−1)mi}vec0(Σ−1(W Σ−1)l) (18) = E d(W Σ −1) dΣ−1 d((W Σ−1)mk−1) d(W Σ−1) d(Tr{(W Σ−1)mk−1}) d((W Σ−1)mk−1) k−2 Y i=0 Tr{(W Σ−1)mi} vec0(Σ−1(W Σ−1)l) (16) (20) = E dΣ−1 dΣ−1(I ⊗ W ) X i+j=mk−1−1 i,j≥0 (W Σ−1)i⊗ (W Σ−1)0j vecI · k−2 Y i=0 Tr{(W Σ−1)mi}vec0(Σ−1(W Σ−1)l) (12) = E 1 2(Ip2 + Kp,p) X i+j=mk−1−1 i,j≥0 (W Σ−1)i⊗ W (W Σ−1)0j vecI · k−2 Y i=0 Tr{(W Σ−1)mi}vec0(Σ−1(W Σ−1)l) . Then, Tr{D2} = 1 2 X i+j=mk−1−1 i,j≥0 E Tr Σ−1(W Σ−1)l(W Σ−1)iI(W Σ−1)jW k−2 Y i=0 Tr{(W Σ−1)mi} +1 2 X i+j=mk−1−1 i,j≥0 E Tr Σ−1(W Σ−1)lW (W Σ−1)0jI(W Σ−1)0i k−2 Y i=0 Tr{(W Σ−1)mi} = 21 2mk−1E Tr{(W Σ−1)l+mk−1} k−2 Y i=0 Tr{(W Σ−1)mi} = mk−1E Tr{(W Σ−1)l+mk−1} k−2 Y i=0 Tr{(W Σ−1)mi} .
The integral D1 is of similar form as C2, hence repeating the same calculations we obtain Tr{C2} = k−1 X i=0 miETr{(W Σ−1)l+mi} k−1 Y j=0 j6=i Tr{(W Σ−1)mj}. Finally, we have p + 1 2 E k−1 Y i=0 Tr{(W Σ−1)mi}Tr{(W Σ−1)l} = n 2E k−1 Y i=0 Tr{(W Σ−1)mi}Tr{(W Σ−1)l} −1 2E k−1 Y i=0 Tr{(W Σ−1)mi}Tr{(W Σ−1)l+1} +1 2 X i+j=l i,j≥0 ETr{(W Σ−1)i}Tr{(W Σ−1)0j} k−1 Y j=0 Tr{(W Σ−1)mj} +l + 1 2 ETr{(W Σ −1)l} k−1 Y i=0 Tr{(W Σ−1)mi} + k−1 X i=0 miETr{(W Σ−1)l+mi} k−1 Y j=0 j6=i Tr{(W Σ−1)mj},
which is equivalent to the statement of the theorem.
If in Theorem 3.1 k = 1 and m1 = l + 1 then the next corollary is obtained.
Corollary 3.1. Let W ∼ Wp(I, n), then for all l ∈ N0
E[Tr{Wl+1}] = p[n − p − 1]l+1+ l X i=0 [n − p + i]l−i i X j=0 E[Tr{Wj}Tr{Wi−j}].
The expectation ETr{Wk}Tr{Wm} in Corollary 3.1 can be determined by for example
choosing k = 2, m1= m and m2= v in Theorem 3.1.
Corollary 3.2. Let W ∼ Wp(I, n). Then, for all v ∈ N and m ∈ N0,
ETr{Wv}Tr{Wm} = p[n − p − 1]vE[Tr{Wm}] +2m
v−1
X
i=0
[n − p + i]v−1−iE[Tr{Wm+i}]
+ v−1 X i=0 [n − p + i]v−1−i i X j=0 E[Tr{Wj}Tr{Wi−j}Tr{Wm}].
In order to obtain an expression for E[Tr{Wl+1}] it was shown that E[Tr{Wj}Tr{Wi−j}]
was needed where i ≤ l as well as j ≤ l. Thus, we gain in power, l instead of l + 1 but instead we got the new problem of finding the expectation of the product of two terms. In Corollary 3.2 the product was considered. Once again we gain in power but from now on expectation of a product of three terms has to be derived. This process can be continued and at the end the following result is needed:
Corollary 3.3. Let W ∼ Wp(I, n), then the following recursive formula holds for all t ∈ N
E(Tr{W })t+1 = np + 2tE(Tr{W })t = (np + 2t)!! (np)!! np.
Proof. The statement of Corollary 3.3 is a special case of Theorem 3.1, when mi = 1 for
all i = 1, . . . , k applied t + 1 times for k = t + 1, . . . , 1.
The chain of moment relations initiated in Corollary 3.1 and ended in Corollary 3.3 shows that any moment of the form E[Tr{Wl}] can be expressed explicitly.
4
Asymptotics of
pk+11nsE
Q
ki=0
Tr{W
mi}
Let W ∼ Wp(I, n). Then, E[Wl] = O(nl) which will serve as a basis for the forthcoming
theorem. As the considered expectation is increasing together with n and p we are interested to understand the asymptotic behavior of its normalized version pk+11nsE
Qk
i=0Tr{Wmi},
where s =Pk
i=0mi and when np→c while n, p → ∞. This means that p and n increase with
the same speed and the criterion is sometimes called the Kolmogorov condition (asymp-totic).
Theorem 4.1. Let W ∼ Wp(I, n), m = {m0, . . . , mk−1}, limnp = c > 0 and
(k)Q(mk, m) := lim n,p→∞ 1 pk+1nsE k Y i=0 Tr{Wmi}.
Then, for all k ∈ N and all mk, m such that m0 = 0, mk∈ N, mi∈ N0 for i = 1, . . . , k − 1
and s =Pk
i=0mi the following holds
(k)Q(mk, m) = (1c+ 1)(k)Q(mk− 1, m) +1cPmk−2
i=1 (k+1)Q(i, {m0, . . . , mk−1, mk− 1 − i}), mk> 2, k ≥ 1,
(1c+ 1)(k)Q(1, m), mk= 2, k ≥ 1, (k−1)Q(mk−1, {m0, . . . , mk−2}), mk= 1, k > 1,
Proof. Let us first notice that the LHS of (9) is a polynomial in n and p of degree k + 1 + s, similarly to the first and third summand on the RHS. The second summand on the RHS is of lower degree. Hence, when p, n → ∞, np → c,
1 pk+1nPki=0mi k−1 X i=0 miETr{Wmk−1+mi} k−1 Y j=0 j6=i Tr{Wmj} → 0. Therefore, (k)Q(mk, m) = 1 −1 c lim n,p→∞ 1 pk+1nPki=0mi−1 ETr{Wmk−1} k−1 Y i=0 Tr{Wmi} +1 cn,p→∞lim X i+j=mk−1 i,j≥0 1 pk+2ni+j+Pk−1 i=0miE Tr{Wi}Tr{Wj} k−1 Y i=0 Tr{Wmi} | {z } =:(k+1)Q(i,{m0,...,mk−1,j}) = (1c + 1)(k)Q(mk− 1, m) +1cPmk−2
i=1 (k+1)Q(i, {m0, . . . , mk−1, mk− 1 − i}), mk> 2, k ≥ 1,
(1c + 1)(k)Q(1, m), mk= 2, k ≥ 1,
(k−1)Q(mk−1, {m0, . . . , mk−2}), mk= 1, k > 1,
1, mk= 1, k = 1.
These limiting results allow us to analyze the asymptotics of the normalized expectation
1 pk+1nsE
Qk
i=0Tr{Wmi}.
5
Moments of the Free Poisson Distribution
A formula for the asymptotic empirical spectral distribution of M = 1nΣ−1/2W Σ−1/2, where W ∼ Wp(Σ, n), has been given in the late 60’s by Marchenko and Pastur (1967). Since
then it has been of constant use due to the development of Random Matrix theory and due to the fact that it arises in Free Probability theory established in Voiculescu (1985, 1991). Frequently the asymptotic empirical spectral distribution of M is called the Marchenko– Pastur law or alternatively, as an analog of the Poisson distribution, the free Poisson law. The support of the distribution is compact and its moments (free moments) are defined as 1pE[Tr{Mk}], k = 1, 2, . . .. Moreover, the free cumulants can be derived by one of the free cumulant-moment relations, see classical results in Nica and Speicher (2006), and a recursive alternative in Pielaszkiewicz, von Rosen and Singull (2014).
Note that the formula for the free moments of M , while p, n → ∞, np → c, satisfy the asymptotic equation derived in Section 4, i.e., 1pE[Tr{Ml}] →(1)Q(l, {∅}), where ∅ denotes
the empty set. Thus, lim p,n→∞ 1 pE[TrM ] = 1, lim p,n→∞ 1 pE[Tr{M 2}] = 1 +1 c, lim p,n→∞ 1 pE[Tr{M l+1}] = 1 +1 c 1 pp,n→∞lim E[Tr{M l}] (11) +1 c l−1 X i=1 lim p,n→∞ 1 p2E[Tr{M i}Tr{Ml−i}].
Moreover, the free moments can be presented non-recursively: lim p,n→∞ 1 pE[TrM ] = 1, lim p,n→∞ 1 pE[Tr{M 2}] = 1 +1 c, lim p,n→∞ 1 pE[Tr{M 3}] = 1 +1 c 2 +1 c = 1 + 3 c + 1 c2, lim p,n→∞ 1 pE[Tr{M 4}] = 1 +1 c 1 +3 c + 1 c2 + 21 c 1 +1 c = 1 + 6 c + 6 c2 + 1 c3.
These relations are the first four free moments of the free Poisson law and the results are in agreement with Oravecz and Petz (1997), where it was stated that
lim p,n→∞ 1 pE[Tr{M k}] = 1 k k X i=1 k i k i − 1 ci.
In the case of c = 1 the moments of the Marchenko-Pastur distribution are given by the Catalan numbers Ck := k+11 2kk. Catalan numbers have various applications, for example
in graph theory as a number of all rooted binary trees with k nodes or number of all ways to triangulate a regular k + 2 sided polygon. Hence, using (11) or computing directly from the asymptotic result from Section 4 we have that
Ck+1 = 2Ck+ k−1 X i=1 CiCk−i = k X i=0 CiCk−i,
which agrees with a recursive result for Catalan numbers. Interesting properties and ap-plications of Catalan numbers, along with results for Fibonacci numbers, can be found in Grimaldi (2011).
6
Implementation of results
The recursive formula given in Theorem 3.1 is a convenient tool for calculating non-asymptotic free moments of random Wishart matrices of fixed size, as well as more general
moment expressions such as E[Qk
i=0Tr{Wmi}] and its normalized version
1 pk+1nPki=0mi E[ k Y i=0 Tr{Wmi}],
where W ∼ Wp(I, n). In Table 2 a few numerical results are presented for m1 = 1, m2 =
4, m3 = 3.
Table 2: Numerical results for p31n8E[Tr{W1}Tr{W4}Tr{W3}] with given n and p such that
c = np = 2 and its asymptotics, when W ∼ Wp(I, n).
n 12 24 48 96 120 200 E[Tr{W1}Tr{W4}Tr{W3}] p3n8 41.03 21.99 17.71 16.38 16.17 15.86 n 400 800 2 · 103 2 · 104 2 · 105 2 · 106 E[Tr{W1}Tr{W4}Tr{W3}] p3n8 15.65 15.56 15.50 15.47 15.4691 15.46878 n → ∞ E[Tr{W1}Tr{W4}Tr{W3}] p3n8 c 5+9c4+25c3+25c2+9c+1 c5 = 15.46875
The asymptotic values are obtained by utilizing Theorem 4.1 in the following way:
(3)Q(3, {1, 4}) = 1 + 1 c (3)Q(2, {1, 4}) + 1 c(4)Q(1, {1, 4, 1}) = 1 + 1 c 2 (3)Q(1, {1, 4}) + 1 c(3)Q(1, {1, 4}) = 1 + 1 c 2 (2)Q(4, {1}) + 1 c(2)Q(4, {1}) = 1 +3 c + 1 c2 (2)Q(4, {1}) = 1 + 3 c + 1 c2 1 +1 c (2)Q(3, {1}) + 2 1 c(3)Q(1, {1, 2}) = 1 + 4 c + 4 c2 + 1 c3 1 +1 c (2)Q(2, {1}) + 1 c(3)Q(1, {1, 1}) + 2 c + 6 c2 + 2 c3 1 +1 c (2)Q(1, {1}) = 1 + 9 c + 25 c2 + 25 c3 + 9 c4 + 1 c5 c=2 = 15.46875.
Table 3 presents the moments pn1mE[Tr{W4}], where W ∼ Wp(I, n). The asymptotic value
is obtained using Theorem 4.1:
lim n,p→∞ 1 pn4E[Tr{W 4} = 1 +6 c + 6 c2 + 1 c3 e.g. = 45, c = 0.5 14, c = 1 2.448, c = 5 The case c = 1 corresponds to the Marchenko-Pastur distribution.
Table 3: Numerical results for 4th moment pn14E[Tr{W4}] for given finite n and p, such
that c = np ∈ {12, 1, 5}, and its asymptotic value.
n 6 60 100 400 103 104 → ∞
c = 12 57.51 46.08 45.65 45.16 45.06 45.01 45
1
pn4E[Tr{W4}], c = 1 20.09 14.50 14.29 14.07 14.03 14.00 14
c = 5 4.847 2.616 2.547 2.472 2.458 2.449 2.448
The values in Table 2 and 3 have been obtained using a function denoted Etr(m, n, p, Q = 1), which was implemented in R according to the algorithm given by flowchart in Figure 1. The blue states and dashed blue lines on flowchart indicate the recursive steps, diamonds stand for the decision states and the red parallelograms are the ending states. Yellow blocks indicate parts of algorithm responsible for calculating summands of formula (3.1). The two last summands are calculated using procedure given in the blocks on the right hand side of Figure 1. While the first summand as well as final summing of obtained results is implemented according to block on the left hand side of flowchart.
The algorithm depends on the following four input parameters: a) m - a k-dimensional vector of powers m1, . . . , mkwhich was used in Theorem 3.1; b) n and p - degrees of freedom
and size of the Wishart matrix, respectively; c) binary Q which describes if we are interested in the normalized expectation (Q = 1) or E Qk
i=0Tr{(W Σ−1)mi} (Q = 0).
Example 6.1. The function Etr((4), 120, 24, 1), counting numerical value of the fourth moment pn14E[Tr{(W Σ−1)4}] with n = 120, p = 24, is executed in R by
> Etr(c(4),120,24,1) [1] 2.530095
The forth moment E[Tr{(W Σ−1)4}] with n = 120, p = 24, can be obtained in R using function Etr((4), 120, 24, 0)
> Etr(c(4),120,24,0) [1] 12591371520
Value of p31n10E[Tr{(W Σ−1)4}Tr{W Σ−1}Tr{(W Σ−1)5}] with n = 120, p = 24, can be
obtained by
> Etr(c(4,1,5),120,24,0) [1] 10.6138
The code can be downloaded from http://www.mai.liu.se/~jolpi20/Etr.R.
Begin m, n, p, Q k = length(m) mk> 0 mk = mk− 1 k > 1 ¯ m = (m1, . . . , mk−1) temp = p · run( ¯m, n, p, Q) temp = 0 temp p mk> 1 i = 1 ¯ m = (m1, . . . , mk−1, i, mk− i)
tempq4[i] = run( ¯m, n, p, Q)
i = i + 1 i ≤ k − 1 temp4 = sum(tempq4) k > 1 i = 1 ¯ m = (m1, . . . , mi+ mk, . . . , mk−1)
tempq3[i] = 2mirun( ¯m, n, p, Q)
i = i + 1 i ≤ k − 1 temp3 = sum(tempq3) mk= 0 temp = n · run( ¯m, n, p, Q) temp = (n + p + mk)run( ¯m, n, p, Q)
temp + temp3 + temp4
YES NO YES NO NO YES YES YES NO NO YES YES NO NO YES NO E[Tr{(W Σ−1)0}QiTr{(W Σ −1)mi}] Pmk−1 i=1 E[Tr{(W Σ −1)i}Tr{(W Σ−1)mk−i}Q jTr{(W Σ −1)mj}] 2P imiE[Tr{(W Σ −1)mi+mk}Q jTr{(W Σ −1)mj}] (n + p + mk)E[Tr{(W Σ−1)mk}Q iTr{(W Σ−1)mi}] or nE[Tr{(W Σ−1)0}Q iTr{(W Σ −1)mi}]
Figure 1: Flowchart of implemented algorithm to calculate E[Qk
References
E. Di Nardo, On a symbolic representation of non-central Wishart random matrices with applications, J. Multivar. Anal. 125 (2014) 121–135.
Y. Fujikoshi, V. V. Ulyanov and R. Shimizu, Multivariate Statistics: High-Dimensional and Large-Sample Approximations (John Wiley & Sons, Hoboken, 2011).
D. H. Glueck and K. E. Muller, On the trace of a wishart, Comunnications in Statistics - Theory and Methods 27 (1998) 2137-2141.
N. R. Goodman, Statistical analysis based on a certain multivariate complex Gaussian distribution (an introduction), Ann. Math. Statist. 34 (1963) 152–177.
R. Grimaldi, Fibonacci and Catalan Numbers : An Introduction (John Wiley & Sons, Somerset, 2011).
A. K. Gupta and D. K. Nagar, Matrix Variate Distributions, Monographs and surveys in pure and applied mathematics; 104 (Chapman & Hall/CRC, Boca Raton, 2000). U. Haagerup and S. Thorbjornsen, Random matrices with complex Gaussian entries, Expo. Math. 21 (2003) 293–337.
P. J. Hanlon, R. P. Stanley and J.R. Stembridge, Some combinatorial aspects of the spectral of normally distributed random matrices, Contemp. Math. 138 (1992) 151–174. J. L. Harer and D. Zagier, The Euler characteristic of the moduli space of curves, Invent. Math. 85(3) (1986) 457–485.
H. Hotelling, A generalized T test and measure of multivariate dispersion, Proc. Second Berkeley Symp. Math. Statist. Prob. (Univ. of California Press, Berkeley, 1951), 23–42. T. Kawasaki and T. Seo, A two sample test for mean vectors with unequal covariance matrices, Technical report 12-19, Hiroshima University (2012).
T. Kollo and D. von Rosen, Advanced Multivariate Statistics with Matrices (Springer, Dordrecht, 2005).
M. Ledoux, A recursion formula for the moments of the Gaussian orthogonal ensemble, Ann. Inst. H. Poincar ˜A c Probab. Statist. 45 No. 3 (2009) 754–769.
G. Letac and H. Massam, All invariant moments of the Wishart distribution, Scand. J. Stat 31 No. 2 (2004) 295–318.
G. Letac and H. Massam, The noncentral Wishart as an exponential family, and its moments, J. Multivar. Anal. 99 No. 1 (2008) 1393–1417.
V. A. Marchenko and L. A. Pastur, Distribution of eigenvalues in certain sets of random matrices, Mat. Sb. (N.S.) 72(114):4 (1967) 507–536.
A. M. Mathai, Moments of the trace of a noncentral wishart matrix, Commun. Stat. Theory Methods 9:8 (1980) 795–801.
D. G. Nel, The h-th moment of the trace of a noncentral Wishart matrix, S. Afr. Statist. J. 5 (1971) 41–52.
A. Nica and R. Speicher, Lectures on the Combinatorics of Free Probability (Cambridge University Press, Cambridge, 2006).
F. Oravecz and D. Petz, On the eigenvalue distribution of some symmetric random matrices, Acta Sci. Math. 63 (1997) 383–395.
J. Pielaszkiewicz, D. von Rosen and M. Singull, Cumulant-moment relation in free probability theory, to appear in Acta Comment. Univ. Tartu. Math. 18(2) (2014). R. Speicher, Multiplicative functions on the lattice of noncrossing partitions and free convolution, Math. Ann. 298 (1994) 611–628.
K. Subrahmaniam, Recent trends in multivariate normal distribution theory: On the zonal polynomials and other functions of matrix argument, Sankhya Ser. A 38 No. 3 (1976) 221–258.
D. Voiculescu, Symmetries of some reduced free product C∗–algebras, Operator algebras and their connections with topology and ergodic theory, in Proc. Conf., Bu¸steni/Rom., Lect. Notes Math. 1132 (1985) pp. 556–588.
D. Voiculescu, Limit laws for Random matrices and free products, Invent. math. 104 (1991) 201–220.
J. Wishart, The generalized product moment distribution in samples from a normal multivariate population, Biometrika 20A (1928) 32–52.
A
Properties of Operator
dXdThe operator defined in (8) has the following properties: dX dX = 1 2(I + Kp,p) for symmetric X, (12) d(aY ) dX = a dY dX, (13) d(Z + Y ) dX = dZ dX + dY dX, (14) d(AXB) dX = B ⊗ A 0, (15) d(AY B) dX = dY dX(B ⊗ A 0 ), (16) dTr{A0X} dX = vecA, (17) dZ dX = dY dX dZ dY, (18) dW dX = dY dX dW dY + dZ dX dW dZ , where W = W (Y (X), Z(X)), (19) dYn dX = dY dX X i+j=n−1 i,j≥0 Yi⊗ (Y0)j , (20) d|X| dX = |X|vec(X −1)0. (21) Proof. Eq. (12) dX dX = X I ∂xij ∂xkl (fl⊗ gk)kl(fj⊗ gi)0 X=X0 = X k,l k6=l ∂xkl ∂xkl (fl⊗ gk) 1 2(fl⊗ gk) 0+X k,l k6=l ∂xlk ∂xkl (fl⊗ gk) 1 2(fk⊗ gl) 0 +X k ∂xkk ∂xkk (fk⊗ gk)(fk⊗ gk)0 = 1 2 X k,l k6=l (fl⊗ gk)[(fl⊗ gk)0+ (fk⊗ gl)0] + X k (fk⊗ gk)(fk⊗ gk)0 = 1 2 X k,l (fl⊗ gk)[(fl⊗ gk)0+ (fk⊗ gl)0] = 1 2(I + Kp,p)
Eq. (13) d(aY ) dX = X I ∂(ayij) ∂xkl (fl⊗ gk)kl(ej⊗ di)0 = aX I ∂yij ∂xkl (fl⊗ gk)kl(ej⊗ di)0 = a dY dX Eq. (14) d(Z + Y ) dX = X I ∂(zij + yij) ∂xkl (fl⊗ gk)kl(ej⊗ di)0 = X I ∂zij ∂xkl (fl⊗ gk)kl(ej⊗ di)0+ X I ∂yij ∂xkl (fl⊗ gk)kl(ej ⊗ di)0 = dZ dX + dY dX Eq. (15). Equivalently dYdX can be given as
dY dX = dvec0Y dvecX, then d(AXB) dX = d =vec0X(B⊗A0) z }| { vec0(AXB) dvecX = dvec0X(B ⊗ A0) dvecX = dvec0X dvecX(B ⊗ A 0) = dX dX(B ⊗ A 0 ) = B ⊗ A0 Eq. (16) d(AY B) dX = dY dX d(AY B) dY = dY dX(B ⊗ A 0 ) Eq. (17) dTr{A0X} dX (18) T r{Y }=vec0IvecY
= d(A 0X) dX d(vec0Ivec(A0X)) d(A0X) = d(A 0X) dX
d(vec0(A0Xvec0I)) d(A0X) = d(A 0X) dX d(A0X) d(A0X)vecI (15) = (I ⊗ A)vecI = vecA
Eq. (18) dZ dX = X I ∂zij ∂xkl (fl⊗ gk)kl(ej⊗ di)0 = X i,j,k,l X m,n ∂zij ∂ymn ∂ymn ∂xkl (fl⊗ gk)kl(hn⊗ km)0(hn⊗ km)nm(ej⊗ di)0 = X i,j,k,l,m,n,o,p ∂zij ∂ymn ∂yop ∂xkl (fl⊗ gk)kl(hp⊗ ko)0(hn⊗ km)nm(ej⊗ di)0 = X k,l,o,p ∂yop ∂xkl (fl⊗ gk)kl(hp⊗ ko)0 X i,j,m,n ∂zij ∂ymn (hn⊗ km)nm(ej⊗ di)0 = dY dX dZ dY Eq. (19). As ∂wij ∂xgh =X m,n ∂ymn ∂xgh ∂wij ∂ymn +X m,n ∂zmn ∂xgh ∂wij ∂zmn
then by the same argument as in proof of chain rule 6. the statement is proven. Eq. (20) (mathematical induction). For n=1 by Eq. 5. dYdX = dXdY(I ⊗ I).
For n = 2 by Eq. 7. dY2 dX = dY dX(Y ⊗ I) + dY dX(I ⊗ Y 0) = dY dX X i+j=1 i,j≥0 Yi⊗ (Y0)j.
Let assume statement is true for n = k − 1. Then dYk dX = d(Y Yk−1) dX prop.7. = dY k−1 dX (I ⊗ Y 0) +dY dX(Y k−1⊗ I) ind.assump. = dY dX X i+j=k−2 i,j≥0 Yi⊗ (Y0)j(I ⊗ Y0) +dY dX(Y k−1⊗ I) = dY dX X i+j=k−2 i,j≥0 Yi⊗ (Y0)j+1 + dY dX(Y k−1⊗ I) = dY dX X i+j=k−1 i,j≥0 Yi⊗ (Y0)j
Eq. (21). Determinant |X| and element of inverse matrix X−1 can be written as |X| = X j xij(−1)i+j|X(ij)|, (X−1)ij = (−1)i+1|X(ji)| |X| ,
where X(ij) denotes the minor of an element xij. Then as for all indeces i, j
∂|X| ∂xij
= (−1)i+j|X(ij)| = (X−1)ij|X|
the statement of (21) holds.
B
Commutation matrix K
p,pThe commutation matrix is defined as Kp,p =
X
I
eie0j⊗ eje0i,
where I = {i, j, k, l : 1 ≤ i ≤ p, 1 ≤ j ≤ p} and ei is i-th column of Iq. Let vec(•) denote
the usual vec-operator. Then, assuming sizes are appropriate,
Tr{AB} = vec0(A0)vec(B), (22)
Tr{A ⊗ B} = Tr{A}Tr{B}, (23)
