arXiv:1707.02564v1 [] 9 Jul 2017

Holonomic Gradient Method for the Distribution
Function of the Largest Root of Complex
Non-central Wishart Matrices
Fadil Habibi Danufane (Kanazawa University)
Katsuyoshi Ohara (Kanazawa University)
Nobuki Takayama (Kobe University)
Abstract
We give a new method to evaluate the cumulative distribution function of the largest root of complex non-central Wishart matrices. We
are motivated by a performance analysis in wireless communication systems; we evaluate the outage probability of a multiple-input-multipleoutput (MIMO) system employing maximal ratio combining (known also
as beamforming systems) and operating over Rician-fading channels, which
is expressed as the cumulative distribution function of the largest eigenvalue of a non-central Wishart matrix, by using a relatively new numerical
method called holonomic gradient method (HGM). As a validity check,
we compare our result with the computation using numerical integration.
It is shown by the numerical computation that the HGM gives a very
accurate result within a much shorter computation time.

1

Introduction

This paper proposes a new method to evaluate numerically the cumulative distribution function (CDF) of the largest root of s×s complex non-central Wishart
random matrices (see, e.g., [8] and its references). Our method is efficient and
accurate. The evaluation is a fundamental problem in statistics and has broad
applications. Among them, the main application in our mind is the performance
analysis of wireless communications. In conventional wireless communications,
a single antenna is used at the source and another single antenna is used at the
destination. In some cases, this gives rise to problems with multipath effects
which can cause a reduction in bandwidth and an increase in the number of
transmission errors. One of the solutions to this problem which is very popular
among telecommunication engineers is the use of two or more antennas, hence
the multiple-input-multiple-output (MIMO) system. The outage probability is
an important statistical measure to assess the quality of service provided by the
system. It is defined as the probability of failing to achieve a specified signal
to noise ratio (SNR) value sufficient for satisfactory reception. It is known that
the outage probability can be expressed as the CDF of the largest eigenvalue of
non-central Wishart matrices.
Our paper based on the two milestone papers on the analysis of MIMO
systems and a novel method in statistics that is called the holonomic gradient
1

method (HGM). The first paper is the determinant formula of the CDF derived
by Kang and Alouini [5]. They also gave numerical evaluations of the CDF with
their determinant formula under several restrictions to analyze MIMO systems
including beamforming systems operating over Rician-fading channels. The second paper is an analysis of zero forcing MIMO systems by Siriteanu, Takemura,
Kuriki, Shin, Koutschan [12]. They first applied the HGM to analyze the performance measures, such as outage probability and ergodic capacity, of the MIMO
system that employs zero forcing detection method under Rician fading. They
evaluated the outage probability by first finding differential equations satisfied
by the moment-generating function of the SNR and taking its inverse Laplace
transformation to obtain differential equations for the probability distribution
function of the SNR before finally numerically integrating the equations. The
HGM is an extrapolation method to numerically evaluate various CDF’s and
normalizing constants in statistics by utilizing a numerical analysis of ordinary
differential equations (ODE). See [3], [7], [11] on the HGM for some important
distributions in statistics and [15] for a comprehensive references on the HGM.
We refer to the book [4, Chapter 6] on some fundamental background on the
HGM.
Comparatively to [12], our paper aims to numerically evaluate the outage
probability of MIMO systems by using the determinant formula which has been
derived and studied theoretically by Kang and Alouini [5]. We will remove
restrictions on the numerical evaluation of the paper with the HGM. However,
we note that the application is more challenging than the case of [12] where
they analyze the zero forcing case and there is no instability of the HGM.
In the past years, there have been several attempts to compute the outage
probability by using methods such as Monte-Carlo or power series expansion
(see, e.g., references in [5]). However, those computations are relatively slow.
This disadvantage becomes more apparent especially when we want to perform
a computation over several threshold values or when a large number of antennas
are involved. Our application of the HGM can perform this computation for
several values and a large number of antennas efficiently and accurately.
The remainder of this paper is organized as follows. Section 2 presents the
known results derived in [5] regarding the formula of outage probability with
respect to a certain MIMO system model. We will see that computation of
the formula will reduce to an evaluation of a certain function that involves
hypergeometric series. Section 3 demonstrates the HGM application to the
hypergeometric function, utilizing a 4th order ODE that it satisfies. We will
also show that the computation by using this method is unstable at points
which are far from the initial value. Section 4 gives an alternative method of
HGM to compute the outage probability which corresponds to 3rd order ODE
instead of the 4th order one that we used previously. Finally, Section 5 shows
some results of numerical evaluation on different cases of antenna’s number or
eigenvalues by using the HGM and numerical integration as a comparison. It
also highlights the better performance of the HGM compared to the integration
in term of computation time for very similar accuracy.

