Holonomic Gradient Method-Based CDF Evaluation for the Largest Eigenvalue of a Complex Noncentral Wishart Matrix

# Holonomic Gradient Method-Based CDF Evaluation for the Largest Eigenvalue of a Complex Noncentral Wishart Matrix

Fadil H. Danufane, Katsuyoshi Ohara, Nobuki Takayama, Constantin Siriteanu F. H. Danufane was with Kanazawa University, Japan. He is now with the Research Center for Electronics and Telecommunications, Indonesian Institute of Sciences, Bandung, Indonesia.K. Ohara is with the Faculty of Mathematics and Physics, Kanazawa University, JapanN. Takayama is with the Department of Mathematics, Kobe University, Japan.C. Siriteanu is with the Graduate School of Information Science and Technology, Osaka University, Japan.
###### Abstract

The outage probability of maximal-ratio combining (MRC) for a multiple-input multiple-output (MIMO) wireless communications system under Rician fading is given by the cumulative distribution function (CDF) for the largest eigenvalue of a complex noncentral Wishart matrix. This CDF has previously been expressed as a determinant whose elements are integrals of a confluent hypergeometric function. For the determinant elements, conventional evaluation approaches, e.g., truncation of infinite series ensuing from the hypergeometric function or numerical integration, can be unreliable and slow even for moderate antenna numbers and Rician -factor values. Therefore, herein, we derive by hand and by computer algebra also differential equations that are then solved from initial conditions computed by conventional approaches. This is the holonomic gradient method (HGM). Previous HGM-based evaluations of MIMO relied on differential equations that were not theoretically guaranteed to converge, and, thus, yielded reliable results only for few antennas or moderate . Herein, we reveal that gauge transformations can yield differential equations that are stabile, i.e., guarantee HGM convergence. The ensuing HGM-based CDF evaluation is demonstrated reliable, accurate, and expeditious in computing the MRC outage probability even for very large antenna numbers and values of .

## I Introduction

### I-a Background and Model

For multiple-input multiple-output (MIMO) wireless communications systems, maximal-ratio combining (MRC) is often considered for analysis and implementation[1][2]. In MIMO MRC, each data symbol is transmitted into and received from the dominant singular mode of the channel matrix . This relatively low-complexity technique111Only knowledge of the dominant right and left singular vectors of is needed, at the transmitter and receiver, respectively. maximizes the symbol-detection signal-to-noise ratio (SNR) and, thus, the diversity gain.

The MRC SNR is proportional to the largest eigenvalue of [1, Eq. (27)]. Therefore, the MRC outage probability, i.e., the probability of failing to achieve an SNR that ensures reliable symbol detection[1][2], is given by the cumulative distribution function (CDF) of .

For complex circularly-symmetric Gaussian-distributed with nonzero mean , i.e., for Rician fading, and zero row correlation, matrix has the complex noncentral Wishart distribution[1][3]. Then, the distribution of has been characterized several times before, for various settings, with expressions that are difficult to evaluate reliably, as discussed next.

### I-B Previous Work and Its Limitations

For the largest eigenvalue of a central Wishart matrix, Hashiguchi et al.[4] computed the CDF expressed as a confluent hypergeometric function with matrix third argument[4, Eq. (4)]. Because here we are interested in the complex-valued noncentral Wishart distribution, we shall start from Kang and Alouini’s paper[1]. For the case when mean has arbitrary rank but zero row and column correlation, their[1, Eq. (2)] is a determinantal expression for the CDF of the largest eigenvalue of , where is the Rician factor. The CDF determinant elements are integrals of the special-case confluent hypergeometric function , where stands for the eigenvalues of , and is the integration variable. However, because for the eigenvalues of the typical normalization yields , larger , , or implies the requirement to evaluate for larger .

The obvious approach to compute the elements of the CDF determinant is numerical integration of native implementations of in software tools such as Mathematica or Matlab222Technical details of such implementations are not publicly available.. Another approach is truncation of the infinite series ensuing from the well-known formula[5, Eq. (13.2.2), p. 322] for . However, as becomes larger, these conventional approaches can lead to numerical inaccuracy (and even divergence) and/or long computation time, as illustrated in this paper. Note that[1] showed MRC outage probability results only for small antenna numbers, i.e., , and small Rician factor, i.e.,  dB.

