Holonomic Gradient MethodBased
CDF Evaluation for the Largest Eigenvalue of
a Complex Noncentral Wishart Matrix
Abstract
The outage probability of maximalratio combining (MRC) for a multipleinput multipleoutput (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 HGMbased 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 HGMbased 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
Ia Background and Model
For multipleinput multipleoutput (MIMO) wireless communications systems, maximalratio 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 lowcomplexity technique^{1}^{1}1Only knowledge of the dominant right and left singular vectors of is needed, at the transmitter and receiver, respectively. maximizes the symboldetection signaltonoise 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 circularlysymmetric Gaussiandistributed 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.
IB 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 complexvalued 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 specialcase 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 Matlab^{2}^{2}2Technical details of such implementations are not publicly available.. Another approach is truncation of the infinite series ensuing from the wellknown 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 arbitraryrank 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 smalltomoderate Rician factor, i.e., dB.
Furthermore, these MIMO performance evaluation limitations are not MRCspecific. As another example, for MIMO multiplexing with zeroforcing (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 symboldetection 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 numerically^{3}^{3}3This 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].
IC 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 largesteigenvalue CDF only for a central Wishart matrix (i.e., for ). On the other hand, the HGMbased evaluations of MIMO ZF for Rician fading from[6][7] were found reliable only for small MIMO deployments (e.g., , ) or for smalltomoderate 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].
ID Contribution: Reliable HGMBased 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 HGMbased approach that avoids the limitations of numerical integration and series truncation, as well as limitations of the HGMbased MIMO ZF evaluations from[6][7]. Our contributions are as follows:

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 stabile^{4}^{4}4Concept defined later..

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 HGMbased evaluation of MRC under Rician fading.

We reveal that the Step2 stabile LODEs can also be obtained by a gauge transformation[15] from the Step1 LODEs, and that stabile LODEs can be obtained algorithmically by gauge transformations of other LODEs, in general.

We illustrate numerically that, unlike conventional integration and series truncation, HGM based on the Step2 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 MonteCarlo simulation.
IE 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., complexconjugate) transpose; is the identity matrix.

stands for the enumeration .

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

denotes an complexvalued circularlysymmetric Gaussian vector with mean and covariance matrix ; an complexvalued circularlysymmetric 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 .
IF 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 pointtopoint MIMO wireless communications link. The transmitter is equipped with an element antenna array that transmits the vector , where is the complexvalued 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 unitnorm transmit combiner (beamformer).
Herein, we assume that the transmitted signal vector encounters a fading radio channel characterized by the matrix distributed as , with transmitside 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
(1) 
where and are normalized as usual, i.e.,
(2) 
Above, is the Rician factor, i.e.,
(3) 
Finally, the receiver frontend adds the complex additive noise vector . Thus, the dimensional received signal vector is
(4) 
We define the persymbol transmit SNR as
(5) 
Remark 1.
For the MIMO signal model in (4), MRC maximizes the symboldetection 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 symboldetection 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 symboldetection SNR does not exceed the threshold value . Clearly, the outage probability depends on the distribution of symboldetection 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 eigenvalues^{5}^{5}5All these eigenvalues are realvalued 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 transmitside 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
(6) 
and the MRC outage probability is
(7) 
Expression (7) confirms that a larger persymbol 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.
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
(8) 
where is an matrix whose th element is given by
(9) 
Recall from Remark 1 that the MRC SNR is given by the dominant eigenvalue of , regardless of the columncorrelation matrix . However, because the dependence between the eigenvalues of and is not straightforward for , herein, we employ Theorem 1 for MRC analysis only for^{6}^{6}6As 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.
Defining and , the determinant entry from (9) can be rewritten as
(10) 
One may then attempt to compute by numerical integration, e.g., in Mathematica. However, we show later that this approach is not timeefficient for realisticallylarge values of . Therefore, next, we derive an alternative expression and discuss its computation.
By employing for the confluent hypergeometric function its wellknown 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
(11) 
which can be proven to converge . One may attempt to compute with truncation
(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 doubleprecision arithmetic^{7}^{7}7The 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 doublearithmetic 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 doublearithmetic 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 doubleprecision representation. Consequently, for , these denominators are treated as ‘infinite’ and subsequent series terms remain constant, which precludes numerical convergence. The doublearithmetic C computation of series truncation (12) with takes about 45 microseconds (but does not achieve numerical convergence). On the other hand, the rationalarithmetic 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 Using Differential Equations for Infinite Series (11)
Iva 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 complexvalued function , where , is holonomic if and only if, , there exists a LODE of some order w.r.t. , as follows
(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
(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]:

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

Recast the holonomic system as a Pfaffian system.

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

Numerically evaluate at desired point , e.g., by the Runge–Kutta method.
IvB Differential Equations and Pfaffian System from Infinite Series (11) for
Theorem 2.
Function represented by the infinite series in (11) is holonomic because it satisfies the following LODEs w.r.t. and , respectively,
(15)  
(16) 
where and . These LODEs can be recast as the Pfaffian system
(17)  
(18) 
where is the dimensional vector
(19) 
with being a generic function^{8}^{8}8 stands for all four linearly independent solutions of the dimensional system of LODEs in (17)., and
(20) 
where and are polynomials in and shown in Table I.
Expression  

Proof:
Shown in Appendix A; it entails 1) byhand deduction of two partial differential equations, and 2) LODEs deduction by computeralgebrabased computations of Gröbner bases. //
IvC 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:
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 doublearithmetic 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 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 stabile^{9}^{9}9Which 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 lowerdimensional LODE system that is stabile.
Va LODE Stabile for Its Solution
Definition 5.
Consider a holonomic function^{10}^{10}10Here, 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 vectorvalued 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 theoretically^{11}^{11}11I.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.
Proof:
See Appendix BA. //
VB Direct Derivation From Integration (10) of New LODE System Stabile for
By differentiating (10) w.r.t. , we obtain the following inhomogeneous equation for :
(21) 
Now, for any integer , there exists constant such that[17]
(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)
(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:
(24) 
Theorem 5.
The dimensional LODE system (24) is stabile for when , .
Proof:
Sketched in Appendix BC. //
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 .
VC Gauge TransformationBased 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.
The proofs of the following two theorems are not included due to manuscript length limitation.
Theorem 6.
Theorem 7.
From a given LODE system that is not stabile for a target function, a lowerdimensional 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 HGMbased performance measure evaluation, also for other MIMO transceiver methods and fading models.
VD 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 C^{12}^{12}12The 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:

Maintain values within bounds of the double representation by applying to the vector from Theorem 4 the following additional gauge transformations:
(27) (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.

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

Use big float^{13}^{13}13I.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.

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 .
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 E54650 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.