2

2

Distribution of the largest root of complex
non-central Wishart matrices

Consider an antenna system which consists of NT and NR antennas at the
transmitter side and the receiver side respectively. This antenna system can be
modelled by using an NT × NR matrix H where the columns are independent
NT -variate complex Gaussian vectors with the covariance matrix Σ and the
mean E(H) = M .
Let s = min(NT , NR ), t = max(NT , NR ), and 0 < φ1 < φ2 < · · · < φs
be s nonzero eigenvalues of S = Σ−1 HH ∗ where the matrix H ∗ denotes the
Hermitian conjugate transpose of H. If NT ≤ NR , then S has complex Wishart
distribution. Kang and Alouini [5] show that if M ∗ Σ−1 M has s nonzero distinct eigenvalues 0 < λ1 < λ2 < · · · < λs , then the outage probability can be
expressed in terms of the cumulative distribution function (CDF) of the largest
eigenvalue φs evaluated at certain threshold value x, that is,
Pr (φs ≤ x) =

exp(−Λ)
det ψ(x).
V · (Γ(t − s + 1))s

(1)

Here, we put Λ = λ1 + · · · + λs ,
Y

V =

(λj − λi )

1≤i<j≤s

and ψ(x) is a square matrix of order s whose (i, j)-entries are given by
Z x
ψ(x)i,j =
y t−i e−y 0 F1 (; t − s + 1; y)dy

(2)

0

where 0 F1 (; a; x) is the confluent hypergeometric function defined as
0 F1 (; a; x)

=

∞
X

xp
(a)p (1)p
p=0

(3)

and (n)p = n(n + 1)(n + 2) · · · (n + p − 1) is the Pochhammer symbol.
For the sake of simplicity, let k = t − i and n = t − s + 1. Then we have
Z x
ψ(x)i,j =
y t−i e−y 0 F1 (; t − s + 1; λj y)dy
0
Z x
=
y k e−y 0 F1 (; n; λj y)dy
0

= Hnk (x, λj ).

(4)

By exchanging the integral and the sum, we obtain the following series expansion
of Hnk (x, λj ) where
Hnk (x, y) =

∞
(−1)q (k + 1)p+q
xk+1 X
xp+q y p .
k + 1 p,q=0 (n)p (1)p (1)q (k + 2)p+q

(5)

Now the formula for the outage probability looks like this:

Pr (φs ≤ x)

=

exp(−Λ)
V · (Γ(t − s + 1))s

3

Hnt−1 (x, λ1 ) · · ·
Hnt−2 (x, λ1 ) · · ·
..
..
.
.
t−s
Hn (x, λ1 ) · · ·

Hnt−1 (x, λs )
Hnt−2 (x, λs )
. (6)
..
.
Hnt−s (x, λs )

Given a threshold value x, we want to compute the value of outage probability, Pr (φs ≤ x). Our strategy is to evaluate the values of all of the function
Hnk (x, λj ) and compute the probability straightforwardly. However, due to its
slow convergence, each value of Hnk (x, λj ) is difficult to compute accurately and
efficiently by using power series expansion.

3

Local Numerical evaluation of Hnk (x, y) by 4th
order differential equation by using holonomic
gradient method