For with arbitrary-rank and uncorrelated columns and rows, Wu et al. have recently expressed the CDF of the dominant eigenvalue of in[2, Eq. (3)]. However, their expression is an infinite series whose coefficients , , expressed in[2, Eq. (4)], are difficult to express beyond . Further, the series was approximated with only its first term in[2, Eq. (19)], and MRC outage probability results were shown only for small antenna numbers, i.e., , and small-to-moderate Rician factor, i.e.,  dB.

Furthermore, these MIMO performance evaluation limitations are not MRC-specific. As another example, for MIMO multiplexing with zero-forcing (ZF) detection, which requires knowledge of only at the receiver, Siriteanu et al. analyzed and evaluated the performance for Rician fading with for a special case in[6] and the general case in[7]. The moment generating function of the ZF symbol-detection SNR, which for stream is proportional to , was characterized in[6, Eq. (19)][7, Eq. (52)] with the confluent hypergeometric function of third argument proportional with [7, Eq. (66)]. The truncation of ensuing infinite series that characterize the ZF performance was shown to diverge numerically333This is because the infinite series for diverges numerically for large values of the third argument[8]. for , , and moderate values in[6][7].

Such limitations are to be expected when evaluating any MIMO transceiver technique for Rician fading because the distribution of matrix , which determines the performance, is then noncentral Wishart, i.e., characterized by a confluent hypergeometric function with matrix argument[3, Eq. (3)] that is inherently difficult to manipulate and compute[9].

### I-C Alternative Computation by the Holonomic Gradient Method (HGM)

As mentioned above, Siriteanu et al.[6][7] evaluated MIMO ZF performance for Rician fading with . First, they derived, by hand in[6] and by computer algebra in[7], linear ordinary differential equations (LODEs) satisfied by infinite series that describe ZF performance measures, e.g., SNR probability density function in[6] and outage probability in[7]. Then, they numerically evaluated these measures by a hybrid approach that solves LODEs for a desired value of by starting from initial conditions computed by truncating infinite series for a value of that is sufficiently small to ensure numerical convergence for the truncation. This approach has become known as the holonomic gradient method (HGM)[6][7]. However, HGM was applied in[4] to compute the largest-eigenvalue CDF only for a central Wishart matrix (i.e., for ). On the other hand, the HGM-based evaluations of MIMO ZF for Rician fading from[6][7] were found reliable only for small MIMO deployments (e.g., , ) or for small-to-moderate values of ( dB).

HGM introductions, applications, and references appear in[6][7][10, Ch. 6][11]. HGM entails concepts, methods, and terminology from differential equations, abstract algebra, and Gröbner basis (or generator) computation, e.g., holonomic function and system, Pfaffian system, ring, ideal, etc., which are used herein but are thoroughly introduced in[10, Ch. 1][12][13]. Gröbner bases applications in signal processing are presented in[14].

### I-D Contribution: Reliable HGM-Based MIMO MRC Evaluation

Herein, for the evaluation of the CDF of , i.e., the MRC outage probability, for uncorrelated Rician fading with of arbitrary rank, we propose a reliable HGM-based approach that avoids the limitations of numerical integration and series truncation, as well as limitations of the HGM-based MIMO ZF evaluations from[6][7]. Our contributions are as follows:

1. From Kang and Alouini’s integral expression of the CDF determinant elements[1, Eq. (2)], using the infinite series known for , we derive an infinite series and then deduce LODEs it satisfies from Gröbner bases obtained by computer algebra. Thereafter, we find that the ensuing HGM is reliable only for small , as for ZF in[6][7]. We reveal the cause of this limitation to be that the used LODEs are not stabile444Concept defined later..

2. From the integral expression, using the differential equation known for , we derive a new LODE system and prove it stabile, thus theoretically ensuring HGM effectiveness. This guarantees reliable HGM-based evaluation of MRC under Rician fading.

3. We reveal that the Step-2 stabile LODEs can also be obtained by a gauge transformation[15] from the Step-1 LODEs, and that stabile LODEs can be obtained algorithmically by gauge transformations of other LODEs, in general.

4. We illustrate numerically that, unlike conventional integration and series truncation, HGM based on the Step-2 stabile LODEs, and further enhanced by gauge transformations, yields reliable evaluation of the CDF of even for very large values of , i.e., of the MRC outage probability for very large , , and Rician -factor. HGM is also much faster than integration and Monte-Carlo simulation.

### I-E Notation

• Scalars, vectors, and matrices are represented with lowercase italics, lowercase boldface, and uppercase boldface, respectively, e.g., , , and ; superscripts and stand for transpose and Hermitian (i.e., complex-conjugate) transpose; is the identity matrix.

• stands for the enumeration .

• stands for ‘proportional to’; stands for ‘approximately equal’.

• denotes an complex-valued circularly-symmetric Gaussian vector with mean and covariance matrix ; an complex-valued circularly-symmetric Gaussian random matrix with mean , row covariance , and column covariance , i.e., a matrix whose vectorized form is distributed as , is denoted herein as , based on the definition from[7, Sec. I.F]; subscripts and identify, respectively, deterministic and random components; subscript indicates a normalized variable; denotes statistical average;

• stands for the th partial derivative with respect to variable of function ; is the theta differential operator w.r.t. ; represents the application of differential operator to function .

• is the confluent hypergeometric function[5, Ch. 13]; is the Pochhammer symbol, i.e., and , , and .

### I-F Paper Organization

Section II introduces the MIMO signal and channel model as well as the MRC method and its SNR. Section III gives the determinantal expression derived by Kang and Alouini in[1] for the CDF of , i.e., the MRC outage probability. For the determinant elements, which are integrals of , we derive an infinite series. Section IV derives LODEs satisfied by this infinite series, by hand and by computer algebra, and then implements HGM. Section V derives by hand alternative LODEs from the integral expression of the determinant elements, and proves them stabile. For an enhanced HGM based on these LODEs and gauge transformations, Section VI shows results for various antenna numbers, sets of eigenvalues of , and values of .

## Ii MIMO Model, MRC SNR, Outage Probability

We consider a point-to-point MIMO wireless communications link. The transmitter is equipped with an -element antenna array that transmits the vector , where is the complex-valued transmitted symbol taken from a constellation with unit average energy, is the energy consumed per transmitting antenna, i.e., is the energy consumed per transmitted symbol, and is the unit-norm transmit combiner (beamformer).

Herein, we assume that the transmitted signal vector encounters a fading radio channel characterized by the matrix distributed as , with transmit-side correlation matrix given by . Thus, we assume Rician fading that is uncorrelated both at the transmitter and the receiver. With its deterministic and random components denoted as and , respectively, we can write the channel matrix as

 H=Hd+Hr=√KK+1Hd,n+√1K+1Hr,n, (1)

where and are normalized as usual, i.e.,

 ∥Hd,n∥2=E{∥Hr,n∥2}=NRNT,which implies E{∥H∥2}=NRNT. (2)

Above, is the Rician factor, i.e.,

 K=∥Hd∥2E{∥Hr∥2}=KK+1∥Hd,n% ∥21K+1E{∥Hr,n∥2}. (3)

 r=√EbNTHwTb+n. (4)

We define the per-symbol transmit SNR as

 Γb=EbN01N% T. (5)
###### Remark 1.

For the MIMO signal model in (4), MRC maximizes the symbol-detection SNR by beamforming with the dominant right and left singular vectors of at the transmitter and receiver, respectively. Consequently, the MRC SNR is given by multiplied with the maximum eigenvalue of matrix — see[1, Sec. III.B] for the simple proof.

In general, the outage probability is defined as the probability that the symbol-detection error probability — instantaneous, i.e., given — exceeds a certain threshold, i.e., . The threshold SNR that yields the threshold error probability is denoted herein with . For instance, for QPSK modulation, a symbol error probability level of can be achieved for  dB. Thus, the outage probability of a detection method can also be defined as the probability that its symbol-detection SNR does not exceed the threshold value . Clearly, the outage probability depends on the distribution of symbol-detection SNR, which is determined by the distribution of the MIMO channel matrix .

In the next section, we shall characterize the MRC outage probability based on Kang and Alouini’s work from[1]. They considered the matrix , with eigenvalues555All these eigenvalues are real-valued and positive because, with probability 1, matrix is Hermitian and has rank . , and the corresponding matrix with eigenvalues , where and . They characterized the CDF of the dominant eigenvalue of based on the eigenvalues of . Nevertheless, they later assumed zero transmit-side correlation in order to express the MRC SNR and its outage probability.