Recently, a relatively new and efficient method to do a numerical evaluation
of multivariate series such as Hnk (x, y) called the holonomic gradient method
(HGM) has been used and demonstrated effectiveness in various problems in
mathematics and statistics (see [3], [7], [11]). This method utilizes the fact that
any holonomic function satisfies linear differential equations with polynomial coefficients and exterpolate the function with the differential equations. Generally,
we can evaluate any holonomic function at some points by numerically solving
its differential equations starting from some initial value where the function is
either known analytically or can be approximated accurately. The HGM is a
new and efficient method that combines the theories of D-modules and Gröbner
basis to do a computation of a holonomic function by numerically solving differential equations satisfied by the holonomic function. The HGM consists of
three steps (see, e.g., [4, Chapter 6]): (1) Finding a holonomic system satisfied
by the function, (2) translating the holonomic system into a Pfaffian system and
finding an initial value, (3) numerically evaluating a certain function at several
points using the obtained Pfaffian system with the Runge-Kutta method (see,
e.g., [9]).

3.1

Holonomic system for Hnk (x, y)

Consider a complex-valued function f (x) of x = (x1 , . . . , xn ) ∈ Cn . We call
f (x) a holonomic function if and only if, for all i ∈ {1, . . . , n}, there exists the
following ODE of order Ni in xi -direction
!
∂ Ni
∂
hi,Ni (x) Ni + · · · + hi,1 (x)
+ hi,0 (x) • f (x) = 0
∂xi
∂xi
where hi,Ni (x), . . . , hi,0 (x) are polynomials in x1 , . . . , xn . We want to show
that the function Hnk (x, y) is holonomic by finding ODE’s in x-direction and
y-direction that it satisfies. Recall that
Hnk (x, y) =

∞
(−1)q (k + 1)p+q
xk+1 X
xp+q y p .
k + 1 p,q=0 (n)p (1)p (1)q (k + 2)p+q

We can compute a system of homogeneous linear partial differential equations
for such a power series whose coefficients consist of Pochhammer symbols by
using a well-known method (e.g. see [6] Section 12). The function Hnk (x, y)

4

satisfies the following system.
{θy (θy + n − 1) + y(θx − θy − k − 1)} • Hnk (x, y) = 0,
{θx (θx − θy − k − 1) + xθx } • Hnk (x, y) = 0.

(7)

Here, we put
θx = x

∂
,
∂x

θy = y

∂
.
∂y

Remember that we need to find ordinary ODE’s to prove the holonomicity
of Hnk (x, y). One way to achieve that is to use Gröbner basis of a left ideal on
a ring of differential operators (see, e.g., [4]).
∂
∂
, ∂y
i be the ring of differential operators whose coeffiLet R = C(x, y)h ∂x
cients are complex-valued rational functions of (x, y) and put
l1

= θy (θy + n − 1) + y(θx − θy − k − 1),

l2

= θx (θx − θy − k − 1) + xθx .

Clearly operators l1 and l2 are elements of the ring R. Denote by I the left
R-ideal generated by l1 and l2 . We identify the left ideal with the system (7).
We can compute Gröbner bases of the left ideal using Gröbner basis package
“yang.rr” for a computer algebra system Risa/Asir (see, e.g., [10] for reference).
It is well-known that Gröbner basis depends on a choice of the monomial
∂
∂
ordering. In our case, the ring R has two variables ∂x
and ∂y
. First, by choosing
∂
∂
the lexicographic ordering with ∂x < ∂y , we obtain a Gröbner basis {l3 , l4 , l5 }
where
l3

=

((θx + x − k − 1)(θx + x − k + n − 2) − xy)θx ,

l4

=

(θy − θx − x + k + 1)θx ,

l5

= θy2 − (y − n + 1)θy + yθx − y(k + 1).

Then we have that the holonomic rank of the system is 4.
Likewise, by choosing the lexicographic ordering with
the Gröbner basis {l6 , l7 } where
l6

(8)
∂
∂y

<

∂
∂x ,

we obtain

= θy4 + (y + 2n − 4)θy3 + (n2 − 5n + 5 − (x + k + n)y) θy2
+(xy 2 − (n − 1)xy − (k + 1)(n − 1)y − (n − 2)(n − 1))θy

l7

+(ky + 1)xy,


1 2
n−1
θy − (k + 1).
= θx + θy + −1 +
y
y

(9)

The differential operator l3 contains only derivatives in the x-direction. Similarly, the operator l6 contains only derivatives in the y-direction. Therefore, we
have the following theorem.
Theorem 1. The function Hnk (x, y) of two variables is holonomic.
Proof. Hnk (x, y) satisfies the following system of ODE’s
l3 • u(x, y)

=

0,

l6 • u(x, y)

=

0.

5

The next task is defining a Pfaffian system for the obtained (homogeneous)
ODE. From the Gröbner bases above, we have a Pfaffian system
∂
F = P2 F ,
∂y

∂
F = P1 F ,
∂x

(10)

where F (x, y) = t (1, θy , θy2 , θy3 ) • u(x, y) and


a1
1 
 0
P1 (x, y) =
xy  a5
a8
a1

3.2

=


a2 −1
0
a3 a2 + 1 −1 
,
a6
a7
n−1 
a9 a10
a11



0
1
0
P2 (x, y) = 
y 0
−b0

1
0
0
−b1

0
1
0
−b2


0
0 

1 
−b3

(k + 1)y

a2

= y−n+1

a3

=

(k + 1)y + n − 1

a5

=

(k + 1)xy 2

a6

= xy 2 − (n − 1)(xy + (k + 1)y + n − 1)

a7

= −(x + n − 1)y + (n − 1)(n − 2)

a8

= −(n − 2)(k + 1)xy 2

a9

=

(k − n + 3)xy 2 + (n − 1)2 (xy + (k + 1)y + n − 1)

a10

= xy 2 + (n − 1)2 (y − n + 2)

a11

= −xy − (n − 1)2

b3

= y + 2n − 4

b2

= n2 − 5n + 5 − (x + k + n)y

b1
b0

= xy 2 − (n − 1)xy − (k + 1)(n − 1)y − (n − 2)(n − 1)
= (ky + 1)xy.

Numerical evaluation

Since we have obtained the Pfaffian system for Hnk (x, y), we want to solve the
initial value problem in y-direction numerically by the Runge-Kutta method.
We try to implement the computation prototype by using Risa/Asir program to
evaluate the function at several points. The program includes these procedures:
• constructing an ODE

∂
∂y F

= P2 F for a vector-valued function F ,

• finding a initial value of F (x, y) by using a truncated series expression,
• applying the Runge-Kutta method to the ODE, advancing the vector valued function from the initial point y0 to endpoint ye .
To measure the accuracy of the evaluation by using the HGM, ideally, we
want to compare the result with some value which is equal or very close to the
exact value. In this case, we chose to use an approximated value of the function
Hnk (x, ye ) obtained by its truncated power series representation.

6

It is easy to prove that the series (5) converges for every point (x, y) in the
complex plane. Therefore, we can approximate Hnk (x, y) by its truncated series
fN (x, y) =

N
xk+1 X
(−1)q (k + 1)p+q
xp+q y p
k + 1 p,q=0 (n)p (1)p (1)q (k + 2)p+q

(11)

using a sufficiently large N .
When a point (x, y) and ε > 0 are given, we will determine an integer N
such that |fN +1 (x, y) − fN (x, y)| < ε |fN (x, y)|.
Table 1 show the computational result of Hnk (x, y) by using HGM and power
series with the parameters k = 2, n = 3, x = 1, and initial point y0 = 1.0 × 10−5
over several values of ye . In our experiment, we set ε = 1.0 × 10−6 and find out
that N = 10000 is sufficiently large enough for any of our test cases.
Table 1: Numerical evaluation for stability test
ye
16
1
0.1
0.01
0.001
0.0001
0.00002

Hnk (x, ye ) (series)

Hnk (x, ye ) (HGM)

3.21694
0.202444
0.164437
0.160982
0.160640
0.160606
0.160603

2.46313
0.377089
0.177692
0.162273
0.160769
0.160619
0.160606