###### Remark 2.

The MRC SNR characterized in Remark 1, is given for by the dominant eigenvalue of as follows

 ρMRC=ΓbK+1ϕs, (6)

and the MRC outage probability is

 Pr(ρMRC≤Γth)=Pr(ϕs≤(K+1)ΓthΓb). (7)

Expression (7) confirms that a larger per-symbol transmit SNR (i.e., a larger transmitted energy) yields a lower outage probability whereas a larger threshold SNR (i.e., a higher quality of service) yields a higher outage probability. It also confirms that the MRC outage probability depends on the distribution of , which is characterized next, based on[1].

## Iii CDF of Largest Eigenvalue of Complex Noncentral Wishart Matrix

We begin with results from[1] that apply for the distribution . Then, for analysis tractability, we specialize to the case with .

###### Definition 1.

If , then the matrix has a complex noncentral Wishart distribution with the following parameters: dimension , degrees of freedom , and noncentrality matrix [1][3][7].

Now we state the result proved by Kang and Alouini in[1].

###### Theorem 1.

If and the matrix has nonzero distinct eigenvalues , then the CDF of the largest eigenvalue of evaluated at certain threshold value is given by the determinantal expression

 Pr(ϕs≤x)=e−(λ1+⋯+λs)[∏1≤i

where is an matrix whose th element is given by

 [Φ(x)]i,j=∫x0yt−ie−y0F1(;t−s+1;λjy)dy. (9)

Recall from Remark 1 that the MRC SNR is given by the dominant eigenvalue of , regardless of the column-correlation matrix . However, because the dependence between the eigenvalues of and is not straightforward for , herein, we employ Theorem 1 for MRC analysis only for666As in[1], where Kang and Alouini first derived the CDF for the dominant eigenvalue of in[1, Eq. (2)], and then assumed zero correlation between the columns of in order to express the MRC outage probability in[1, Eq. (29)]. . Then, Theorem 1 characterizes the CDF of the dominant eigenvalue which enters the MRC SNR as in (6), and determines the MRC outage probability as in (7).

###### Remark 3.

Our model and assumptions yield , i.e., the eigenvalues , , of satisfy . Consequently, even a relatively small MIMO deployment with , in moderate Rician fading with  dB (i.e., ) yields . Much larger values for ’s ensue with increasing , and . For instance, a truly massive MIMO system with in strong Rician fading with  dB, yields .

###### Objective.

Given the eigenvalues , , and a threshold value , we seek the ability to compute the value of the CDF by first evaluating the function in (9) and then the determinant in (8). Finally, scaling as in (7) yields the MRC outage probability.

Defining and , the determinant entry from (9) can be rewritten as

 Hkn(x,λ)=∫x0yke−y0F1(;n;λy)dy. (10)

One may then attempt to compute by numerical integration, e.g., in Mathematica. However, we show later that this approach is not time-efficient for realistically-large values of . Therefore, next, we derive an alternative expression and discuss its computation.

By employing for the confluent hypergeometric function its well-known infinite series [5, Eq. (13.2.2), p. 322] and then changing integration and summation order in (10) we have obtained the following infinite series

 Hkn(x,λ)=xk+1k+1∞∑p=0∞∑q=0(−1)q(k+1)p+q(n)p(1)p(1)q(k+2)p+qxp+qλp, (11)

which can be proven to converge . One may attempt to compute with truncation

 Hkn(x,λ)≈Hkn,N(x,λ)=xk+1k+1N∑p=0N∑q=0(−1)q(k+1)p+q(n)p(1)p(1)q(k+2)p+qxp+qλp, (12)

with determined as follows: given the pair and a positive scalar , we select as the smallest integer that yields numerical convergence, i.e., satisfies .

Thus, we have attempted to compute by implementing the series truncation (12) in C with double-precision arithmetic777The series truncation computations presented here have been executed on the machine described in Section VI on page VI., and compared for correctness and execution time against computations in Risa/Asir with rational arithmetic[16]. Thus, we have found that the double-arithmetic C implementation of series truncation (12), with , can compute accurately for taking increasing values in the interval ; numerical convergence within is then achieved for values that increase in the interval , but the computation takes negligible time (i.e., less than a microsecond).

However, for larger values of and (which may arise in practice, cf. Remark 3), the double-arithmetic C implementation of series truncation (12) can experience computational difficulties. For example, we have found that the computation of begins to break down for because denominators in the terms of series (12) reach the upper limit of double-precision representation. Consequently, for , these denominators are treated as ‘infinite’ and subsequent series terms remain constant, which precludes numerical convergence. The double-arithmetic C computation of series truncation (12) with takes about 45 microseconds (but does not achieve numerical convergence). On the other hand, the rational-arithmetic Risa/Asir implementation of series truncation (12) successfully achieves numerical convergence within for , but takes milliseconds.

These limitations of integration and truncation have motivated us to seek more reliable and efficient alternative approaches. Thus, we propose below new evaluation approaches based on differential equations and HGM. First, in Section IV, we derive a -dimensional system of LODEs for from the infinite series in (11) based on Gröbner basis computation, but find the ensuing HGM ineffective. Then, in Section V, we use the integral relationship (10) to derive a -dimensional system of LODEs for and propose a series of HGM enhancements that lead to accurate CDF evaluation.

## Iv HGM for Hkn(x,λ) Using Differential Equations for Infinite Series (11)

### Iv-a Holonomic Functions and the Holonomic Gradient Method (HGM)

HGM is a relatively new approach to numerically evaluate infinite series, such as our series for from (11). It relies on the notion of holonomic function and on the computation of Gröbner bases.

###### Definition 2.

A complex-valued function , where , is holonomic if and only if, , there exists a LODE of some order w.r.t. , as follows

 {pi,Ni(x)∂Ni∂xNii+⋯+pi,1(x)∂∂xi+pi,0(x)}∙f(x)=0, (13)

where are polynomials. The set of LODEs w.r.t. all variables is called the holonomic system satisfied by .

###### Definition 3.

A Pfaffian system is a set of LODEs of the form

 ∂∂xif(x)=Pi(x)f(x),i=1,…,n, (14)

where is a function vector, and matrices satisfy the compatibility (integrability) condition , [10, p. 296].

Pfaffian LODE systems such as (14) are desirable because they can be solved numerically, e.g., by Runge–Kutta method, given the initial condition value for some .

###### Definition 4.

HGM for the computation of a function entails the following steps[10, Ch. 6]:

1. Derive a holonomic system for , e.g., by hand or by computer algebra.

2. Recast the holonomic system as a Pfaffian system.

3. Determine an initial condition for some , e.g., by series truncation.

4. Numerically evaluate at desired point , e.g., by the Runge–Kutta method.

### Iv-B Differential Equations and Pfaffian System from Infinite Series (11) for Hkn(x,λ)

###### Theorem 2.

Function represented by the infinite series in (11) is holonomic because it satisfies the following LODEs w.r.t.  and , respectively,

 {[(θx+x−k−1)(θx+x−k+n−2)−xλ]θx}∙Hkn(x,λ)=0, (15) {θ4λ+(λ+2n−4)θ3λ+[n2−5n+5−(x+k+n)λ]θ2λ+[xλ2−(n−1)xλ −(k+1)(n−1)λ−(n−2)(n−1)]θλ+(kλ+1)xλ}∙Hkn(x,λ)=0, (16)

where and . These LODEs can be recast as the Pfaffian system

 ∂∂xf(x,λ) = P(x,λ)f(x,λ), (17) ∂∂λf(x,λ) = Q(x,λ)f(x,λ), (18)

where is the -dimensional vector

 f(x,λ)=(1θλθ2λθ3λ)T∙u(x,λ), (19)

with being a generic function888 stands for all four linearly independent solutions of the -dimensional system of LODEs in (17)., and

 P(x,λ)=1xλ⎛⎜ ⎜ ⎜ ⎜⎝a1  a2  −1  00  a3  a2+1  −1a5  a6  a7  n−1a8  a9  a10  a11⎞⎟ ⎟ ⎟ ⎟⎠,Q(x,λ)=1λ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝0 1 0 00 0 1 00 0 0 1−b0 −b1 −b2 −b3 ⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (20)

where and are polynomials in and shown in Table I.

###### Proof:

Shown in Appendix A; it entails 1) by-hand deduction of two partial differential equations, and 2) LODEs deduction by computer-algebra-based computations of Gröbner bases.      //

### Iv-C Performance of HGM w.r.t. λ Based on LODE System (18)

Recall from Remark 3 that large values for , and can yield large magnitudes for the eigenvalues , , of . Thus, we have attempted to evaluate the function over a range of values for with HGM w.r.t. , based on system (18). The numerical procedures and settings are summarized on the second line in Table II and relevant footnotes on page II.

Our HGM w.r.t.  based on system (18) for comprises the following steps:

• Compute at some initial value (i.e., the initial condition), by truncation (12).

• Solve the LODE system in (18), and thus advance from initial value to the value of interest , by the Runge–Kutta method.

We have attempted to compute by the above HGM procedure from the initial condition obtained by series truncation, and compared with results obtained (upon numerical convergence) from the double-arithmetic C implementation of series truncation (12). We have found that they agree closely when is near the initial value , but they begin to disagree as becomes larger, e.g., . We conclude that the application of HGM based on the -dimensional LODE system from (18) is reliable only locally, i.e., within a narrow range around small . The reason is explained further below.

## V HGM for Hkn(x,λ) Using Stabile LODEs

We have seen above that HGM based on (18) is reliable only locally. This phenomenon can be understood based on the notion, defined below, of a LODE (system) that is stabile999Which is different in meaning and spelling than the conventional notion of (Lyapunov) stable ODE. for its solution. In this section, we first directly derive a new LODE system stabile for from integration (10), and then propose gauge transformations as a general means to obtain from a given LODE system that is not stabile a lower-dimensional LODE system that is stabile.