We can see from the table that the results by the HGM agree with the series
expansion up to the sixth digit when the final point ye = 0.00002 = 2.0 × 10−5
is very close to the origin and the accuracy degrades really fast as ye becomes
larger. We conclude that the application of the HGM by the 4th order ODE
works only in a short range, that is, it can exterpolate the value at y = y0 within
its small neighborhood.
From the general theory of linear differential equations, the holonomic rank
of the system gives the dimension of the solution space. Since the solution space
of (10) is 4-dimensional, there exists a linearly independent solution F1 on the
corresponding solution to Hnk (x, y). The numerical instability happens because
the first entry of F1 is growing faster than the function Hnk (x, y) does. Consequently, as the value y becomes larger, this function becomes more dominant
than Hnk (x, y) and the evaluation of Hnk (x, y) by using the HGM gives very
significant error compared to its exact value.

4

Global Numerical evaluation of Hnk (x, y) by a
modified HGM

We have seen that a naive application of the HGM works well only locally. We
want to give an HGM for our problem which works in a long range evaluation.
In other words, we want to give a modified HGM which utilize one or more
inhomogeneous differential equations instead of the homogeneous ones.
7

4.1

A differential equation in x-direction

We will utilize the following inhomogeneous differential equation:
∂ k
H (x, y)
∂x n

= xk
= xk

∞ X
∞
X
xp+q y p (−1)q (k + p + q + 1)
p=0 q=0
∞
q
X
q=0

k −x

= x e

(n)p (1)p (1)q (k + p + q + 1)

∞
x (−1)q X (xy)p
(1)q p=0 (n)p (1)p

0 F1 (; n; xy).

(12)

The function 0 F1 (; n; z) itself satisfies the following differential equation:
{θz (θz + n − 1) − z} • 0 F1 (; n; z) = 0.

(13)

Substituting z = xy, the equation (13) becomes
{θx (θx + n − 1) − xy} • 0 F1 (; n; xy) = 0.

(14)

Using equations (12) and (13), we can construct an ODE in the x-direction
∂
G = QG,
∂x
where

Hnk (x, y)
G(x, y) =  0 F1 (; n; xy)  ,
θx • 0 F1 (; n; xy)




0
Q(x, y) =  0
0

xk e−x
0
xy


0
1/x  .
n−1

The following proposition is well-known:
Proposition 1. (well-known, e.g. [1]) Suppose that n > 0. There exists a
non-zero constant c such that
0 F1 (; n; z)

∼ c exp(2z 1/2 ) + o(exp(z 1/2 ))

holds when z → +∞.
Note that the function 0 F1 (; n; z) grows by the factor of exp(2z 1/2 ) as z
approaches infinity. In our numerical analysis, we want to avoid the exponential
growth. For that
reason, we need to apply a suitable Gauge transformation
√
to factor out e2 xy from the function 0 F1 (; n; xy) and find a new ODE that
corresponds to the new function.
√
Let f (x, √
y) = 0 F1 (; n; xy). Based on Proposition 1, we consider e−2 xy f (x).
√
Letting t = x and z = y, we define a new function
u(t) = e−2tz f (t2 ).

(15)