### V-a LODE Stabile for Its Solution

###### Definition 5.

Consider a holonomic function101010Here, the generic variable symbol stands for either of the variables of . that satisfies a LODE of order , and denote with its linearly independent solutions. Then, let be the dominant solution for , i.e., , . We refer to the LODE as stabile for if . (Note that a LODE is stabile or not regardless of the selected set of linearly independent solutions.) The notion of stabile LODE is defined analogously in the case of a vector-valued function, by replacing with a vector norm .

For example, the solution space of LODE is spanned by , which is the dominant solution, and by . Consider the general LODE solution . Then, if the LODE is stabile for ; otherwise, the LODE is not stabile for .

###### Remark 4.

Typically, the initial condition employed in evaluating function by HGM is affected by some error (e.g., due to series truncation). Then, HGM using a stabile LODE can theoretically111111I.e., for computation on a machine with unlimited computational power and arithmetic precision. still yield the function accurately as increases, whereas HGM using a LODE that is not stabile may yield a solution of much larger magnitude than . Thus, the latter can be effective in evaluating only locally, for around some small .

###### Theorem 3.

The LODE (16) and ensuing -dimensional LODE system (18) are not stabile for when . On the other hand, LODE (15) is stabile whereas ensuing -dimensional LODE system (17) is not stabile when .

###### Proof:

See Appendix B-A.      //

###### Remark 5.

The fact that LODE (16) and system (18) are not stabile explains the ineffectiveness of the HGM w.r.t.  from Section IV.

### V-B Direct Derivation From Integration (10) of New LODE System Stabile for Hkn(x,λ)

By differentiating (10) w.r.t. , we obtain the following inhomogeneous equation for :

 ∂∂xHkn(x,λ)=xke−x0F1(;n;xλ). (21)

Now, for any integer , there exists constant such that[17]

 limz→∞|0F1(;n;z)−cne2√z(√z)−n+1/2|e2√z(√z)−n+1/2→0. (22)

Thus, for large , increases with as rapidly as . We can avoid the exponential growth by working instead with the function (in which we assume to be a given constant)

 v(x)=e−2√xλ0F1(;n;xλ)=exe−2√xλx−k−1θx∙Hkn(x,λ), (23)

where the second equality follows from and (21).

###### Remark 6.

Using the notation and , with constant, functions , , and become , , and , respectively. Because for , from (23) we can write .

The following theorem gives the LODE system w.r.t.  satisfied by and .

###### Theorem 4.

The vector valued function satisfies the following -dimensional LODE system:

 ∂∂φg(φ)=1φ⎛⎜⎝02e−φ2+2φψφ2(k+1)00010−2(2n−1)φψ−[4φψ+2(n−1)]⎞⎟⎠g(φ). (24)