The annihilating differential equation for u(t) is
(θt2 + (4tz + 2(n − 1)θt + 2tz(2n − 1)) • u(t) = 0.

(16)

Also, from (12) and (15), we know that
2
∂ k 2 2
Hn (t , z ) = 2t2k+1 e−t +2tz u(t).
∂t

(17)

Using the system of differential equations formed by (16) and (17), we can
construct a new ODE.
8

Theorem 2. The vector-valued function H(t, z) = t (Hnk (t2 , z 2 ), u(t), θt • u(t))
satisfies an ODE
∂
H = RH,
∂t
where


−t2 +2tz 2(k+1)
t
0
1  0 2e
.
R(t, z) =
0
0
1
t
0
−2tz(2n − 1)
−4tz + 2(n − 1)
Note that, unlike the one that we obtained in Section 3, the vector-valued
function H is 3-dimensional. In Section 5, we will use the ODE above to evaluate
the values of Hnk (x, y) at several points y = λj and see that this ODE works
more stably than the previous one.

5

Numerical experiment

In this section, we use the ODE from Theorem 2 to evaluate numerically some
values of Hnk (x, y) at several points of (x, y) by using the HGM and calculate the
resulting outage probability. We also perform evaluation on the equation (2) by
using numerical integration and compare the result with the one we obtain by
using the HGM.
√
√
Let t = x and z = y. Assume that z is fixed. Consider the sequence
t0 , t1 , t2 , . . . , tn where t0 < t1 < t2 < . . . < tn = t and t0 is very small. First, we
compute the initial vector-valued function H(t0 , z) = Hnk (t20 , z 2 ), u(t), θt • u(t0 )
at the initial point t0 by using truncated series whose the method was explained in Section 3.2. Using the value H(t0 , z), we want to approximate the
value of H(t, z) by the Runge-Kutta method, utilizing the ODE that we obtained in Theorem 2. We repeat the same procedure to evaluate the values of
H(t1 , z), . . . , H(tn , z). Note that when evaluating H(tk , z), we use the value of
H(tk−1 , z) as the new initial value.

5.1

Setting

We use the GNU Scientific Library [2] to write a program in the C language that
performs computation by using the HGM and use the functions “NIntegrate”
and “Hypergeometric0F” of Wolfram Mathematica [13] to perform numerical
integration. More precisely, we evaluate Hnk (x, y) by
NIntegrate[t^k*Exp[-t]*Hypergeometric0F1[n,t*y],{t,0,x}];
Both computations are run on Raspberry Pi 3 Model B with the specification
according to Table 2. We choose this specific system as our computing platform
because of the free access to Wolfram Mathematica.
Table 2: System specification for numerical evaluation
System component

Specification

Operating system
CPU
Memory

Raspbian Linux
ARM Cortex-A53, 1.2GHz, 4 cores
1GB

9

We compute the values of Pr (x) at several points x and several sets of
eigenvalues. The number of eigenvalues represents the number of transmitter
antennas of a MIMO system. In computation, we use discrete points x ∈ {10z |
z = 1.0, 1.1, . . . , 2.0} the values of x are represented in the logarithmic unit dB.

5.2

Computations with different number of receiver antennas

This subsection gives the computational result with different numbers of receiver antennas. We want to see how the HGM compares to numerical integration in term of accuracy and computation time and observe how the number
of receiver antennas affects the evaluation. We present here the cases where
(NT , NR ) = (5, 5), (5, 6), (5, 7), (5, 8), (5, 9) because we want to see how well our
method works with a small number of antennas. Every case in this section use
the set of eigenvalues λ = {0.1, 0.2, 0.3, 0.4, 0.5} and step size h = 0.0001.

Figure 1: Plot of Pr (x) with various number of (NT , NR )

Table 3: Elapsed time for evaluation using different number of receiver antennas
(NT , NR )

Integration

HGM

(5,5)
(5,6)
(5,7)
(5,8)
(5,9)

18.97s
20.19s
20.45s
20.21s
20.25s

1.588538s
1.744283s
1.753809s
1.764655s
1.762677s

From Table 3, we can see that to achieve significantly accurate computational results, the HGM requires a much shorter time compared to numerical
integration. On the other hand, increasing the number of antennas seemingly
also increases the computation time for both methods although not significantly.
Note that, because HGM uses an iterative method to evaluate Pr (x), the table
does not clearly show the true advantage of this method over numerical inte-

10

gration.
To better understand this particular aspect, consider the case where (NT , NR )
= (5, 5). In this case, the sample x is a point of the set {10z | z = 1.0+0.1n, 0 ≤
n ≤ 10} of 11 elements. Suppose that we want to increase the number of samples by adding more middle points x at which we want to evaluate the function
Pr (x), like as {10z | z = 1.0 + 0.05n, 0 ≤ n ≤ 20} of 21 elements. In this case,
the HGM will also give an elapsed time around 1.6 seconds, not very far from
the time shown in Table 3. However, the numerical integration gives an elapsed
time of around 38 seconds, twice the computation time of the previous test.
This phenomenon happens because, when using HGM, we actually evaluate
the function Pr (x) at many points in the interval 10−5 ≤ x ≤ 102 . Therefore,
increasing the number of samples, as long as it is within the iteration and the
final point is fixed, will not affect the elapsed time significantly. On the other
hand, the numerical integration requires us to evaluate the function as many
times as the number of samples. Thus, if we increase the number of sample
points, the elapsed time for the numerical integration will increase, too.
Table 4: Computational result for the case (NT , NR ) = (5, 5)
x in dB

Pr (x) (Integration)

Pr (x) (HGM)

2
1.9
1.8
1.7
1.6
1.5
1.4
1.3
1.2
1.1
1

1
1
1
1
1
0.999952
0.997096
0.953502
0.743021
0.369936
0.100486

1
1
1
1
1
0.999952
0.997101
0.953567
0.743085
0.369961
0.100884

11

Table 5: Computational result for the case (NT , NR ) = (5, 6)
x in dB

Pr (x) (Integration)

Pr (x) (HGM)

2
1.9
1.8
1.7
1.6
1.5
1.4
1.3
1.2
1.1
1

1
1
1
1
0.999999
0.999812
0.991243
0.894786
0.569069
0.192951
0.0305385

1
1
1
1
0.999999
0.999813
0.991256
0.894912
0.569151
0.19297
0.0307046

Table 6: Computational result for the case (NT , NR ) = (5, 7)
x in dB

Pr (x) (Integration)

Pr (x) (HGM)

2
1.9
1.8
1.7
1.6
1.5
1.4
1.3
1.2
1.1
1

1
1
1
1
0.999997
0.999359
0.97728
0.797475
0.383415
0.0819006
0.00705474

1
1
1
1
0.999997
0.999362
0.977313
0.79768
0.3834992
0.081912
0.00710467

12

Table 7: Computational result for the case (NT , NR ) = (5, 8)
x in dB

Pr (x) (Integration)

Pr (x) (HGM)

2
1.9
1.8
1.7
1.6
1.5
1.4
1.3
1.2
1.1
1

1
1
1
1
0.999989
0.998067
0.948626
0.662382
0.223469
0.0280927
0.00124276

1
1
1
1
0.999987
0.998068
0.948646
0.662405
0.223537
0.0280977
0.00124381

Table 8: Computational result for the case (NT , NR ) = (5, 9)
x in dB

Pr (x) (Integration)

Pr (x) (HGM)

2
1.9
1.8
1.7
1.6
1.5
1.4
1.3
1.2
1.1
1

0.999984
1.00001
0.999994
1.00001
0.999961
0.994812
0.897575
0.50403
0.111381
0.00776535
0.000167867

1
1
1
1
0.999956
0.99482
0.897611
0.504057
0.111427
0.00776711
0.000168029

13

5.3

Computations with different set of eigenvalues λ

This subsection gives the computational result with several sets of eigenvalues.
We consider the MIMO system of 10 × 10 antennas. The first case uses the
eigenvalues λ = {1, 2, . . . , 10}, while the second case uses λ = {10, 20, . . . , 100}.
We also show the elapsed time of both computations.
Table 9: Computational result for the case λ = {1, 2, . . . , 10}
x in dB

Pr (x) (Integration)

Pr (x) (HGM)

2
1.9
1.8
1.7
1.6
1.5
1.4
1.3
1.2
1.1
1

1.00001
1.00001
0.99524
0.781498
0.148478
0.00203227
1.36179e-06
5.18699e-11
1.58216e-16
3.6096e-30
9.52818e-45

0.999998
0.999994
0.995231
0.781594
0.14857
0.0020352
1.36963e-06
5.21756e-11
1.56884e-16
5.64982e-23
3.77131e-30

Table 10: Computational result for the case λ = {10, 20, . . . , 100}
x in dB

Pr (x) (Integration)

Pr (x) (HGM)

2
1.9
1.8
1.7
1.6
1.5
1.4
1.3
1.2
1.1
1

4.6899e-07
1.28709e-15
1.02877e-27
9.98029e-43
6.28334e-60
1.42022e-78
4.01919e-98
5.1249e-118
7.29669e-138
2.20095e-157
2.43162e-176

4.70096e-07
1.3078e-15
1.05985e-27
1.03964e-42
6.66723e-60
1.44469e-78
4.07015e-98
5.2212e-118
7.29974e-138
2.25887e-157
2.73348e-176

We can see from Table 11 that the efficiency of the HGM, when compared
to the numerical integration, is more apparent when we increase the number
of antennas. This further observation emphasizes that the HGM has potential
to become a reliable numerical tool to evaluate the performance of MIMO systems, considering the trend that increasing the number of antennas improves
the quality of transmission.
14

Table 11: Computation time for different set of eigenvalues

6

Eigenvalues

Integration

HGM

{1,2,. . . ,10}
{10,20,. . . ,100}

140s
132.51s

6.666501s
6.615950s

Conclusions

In this paper, we treated the outage probability of a MIMO system in terms of
the CDF of the largest eigenvalue of non-central Wishart matrices. The CDF
involves the determinant of a certain matrix whose entries are integrals containing a hypergeometric function. We then derived the infinite series expression
of those entries from their integral representation. Using the resulting formula,
Hnk (x, y), we tried to apply an HGM to numerically evaluate its value at several
points of x and y’s. Our first try was using a 4th order homogeneous ODE in y
direction which annihilates Hnk (x, y) to construct Pfaffian system out of it. We
also compared our numerical computation result with the values computed by
using infinite series to briefly check how accurate our numerical values are. We
found that the HGM in the y-direction works in a short range because of the
instability of the ODE’s. Our alternative strategy is to use a modified HGM,
utilizing an inhomogeneous ODE for Hnk (x, y) instead of the homogeneous one.
This time, we compared the numerical values and the computation time with the
ones obtained by using the numerical integration on Mathematica. As the result,
the computation gives a very accurate result within a much shorter time than
numerical integration. We also explained that this advantage becomes more
apparent as the number of samples increases. Based on these observations, we
conclude that the modified HGM performs better to compute the outage probability of a given MIMO system in terms of accuracy and computation time than
the numerical integration.

References
[1] A. Dixit, V. H. Moll, The integrals in Gradshteyn and Ryzhik, Part 28:
The confluent hypergeometric function and Whittaker functions. Scientia,
Series A: Mathematical Sciences 26 (2015), 49–62.
[2] GSL — GNU Scientific Library, https://www.gnu.org/software/gsl/,
1996–2017.
[3] H. Hashiguchi, Y. Numata, N. Takayama, A. Takemura, Holonomic gradient method for the distribution function of the largest root of a Wishart
matrix, Journal of Multivariate Analysis 117 (2013), 296–312.
[4] T. Hibi (eds), Gröbner Bases — Statistics and Software Systems, Springer,
2013.
[5] M. Kang, M. Alouini, Largest Eigenvalue of Complex Wishart Matrices and
Performance Analysis of MIMO MRC Systems, IEEE Journal on Selected
Areas in Communications 21 (2003), 418–426.
15

[6] T. Kimura, Hypergeometric functions of two variables, Seminar Note in
Math., Univ. of Tokyo, 1973.
[7] T. Koyama, H. Nakayama, K. Nishiyama, N. Takayama, Holonomic Gradient Descent for the Fisher-Bingham Distribution on the d-dimensional
Sphere, Computational Statistics 29 (2014), 661-683
[8] G. Letac, H. Massam, A tutorial on non central Wishart distributions,
http://www.math.univ-toulouse.fr/~letac/Wishartnoncentrales.
pdf, 2004.
[9] W. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical
Recipes, third edition, Cambridge University Press, 2007.
[10] Risa/Asir developing team, Risa/Asir, a computer algebra system, http:
//www.math.kobe-u.ac.jp/Asir/, 1994–2017.
[11] T. Sei, A. Kume, Calculating the Normalising Constant of the Bingham
Distribution on the Sphere using the Holonomic Gradient Method, Statistics and Computing 25 (2015), 321-332
[12] C. Siriteanu, A. Takemura, S. Kuriki, H. Shin, C. Koutschan, MIMO ZeroForcing Performance Evaluation Using the Holonomic Gradient Method,
IEEE Transactions on Wireless Communications 14 (2015), no. 4, 2322–
2335.
[13] Wolfram Research, Inc., Mathematica, 1993–2017.
[14] R. Vidunas, A. Takemura, Differential relations for the largest root distribution of complex non-central Wishart matrices, arXiv:1609.01799
[15] References for HGM, http://www.math.kobe-u.ac.jp/OpenXM/Math/
hgm/ref-hgm.html

16