###### Proof:

See Appendix B-B; it uses differential equation[5, Eq. (15.10.1)].      //

###### Theorem 5.

The -dimensional LODE system (24) is stabile for when , .

###### Proof:

Sketched in Appendix B-C.      //

Next, we reveal that the stabile -dimensional LODE system (24) can be deduced systematically from the -dimensional LODE system (17), which is not stabile, and that this is possible in general. Finally, Section VI demonstrates with numerical results that the HGM w.r.t.  based on (24) yields reliable computation even for very large values of .

### V-C Gauge Transformation-Based Deduction of Stabile LODEs

Above, we derived the stabile LODE (24) from the inhomogeneous differential equation for in (21) after accounting for the asymptotic behavior of . Nevertheless, as we state below, in general, a stabile LODE can be found by applying a suitable gauge transformation to an available LODE. Gauge transformations are widely used in mathematical physics[15]. Our definition for them here is analogous but simpler.

###### Definition 6.

Consider the LODE

 f′(x)=ddxf(x)=P(x)f(x), (25)

where is a matrix-valued function and is a vector-valued function. Let be an invertible matrix and be the vector that satisfies , so that and, hence,

 h′=(G−1PG−G−1G′)h. (26)

Then, is called the gauge transformation that yields LODE (26) from LODE (25).

The proofs of the following two theorems are not included due to manuscript length limitation.

###### Theorem 6.

A gauge transformation yields the stabile -dimensional LODE system (24) from the -dimensional LODE system (17), which is not stabile, by block upper triangulation.

###### Theorem 7.

From a given LODE system that is not stabile for a target function, a lower-dimensional stabile LODE system can be derived algorithmically by gauge transformations.

###### Remark 7.

Gauge transformations can help convert LODEs obtained by hand and computer algebra into stabile LODEs that then guarantee the numerical convergence of HGM-based performance measure evaluation, also for other MIMO transceiver methods and fading models.

### V-D Further Enhancements to the Proposed HGM

Recall from Remark 4 that if a LODE system is stabile for a function then it can be computed by HGM for arbitrary values of its variables, theoretically. Practically, we have found that HGM w.r.t.  based on the stabile LODE system (24) leads for large to exceeding the largest number representable in double arithmetic in C121212The largest double is , whereas numerical integration in Mathematica yields .. Nevertheless, we shall demonstrate that we can overcome this limitation by enhancing HGM with the following modifications:

1. Maintain values within bounds of the double representation by applying to the vector from Theorem 4 the following additional gauge transformations:

 G2 = ⎛⎜⎝exp(−φ2+2φψ)00010001⎞⎟⎠,forφ<ψ, (27)
 G3 = ⎛⎜⎝exp(ψ2)00010001⎞⎟⎠,forφ≥ψ. (28)

On the one hand, in matrix , the exponent takes its maximum value when and becomes when or . On the other hand, matrix is constant because is considered constant. While these transformations have been selected heuristically, they have been found reliable. For example, and yield , i.e., within the double representation range.

2. Use numerical integration of (10) to obtain initial conditions, as series truncation no longer converges in a reasonable time for large .

3. Use big float131313I.e., arbitrary precision. Note that, as required precision increases, computation time can also increase. in computing the determinant of the matrix with elements . Relevant experiments have been implemented in Risa/Asir.

4. Estimate the numerical errors arising in the evaluation of the determinant by adding random errors of the same size as the numerical error of to each element of .

Numerical results from HGM enhanced as above are shown in Sections VI-B and VI-C.

## Vi Numerical Results

Below, we describe our results from evaluating the CDF of the dominant eigenvalue of for several antenna number pairs and eigenvalue sets . Recall that this CDF yields the MRC outage probability, after accounting for the scaling in (7). We shall find that, unlike conventional methods, our HGM is consistently reliable and efficient for CDF computation, even in extreme conditions, i.e., very large ’s (or , , and ).

We employed a computer with the Debian “wheezy” operating system, a Xeon E5-4650 processor at 2.7 GHz, and 256 GB of memory. Table II summarizes the numerical procedures and settings employed for CDF computation, whereas footnotes therein describe the software packages used for the implementation of these procedures as well as further settings.