extfieldsofstudy: Mathematics Physics

text:
Random Matrix Models, Double-Time Painlev\'e Equations, and Wireless Relaying

This paper gives an in-depth study of a multiple-antenna wireless communication scenario in which a weak signal received at an intermediate relay station is amplified and then forwarded to the final destination. The key quantity determining system performance is the statistical properties of the signal-to-noise ratio (SNR) \gamma\ at the destination. Under certain assumptions on the encoding structure, recent work has characterized the SNR distribution through its moment generating function, in terms of a certain Hankel determinant generated via a deformed Laguerre weight. Here, we employ two different methods to describe the Hankel determinant. First, we make use of ladder operators satisfied by orthogonal polynomials to give an exact characterization in terms of a"double-time"Painlev\'e differential equation, which reduces to Painlev\'e V under certain limits. Second, we employ Dyson's Coulomb Fluid method to derive a closed form approximation for the Hankel determinant. The two characterizations are used to derive closed-form expressions for the cumulants of \gamma, and to compute performance quantities of engineering interest.

I. INTRODUCTION
Over the past decade, multiple-input multiple-output (MIMO) antenna systems employing space-time coding have revolutionized the wireless industry. Such systems are well-known to offer substantial benefits in terms of channel capacity, as well as improved diversity and link reliability, and as such these techniques are being incorporated into a range of emerging industry standards. More recently, relaying strategies have been also proposed as a means of further improving the performance of space-time coded MIMO networks. Such methods are particularly effective for improving the data transmission quality of users who are located near the periphery of a communication cell.
Various MIMO relaying strategies have been proposed in the literature; see, e.g., Refs. 1-6, offering different trade-offs in various factors such as performance, complexity, channel estimation requirements, and feedback requirements. In Ref. 7, a communication technique was considered which employed a form of orthogonal space-time block coding (OSTBC), along with non-coherent amplify-and-forward (AF) processing at the multi-antenna relay. This method has the advantage of operating with only low complexity, requiring only linear processing at all terminals, achieving high diversity, and not requiring any short-term channel information at either the relays or the transmitter. For this system, an important performance measure-the so-called outage probability-which gives a fundamental probabilistic performance measure over random communication channels, was considered in Refs. 7 and 8, based on deriving an expression for the moment generating function of the received signal to noise ratio (SNR). The outage probability could then be computed by performing a subsequent numerical Laplace Transform inversion. The exact moment generating function results presented in Refs. 7 and 8, however, are quite complicated functions involving a particular Hankel determinant, and they yield little insight. Moreover, for all but small system configurations, the expressions are somewhat difficult to compute, particularly when implementing the Laplace inversion.
In this paper, we significantly expand and elaborate upon the results of Refs. 7 and 8 by employing analytical tools from random matrix theory and statistical physics. Technically, the challenge is to appropriately characterize the Hankel determinant which arises in the moment generating function computation. In general, the Hankel determinant takes the form, generated from the moments of a certain weight function w(x), 0 ≤ x < ∞, µ k := ∞ 0 x k w(x)dx , k = 0, 1, 2, . . ..
Note that the above is merely one of several possible representations for D n , with equivalent forms regularly found in the fields of statistical physics and random matrix theory.
For the problem at hand, we are faced with a Hankel determinant that is generated from a weight function of the form, which is a "two-time" deformation of the classical Laguerre weight, x α e −x . The parameters in this weight, satisfy the following conditions: For our problem both α and N s are integers, where α := |N R − N D | is the absolute difference between the number of relay and destination antennas, N R and N D respectively. The parameters α and N s are determined by the model, but in fact can be extended to take non-integer values such that α > −1 and N s > 0. 9 We will apply two different methods to characterize this Hankel determinant. For the second method, we derive an approximation for the Hankel determinant using the general linear statistics theorem obtained in Ref. 21, based on Dyson's Coulomb Fluid models. 15,[22][23][24][25] These results are essentially the Hankel analogue of asymptotic results for Toeplitz determinants. 26 Once again, these techniques have been used extensively, particularly in statistical physics, though only very recently have they been applied to address problems in communications and information theory. 13,27 It is important to note that in addition to the two methodologies advocated above, there exists other integrable systems approaches which can be used for characterizing Hankel determinants. For example, D n may be written in an equivalent matrix integral formulation, from which a 'deform-and-study' or isomonodromic deformation approach may be adopted. Essentially, this idea involves embedding D n into a more general theory of the τ-function, 28,29 using bilinear identities and linear Virasoro constraints. For details, see Refs. [30][31][32]. Yet another approach is to characterize D n through the use of Fredholm determinants as employed by Tracy and Widom. 10 Both these exact, non-perturbative integrable systems methods have their own specific advantages and disadvantages, and in most cases lead to non-linear ordinary differential or partial differential equations (ODE/PDEs) satisfied by the Hankel determinant. However, these equations are usually of higher order (equations of Chazy type usually appear), from which first integrals have to be found to reduce to a second order ODE. The advantage of the ladder operator method is that closed-form second order equations are directly obtained for a quantity related to the Hankel determinant, bypassing the need to find a first integral, or evaluate any multiple integral.
The rest of this paper is organized as follows. In Sections I A and I B, we present a brief discussion of the model which underpins the MIMO-AF wireless communication system of interest, and some basic measures which are employed to quantify system performance. Then, in Sections I C and I D, we go on to introduce the moment generating function and cumulants of interest, and pose the key mathematical problems to be dealt with in the remainder of the paper.
In Section II, we establish an exact finite n characterization of the Hankel determinant D n , employing the theory of orthogonal polynomials and their ladder operators.
In Section III, we introduce the Coulomb Fluid method, where we first give a brief overview of its key elements, following Refs. 25, 33, 13 and 21. In Section III B, we compute an asymptotic (large n) approximation for the Hankel determinant D n , and thus a corresponding characterization for the moment generating function of interest. In Section III C, the Coulomb Fluid representation for the moment generating function of the received SNR is shown to yield extremely accurate approximations for the error performance of OSTBC MIMO systems with AF relaying, even when the system dimensions are particularly small.
We also employ our analytical results to compute closed-form expressions for the cumulants of the received SNR, first via the Coulomb Fluid method in Section III D, and then presenting a refined analysis based on Painlevé equations in Sections IV and V.
Subsequently, we give an asymptotic characterization of the moment generating function, valid for scenarios for which the average received SNR is high, deriving key quantities of interest to communication engineers, including the so-called diversity order and array gain. These results are, once again, established via the Coulomb Fluid approximation in Section VI, and subsequently validated with the help of a Painlevé characterization in Section VII.

A. Amplify and Forward Wireless Relay Model
Here we briefly recall the background for the model, having been developed in Ref. 7. The dual hop MIMO communication system features a source, relay and destination terminal, having N s , N R and N D antennas respectively. In the process of transmitting a signal, each transmission period is divided into two time slots. In the first time slot, the source transmits to the relay. The relay then amplifies its received signal subject to an average power constraint, prior to transmitting the amplified signal to the destination terminal during the second time-slot. We assume that the source and destination terminals are sufficiently separated such that the direct link between them is negligible (i.e., all communication is done via the relay).
Let H 1 ∈ C N R ×N s and H 2 ∈ C N D ×N R represent the channel matrices between the source and relay, and the relay and destination terminals respectively. Each channel matrix is assumed to have uncorrelated elements distributed as C N (0, 1 A situation is considered where the source terminal transmits data using a method called OS-TBC encoding. 34 In this situation, groups of independent and identically distributed complex Gaussian random variables (referred to as information "symbols") s i , i = 1, . . ., N are assigned via a special codeword mapping to a row orthogonal matrix X = (x 1 , . . ., x N P ) ∈ C N s ×N P satisfying the power constraint E x k 2 =γ, where N P is the number of symbol periods used to send each codeword.
Since it takes N P symbol periods to transmit N symbols, the coding rate is then defined as The received signal matrix at the relay terminal at the end of the first time slot, Y = (y 1 , . . . , y N P ) ∈ C N R ×N P , is given by 1 The notation C N (µ, σ 2 ) represents a complex Gaussian distribution with mean µ and variance σ 2 .
where N ∈ C N R ×N P has uncorrelated entries distributed as C N (0, 1), representing normalized noise samples at the relay terminal. The relay then amplifies the signal it has received by a constant In the above,b is the total power constraint imposed at the relay, i.e., E Gy k 2 ≤b, whilstγ represents the average received SNR at the relay. The received signal R ∈ C N D ×N P at the destination terminal at the end of the second time slot is then given by where W ∈ C N D ×N P has uncorrelated entries distributed as C N (0, 1), representing normalized noise samples at the destination terminal.
Based on the above relationship, standard linear OSTBC decoding can be applied (see, e.g., Ref. 34). This results in decomposing the matrix model (9) into a set of parallel non-interacting single-input single-output relationships given bỹ where η l is distributed as C N (0, 1), with ||H|| F representing the Frobenius norm (or matrix norm) ofH.
The quantity that we are interested in, the instantaneous SNR for the lth symbol γ l , can then be written as

B. Wireless Communication Performance Measures
Here we recall some basic measures which are employed to quantify the performance of wireless communication systems. One of the most common measures is the so-called symbol error rate (SER), which quantifies the rate in which the transmitted symbols (or signals) are detected incorrectly at the receiver.
where M γ (·) is the moment generating function of the instantaneous SNR γ at the receiver, and Θ = π(M − 1)/M and g MPSK = sin 2 (π/M) are modulation-specific constants. As also presented in Ref. 36, the SER may be approximated in terms of M γ (·) but without the integral, via the following expression: In addition to the SER, another useful quantity is the so-called amount of fading (AoF), which serves to quantify the degree of fading (i.e., the level of randomness) in the wireless channel, and is expressed directly in terms of the cumulants of γ. Specifically, the AoF is defined as follows: where κ 1 and κ 2 are the mean and variance of γ respectively. As shown in Ref. 37, for example, this quantity can be used to describe the achievable capacity of the wireless link when the average SNR is low.
For our AF relaying system under consideration, it is clear that a major challenge is to characterize the moment generating function and also the cumulants of the instantaneous SNR given in (10). This is the key focus of the paper.

C. Statistical Characterization of the SNR γ
Here we introduce the moment generating function and cumulants of interest, and pose the key mathematical problems to be dealt with in the remainder of the paper.
Defining the positive integers, where n = 0 due to physical considerations, it was shown in Ref. 7 that the moment generating function of the instantaneous SNR (10) may be expressed as the multiple integral, From the above, we may then compute the l th cumulant of γ via Let We now consider the moment generating function (15) as a function of two variables, (T,t), or (s,t), since T = t/(1 + cs). We may then write the moment generating function (15) as We see that this involves the weight with the parameters T and t which is a deformation of the classical generalized Laguerre weight w (α) Note that w AF (x,t,t) = w (α) Lag (x). The multiple integral representation (21) is expressed as a ratio of Hankel determinants: where µ j (T,t) is the jth moment of the weight For N s ∈ N, the moments of the weight w AF are expressed in terms of the Kummer function of the second kind, The Hankel determinant for T = t, namely, can be regarded as a normalization constant, so that M γ (t,t) = 1, and its closed form expression is well-known (see Ref. 38).
From a simple change of variables, the cumulants can be calculated by transforming (16) using (20) as Note that an equivalent expression to the result given in (26) has been previously reported in Ref. 7. This result, whilst having numerical computation advantages when the system dimensions are small, has a number of shortcomings. Most importantly, its complexity does not lead to useful insight into the characteristics of the probability distribution of the SNR γ, and it does not clearly reveal the dependence on the system parameters. Moreover, in this form, the expression is not amenable to asymptotic analysis, e.g., as the number of antennas grow sufficiently large, and in such cases its numerical computation becomes unwieldy.
In the following, we will seek alternative simplified representations with the aim of overcoming these shortcomings. The key challenge is to characterize the Hankel determinant in (24 where z := t/T and M Shot-Noise (T, z) is the Hankel determinant generated via Of course, T , z and N s have different interpretations in this situation.
We mention here that the operator theory approach of Ref. 40

D. Alternative Characterization of the SNR γ
An alternative characterization for the moment generating function that will also prove to be useful is derived as follows. From the definitions of m and n, there are two possible choices for the parameter t, We first consider the sub-case N R ≤ N D . To this end, starting with (21) and with the change of variables, x i → nx i , i = 1, . . . , n, we obtain 2 where Note that in terms of s, T ′ can be written as Equivalently, we can write (32) as where and Remark 2. We introduce the variables T ′ and t ′ in order to account for the n-dependence of the variables T and t. This is important since the above representation will be useful for deriving a large n approximation for the moment generating function based on the Coulomb Fluid linear statistics approach in Section III.
For the sub-case N R > N D , we have from (31) that t is instead given by where β = m n − 1. Hence equations (37)- (39) and the results of Sections III-V are valid for the case N R > N D upon transforming t ′ and T ′ to t † and T † via and respectively. In this case, we may write T † in terms of the variable s as

FRAMEWORK
Using the ladder operator framework, and treating T and t as independent variables, the result below is proved in Appendix A. This gives a PDE satisfied by log M γ (T,t) in the variables T and t.
Theorem 1. Let the quantity H n (T,t) be defined through the Hankel determinant D n (T,t) as where the second equality follows from (24). Then H n (T,t) satisfies the following PDE: 3 : 3 Dropping the T and t dependence notation for the sake of brevity. where and We have that v = t and s = 1 c ( t T − 1) and consequently, Under the change of variables, the original PDE (46) becomes the following PDE in the variables (s, v): where A * 1 (H n ) and A * 2 (H n ) obtained from a corresponding transformation are and We note that in the PDE (54), the only higher order derivatives are the mixed derivatives with respect to s and v, and second order derivatives with respect to v.
It is of interest to note that the Hankel determinant in (21) or (24) has the following asymptotic forms: where the generating weight is the "first-time" deformation of the Laguerre weight In previous work by the authors, 13 the Hankel determinant D n [w dLag (·, T, α, N s )] was shown to arise in the analysis of the moment generating function of the single-user MIMO channel capacity, and was shown to involve a Painlevé V differential equation.

III. COULOMB FLUID METHOD FOR LARGE n ANALYSIS
In this section, we make use of the Coulomb Fluid method, which is particularly convenient when the size of the matrix n is large, to describe the MIMO-AF problem. The idea is to treat the eigenvalues as identically charged particles, with logarithmic repulsion, and held together by an external potential. When n, in this context the number of particles, is large, this assembly is regarded as a continuous fluid, first put forward by Dyson,[22][23][24] where the eigenvalues were can be found in Ref. 21 and will be used extensively in this paper. The main benefit of this approach, based on singular integral equations, is that it leads to relatively simple expressions for characterizing our moment generating function.
In Section III A, for completeness, we give a brief overview of the key elements of the Coulomb

A. Preliminaries of the Coulomb Fluid Method
We start by considering the ratio in the moment generating function expression (37), noting that it is of the form where and F n := − log Z n is known as the free energy.
This expression embraces the moment generating function expression (37) with appropriate selection of the functions v 0 (x) and f (x, T ′ ,t ′ ). A key motivation for writing our problem in this form is that it admits a very intuitive interpretation in terms of statistical physics, as originally observed by Dyson (see Refs. [22][23][24]. Specifically, if the eigenvalues x 1 , . . . , x n are interpreted as the positions of n identically charged particles, then Φ(x 1 , . . . , x n ) is recognized as the total energy of the repelling charged particles, which are confined by the external potential nv 0 (x). The function f (x, T ′ ,t ′ ) acts as a perturbation to the system, resulting in a modification to the external potential. The quantity F n (T ′ ,t ′ ) may be interpreted as the free energy of the system under an external perturbation f (x, T ′ ,t ′ ), with F n (t ′ ,t ′ ) the free energy of the unperturbed system.
For sufficiently large n, the system of particles, following Dyson, may be approximated as a continuous fluid where techniques of macroscopic physics and electrostatics can be applied. For large n we expect the external potential nv 0 (x) to be strong enough to overcome the logarithmic repulsion between the particles (or eigenvalues), and hence the particles or fluid will be confined within a finite interval to be determined through a minimization process. For this continuous fluid, we introduce a macroscopic density σ (x)dx, referred to as the equilibrium density. Since v 0 (x) is convex for x ∈ R, this density is supported on a single interval denoted by (a, b), to be determined later (see Ref. 33 for a detailed explanation). The equilibrium density is obtained by minimizing the free-energy functional: With Frostman's Lemma, 42 (p. 65) the minimizing σ (x)dx can be characterized through the integral equation where x ∈ [a, b] and A is the Lagrange multiplier for the normalization condition (64), which can be interpreted as the chemical potential of the fluid. 25,33 Noting that the integral equation above has a logarithmic kernel, taking a derivative with respect to x ∈ (a, b) converts it into a singular where P denotes Cauchy principal value.
Noting the structure (in n) of the left-hand side of (66), it is clear that σ (·) must take the general form: where σ 0 (x)dx is the density of the original system in the absence of any perturbation, while σ c (x, T ′ ,t ′ ) represents the deformation of σ 0 (x) caused by f (x, T ′ ,t ′ ). Furthermore, to satisfy Substituting (67) into (66), and comparing orders of n, we see that Following Ref. 33, where the choice for the solution for σ 0 has been extensively discussed based on the theory described in Ref. 43; the solution to (69) subject to the boundary condition σ 0 (a) = σ 0 (b) = 0 reads together with a supplementary condition, The solution to (70) Finally, the normalization condition The end points of the support of the density σ 0 (x), a, and b, are determined by (72)  With the above results, for sufficiently large n, we may approximate the ratio (60) as (see

Ref. 21 for more details)
where In the sequel, we will find explicit solutions for these quantities.
From here, the l-th cumulant can be extracted from the formula The objective of the subsequent analysis is to evaluate the quantities S AF 1 and S AF 2 . As we shall see, we will need to solve numerous integrals which are not readily available. Thus, to aid the reader, we have compiled these integrals in Appendix B.
We start by considering the end points of the support a and b. These are determined by equations (72) and (74). With the integral identities (B1)-(B3) in Appendix B, we obtain, after a few easy steps, which leads to The limiting density σ 0 (x) in (71) can be computed using the integral identity (B2) as which is the Marĉenko-Pastur law. 13,46,47 Meanwhile, with the aid of the integral identity (B4) in Using σ 0 (x) and invoking the integral identities (B5)-(B7) in Appendix B, gives Moreover, using σ c (x, T ′ ,t ′ ) and invoking the integral identities (B2), (B8) and (B12) in Appendix B gives (87)

C. SER Performance Measure Analysis Based on Coulomb Fluid
Combining (78)  The curves labeled "Coulomb (Exact SER)" were generated by substituting our Coulomb Fluid approximation into the exact expression of the SER (11), whilst the curves labeled "Coulomb (Approx SER)" are generated by substituting our Coulomb Fluid approximation into the approximate expression for the SER (12).
From these curves, it is remarkable that the Coulomb Fluid based approximations, derived under the assumption of large n, very accurately predict the SER even for small values of n. This is evident, for example, by examining the set of curves corresponding to the configuration (2, 3, 2), for which n = 2. In fact, even for the extreme cases with n = 1, the Coulomb Fluid curves still yield quite high accuracy.
We note that the results in Fig. 1  The key difference therein was that the analytic curves were based on substituting into either (11) or (12) an equivalent determinant form of the moment generating function to that given in (25) and (26). That representation, involving a determinant of a matrix with elements comprising sums of Kummer functions, is far more complicated than the Coulomb Fluid representation, which is a simple algebraic equation involving only elementary functions. The fact that the Coulomb Fluid representation yields accurate approximations for the SER for small as well as large numbers of antennas makes it a useful analytical tool for studying the performance of arbitrary MIMO-AF systems.

D. Coulomb Fluid Analysis of Large n Cumulants of SNR
In this subsection we use the Coulomb Fluid results to study the asymptotic cumulants of the SNR. We suppose β is fixed, and distinguish between two cases, β = 0 and β > 0.

Case 1: β = 0
Here N R = N D . Then from (83), a = 0 and b = 4, leading to and From the closed-form approximation of the moment generating function (78) the cumulants can be extracted using (81).
Alternatively, recall that the variable T ′ is related to s via where c =γ/(RN s ). Therefore we expand log M γ (t ′ /(1 +cs),t ′ ) obtained from the Coulomb Fluid formalism (88) and (89) in a Taylor series about s = 0. From (78), where the first few cumulants are and κ 6 (t ′ ) (96)

Case 2: β > 0
In this case, the end points of the support, a and b, are given by (83). Hence, going through the same process, the first five cumulants are and The accuracy of these approximations is demonstrated in Fig.2, where they are compared with simulation results generated by numerically computing the exact cumulants via Monte Carlo methods. As before, different antenna configurations are shown, as represented by the form . From these curves, the accuracy of our Coulomb Fluid based approximations is quite remarkable, even for small values of n.
We note that exact expressions for κ 1 and κ 2 (i.e., the mean and variance) were derived pre-

IV. POWER SERIES EXPANSION FOR H n (s, v)
In this section we seek power series solutions for (54), to determine the cumulants. We begin by removing the square roots of the PDE and arrive at the following lemma. and then the PDE (54) may be rewritten in an equivalent square-root-free form, as where A 1 and A 2 are given by (55) and (56) respectively. Now, recall from (51) that the function H n (s, v) is related to the moment generating function (21) through and note that the logarithm of the moment generating function has the expansion, (105), H n (s, v) has the following expansion in s, where ′ denotes differentiation with respect to v, and a j (v) is related to the jth cumulant κ j (v) by A. Analysis of κ 1 In this subsection, the computation of the mean κ 1 is presented. To this end, upon substituting (107) into the PDE (104), keeping the lowest power of s in the resulting expansion gives a highly nonlinear ODE expressed in the factored form where Ψ (1) 4 We introduce a j in order to make a clearer comparison to the Coulomb Fluid model and 7 without having to exactly specify a value for the constant c, since it falls out of the subsequent expansion. and Ψ (1) Hence, (109) is equivalent to either Ψ (1) (1) 1 = 0 has two types of solutions, the first of the form a 1 (v) = c 1 v + c 2 where c 1 and c 2 are constants to be determined. The second is more complicated but can be succinctly written as where c 1 and c 2 are integration constants, and f and g are algebraic functions (we omit these for brevity).
Whilst these are valid mathematical solutions, it is clear that by taking n large, they differ drastically from that predicted by the Coulomb Fluid method and, for small n, they differ with the solutions obtained in Ref. 7. This suggests that these are not the solutions of interest to our problem; thus, we set them aside and henceforth focus on Ψ (1) 2 = 0. It becomes evident at this point that we require initial conditions for the cumulants κ j (v)/(c j N s ), since the PDE (104) gives rise to a system of ODEs satisfied by κ j (v)/(c j N s ). However, (104) is not supplemented by any initial conditions from which one deduces the initial conditions at t ′ = 0 of the system of ODEs.
As we shall see, results from the Coulomb Fluid analysis (97)-(99) are crucial, as these will specify initial conditions at t ′ = 0. Consequently, the Coulomb Fluid results are in fact leading order approximations to the exact results. Without the Coulomb Fluid analysis, each cumulant κ j (v)/(c j N s ) would carry unknown constants of integration.
We now examine the equation Ψ (1) 2 = 0 in the large n region, which leads to an asymptotic characterization of a 1 (v) and thus κ 1 . To proceed, we scale the variable v by v = nt ′ . Also, note that α = n(m/n − 1) ≡ nβ . We find that a 1 (t ′ ) satisfies the following ODE: For large n, keeping the leading order term, we arrive at whose solutions are where C 1 is a constant of integration. We retain the solution that is bounded as t ′ → ∞, Comparing this with the corresponding result (97) obtained from the Coulomb fluid approximation, we see that κ CF 1 (0) = cN s n, which suggests a 1 (0) = n and hence C 1 = n(1 + β /2). Therefore, the first cumulant κ 1 to O(n) becomes κ Large n B. Analysis of κ 2 Having characterized κ 1 , we now turn to the variance κ 2 . To this end, equating the coefficients of the next lowest power of s to zero in the expansion of the PDE (104), results in an ODE which can be factored into the following form The ODE Ψ (1) 1 = 0, makes a reappearance, which we set aside, leaving the ODE Ψ where Ψ 1 , Ψ 2 and Ψ 3 are given by and Clearly, Ψ 1 depends on a 1 only, whilst Ψ 2 and Ψ 3 depend on both a 1 and a 2 . To extract the large n behavior of κ 2 , we replace v by nt ′ in (118) and make use of the large n formula for a 1 (t ′ ) in (116), resulting to a highly non-linear ODE satisfied by a 2 (t ′ ). To proceed further, keeping the highest powers of n, a first order equation is obtained for a 2 (t ′ ), Integrating this leads to an expression for a 2 (t ′ ) which, again yields an unknown integration constant. This constant can be determined through a comparison with the Coulomb Fluid results for κ 2 (t ′ ) in (98) at t ′ = 0, giving a 2 (0) = n.
The above analysis gives the leading order characterization of the variance, C. Beyond κ 1 and κ 2 The same procedure for computing κ 1 and κ 2 easily extends to higher cumulants. Here we give the leading order formula for κ 3 , κ 4 and κ 5 , and κ Large n where c j (s,t) := t α+ j+1 Γ(α + j + 1) which was derived in Ref. 7. For small fixed integer values of n and N s , e.g., n = 2, 3 and N s = 10, and n = 4, 5 and N s = 1, the above determinant can computed without much difficulty, from which the cumulants follow. It can be seen that (127), satisfied third order ODEs for a 1 (t) given by (111) and a 2 (t) given by (118). Similar results hold for the higher cumulants, which provide a consistency check.

V. LARGE n CORRECTIONS OF CUMULANTS OBTAINED FROM COULOMB FLUID
We have shown that the PDE (104) satisfied by may be used to generate a series of non-linear ODEs that are satisfied by the cumulants κ j (v).
Under a large n assumption, where v = nt ′ , the first few of these ODEs are approximated as first order ODEs for κ Large n l (t ′ ), whose solutions matched exactly with that obtained from the Coulomb Fluid analysis for the cumulants κ CF l (t ′ ). We give a method in this section where the non-linear ODEs generated from the PDE (104) are employed systematically to obtain "correction terms" to the Coulomb Fluid results.

A. Large n expansion of κ 1
We assume the first cumulant or κ 1 has the following large n expansion Substituting the above into (113), and setting the coefficients of n j to zero, a system of first order ODEs for e k (t ′ ) is obtained. These are are solved successively starting from e −1 (t ′ ), followed by e 0 (t ′ ) and so on.
As a result e −1 (t ′ ) is found to satisfy the ODE and has two solutions both satisfying the initial condition e 1 (0) = 1.
We retain the solution that is finite at infinity to match with that obtained from the Coulomb Fluid computation, namely, (97).
Continuing, we set the coefficient of n 3 to zero, which implies, Setting the coefficient of n 2 to zero gives rise to an ODE involving e ′ −1 (t ′ ), e ′′ −1 (t ′ ), e ′′′ −1 (t ′ ), e ′ 0 (t ′ ) and e ′ 1 (t ′ ). Simplifying, with e −1 (t ′ ) given by (132) and e ′ 0 (t ′ ) given by (133), we find that e 1 (t ′ ) satisfies: Setting the coefficient of the next lowest power of n to zero gives rise to an ODE involving e , and e ′ 2 (t ′ ). From the expression of e −1 (t ′ ) and e ′ 0 (t ′ ) given by (132) and (133) respectively, we find Continue with this process, the next three terms e 3 (t ′ ), e 4 (t ′ ) and e 5 (t ′ ) are found to satisfy the following equations: where l Solving ODEs (133)-(138) with initial conditions e k (0) = 0, we obtain, where A (1)

B. Large n Expansion of κ 2
In computing a series expansion for the variance, we proceed in a similar way as was just done for the mean. First, we substitute the expansion for a 1 (t ′ ) from (141) and into (118). Equating the coefficients of n k to zero, a system of first order ODEs for f k (t ′ ) are found.
In this case, setting the coefficient of n 12 to 0 results in a first order ODE for f −1 (t ′ ): The solution of this, with the initial condition f −1 (0) = 1, reads Continuing the process, we obtain an ODE for f 0 (t ′ ): The solution with the condition f 0 (0) = 0 reads from which it can be immediately seen that and we recover κ CF 2 (t ′ ) found previously. The ODEs satisfied by f k (t ′ ), k = 1, 2, 3, 4 are reported in Appendix C, equations (C1)-(C4).
These will allow us to compute f j (t ′ ) with appropriate initial conditions. Briefly, the idea is that to determine the ODE satisfied by the k-th correction term f k (t ′ ), for k odd, the preceding 1 2 (k + 1) ODEs satisfied by f j (t ′ ), j = −1, 1, 3 . . ., k are employed. For even k, the previous k 2 ODEs satisfied by f j (t ′ ), j = 0, 2, 4 . . ., k are employed. Going through the procedure described, we find that κ 2 (t ′ ) has the following large n expansion, where A (2) C. Beyond κ 1 and κ 2 The procedure adopted above for computing the large n expansion series for κ 1 (t ′ ) and κ 2 (t ′ ) easily extends to the higher cumulants. By way of example, here we focus on the third cumulant κ 3 . In this case, an asymptotic expansion for a 3 (t ′ ) is assumed, along with the initial conditions, g −1 (0) = 2 and g k (0) = 0 for k = 0, 1, 2, 3, . . ..
In this case, we find that the large n expansion of the third cumulant reads For ease of reference, the ODEs satisfied by g k (t ′ ), k = 1, 2, 3, 4, are placed in Appendix C, In summary, this entire section has shown that by expanding the cumulants κ 1 (t ′ ), κ 2 (t ′ ) and κ 3 (t ′ ) into asymptotic series in n, the Coulomb Fluid results are recovered as the leading order contributions in the large n scenario. It is also seen that, by examining the finite-n correction terms for each cumulant, no terms of O(n 2 ) or higher are present within the expansions.

VI. ASYMPTOTIC PERFORMANCE ANALYSIS BASED ON COULOMB FLUID
In this section, we return to the analysis of the moment generating function, and consider the high SNR scenario (i.e., asγ → ∞). To this end, we will study the Coulomb Fluid based approximation derived in Section III (to be complemented in Section VII through analysis based on Painlevé equations), where the variables T ′ and t ′ are taken to be dependent onγ, namely, Note that asγ → ∞, To obtain the desired high SNR expansion, we compute the moment generating function M γ as s → ∞ andγ → ∞. The Coulomb Fluid based representation (78), when expressed in terms of where S AF 1 and S AF 2 are given by (86) and (87) respectively. Our goal is to compute the expansion of M γ (s) as s → ∞ andγ → ∞, which turns out to be an expansion in (γs) −1 . It turns out that the cases β = 0 and β = 0 behave fundamentally differently, and as such, these are treated separately.

A. The Case of β = 0
With β = 0, we have N R = n. An easy computation shows that M γ admits the following expansion: where A 0 , A 1 , A 2 , A 3 , . . . are constants independent of s andγ. The leading exponent d is given by whilst the first few A ℓ are We have refrained from presenting A k , k ≥ 4, as these are rather long.

High SNR Analysis of the Symbol Error Rate (SER)
Based on (11), the SER of MPSK modulation can be expanded at high SNR using (165), resulting in where Considering the first order expansion, following Ref. 48, we may write where we identify as the so-called diversity order, and identify the factor Hence, using (176), G a admits the simplified form The high SNR results above are illustrated in Fig. 3. The "Simulation" curves are based on numerically evaluating the exact SER relation (11); the "Coulomb Fluid (Exact)" curves are based on substituting (164) into (11) and numerically evaluating the resulting integral; the "Coulomb Fluid (Leading term only)" curves are based on (174); whilst "Coulomb Fluid (Leading 4 terms)" curves are based on the first four terms of (172). The leading-order approximation is shown to give a reasonably good approximation at high SNR, whilst the additional accuracy obtained by including a few correction terms is also clearly evident.

High SNR Analysis of the Probability Density Function of γ
With the moment generating function expansion given above in (165), we may also readily obtain an approximation for the probability density function (PDF) of γ, denoted f γ (x), by direct Laplace Transform inversion. In particular, we obtain For the case of very largeγ, with the leading term gives an approximation for the PDF deep in the left-hand tail. Of course, with the inclusion of more terms, a more refined approximation is obtained.

B. The Case of β = 0 (with N R < N D )
For the situation where β = 0, two sub-cases arise. This first is N R < N D , for which N R = n; the second is N R > N D , for which N R = (1 + β )n.
In the following, we will focus on the first sub-case, N R < N D . It turns out however, that our results also apply for N R > N D upon transforming the quantity ν to ν * by where ν * (> 0) is interpreted as the (fixed) scaling factor betweenb andγ, i.e., The moment generating function given by (164) admits an expansion distinct from the β = 0 case which does not have fractional powers of 1/(γs), reading and In this situation, the sub-leading terms are very complicated, however, the jth term in the expansion can be written in the following form: where E j and F j also depend upon ν, N s , n and β .
In A 2 , E 2 , and F 2 are given by In A 3 , E 3 , and F 3 are given by In A 4 , E 4 , and F 4 are given by

High SNR Analysis of the Symbol Error Rate (SER)
Based on (11), the SER of MPSK modulation can be expanded at high SNR using (183) into where Note that here, in contrast to (173), the exponent r is a positive integer. As such, (195) admits the following closed-form solution: (Ref. 49, 2.513.1) As before, a first-order approximation is of key interest, giving where we identify the diversity order and the array gain The result for G d above is consistent with a result obtained via a different method in Ref. 8, whilst the expression for G a appears new. The high SNR results above, for the case β = 0, are illustrated in Fig. 4. As before, the "Simulation" curves are based on numerically evaluating the exact SER relation (11), and the "Coulomb Fluid (Exact)" curves are based on substituting (164) into (11) and numerically evaluating the resulting integral. Moreover, the "Coulomb Fluid (Leading term only)" curves are based on (197), whilst the "Coulomb Fluid (Leading 5 terms)" curves are based on the first five terms of (194). Again, the leading-order approximation is shown to give a reasonably good approximation at high SNR, whilst the additional accuracy obtained by including a few correction terms is also evident.

High SNR Analysis of the Probability Density Function of γ
As before, based on the moment generating function expansion (183), applying for β = 0, we can immediately take a Laplace inversion to obtain the following high SNR representation for the PDF of γ, We are mainly interested in the leading order term, and so we write the above as Note that despite the similarity with (180), interestingly, these results do not coincide upon taking β → 0 in (201), due to the differences in d and A 0 . This seems to indicate that the double asymptoticsγ → ∞ and β → 0 are non-commutative.

VII. CHARACTERIZING A 0 THROUGH PAINLEVÉ V
In this section, we obtain the leading term of the large s and largeγ expansion (183) from a Painlevé V differential equation, thus demonstrating the accuracy of the Coulomb Fluid approximation. We will focus in this section on the case β = 0.
Recall that the moment generating function of SNR γ, regarded a function of s and t, given by (21), reads, where K n,α is a normalization constant in (27). We also recall that c and t are given by respectively. So asγ → ∞, we note that A simple computation shows that M γ (s,t) admits the following expansion for largeγs: 1 (γs) nN s +1 nD n w dLag (·,t, α − N s , N s ) where w dLag (x,t, α − N s , N s ) is the deformation of the classical Laguerre weight, i.e., and The condition α − N s > −1 is required for the validity of the orthogonality relation with respect to the generalized Laguerre weight w Ref. 9). This, in turn, will ensure the validity of the first two terms in (204), while for the multiple integral in the third term to converge, the condition α − N s > 0 is required.
For the problem at hand, α and N s are integers, satisfying α ≥ 0 and N s > 0. Therefore α ≥ N s implies that the result for A 0 presented below is valid for α > 0.
Comparing the expansion (204) with (183), we see that the diversity order is and We see that A 0 is up to a constant the Hankel determinant which generates a particular Painlevé V and shows up in the single-user MIMO problem studied in Ref. 13.
We obtain A 0 , through a large n expansion of for α > 0. We will see that A 0 precisely matches that obtained in (185).
From Ref. 13 we learned that the logarithmic derivative of D n [w dLag (·,t, α − N s , N s )] with respect to t, satisfies the Painlevé V: where ′ denote derivative w.r.t. t.
We restrict to the case where N R < N D , for which N R = n. The case where N R > N D can be considered in a similar fashion as outlined in Section VI B. Setting α = nβ in the above equation, where β = m n − 1 is a fixed positive number, and with the change of variable an easy computation shows that Y n (ν) := H n (n/ν) satisfies ν 4 where ′ denotes derivative with respect to ν.
Assuming that the coefficient of n k is zero, we find that the equation c −2 = 0 gives us With MAPLE, the solutions of the above ODE for p −1 (ν) are found to be where C 2 is a constant of integration. 5 .
The equation c −1 = 0 is a coupled differential equation involving both p −1 (ν) and p 0 (ν) given by With p −1 (ν) given by (217), chosen to match the result from the Coulomb Fluid (185), we find that the first order ODE in p 0 (ν) has the solution, We disregard the second and third solutions for p −1 (ν); as these would lead to p 0 (ν) which do not generate the A 0 in agreement with that obtained from the Coulomb Fluid method.
Using a method similar to the cumulant analysis of Section V, we can also compute the first correction term to A 0 (i.e., the quantity q 1 (ν)) using (222), which reads Higher order corrections could also be obtained in a similar way.

VIII. CONCLUSION
In this paper, we have introduced two methods for characterizing the received SNR distribution in a certain MIMO communication system adopting AF relaying. We showed that the mathe- i.e., where h n is the square of the L 2 norm of P n (x). The Hankel determinant (1) is reduced to the following product: Our convention is to write P n (x) as P n (x) := x n + p 1 (n)x n−1 + p 2 (n)x n−2 + · · · + P n (0).
Hence, this implies that properties of the Hankel determinant may be obtained by characterizing the class of polynomials which are orthogonal with respect to w AF (x, T,t), over [0, ∞). It is clear that the coefficients of the polynomial P n (x), p i (n), will depend on T , t, α and N s ; for brevity, we do not display this dependence.
The main aim is to determine these unknown recurrence coefficients α n and β n from the given weight. Substituting (A4) into the three term recurrence relation results in where p 1 (0) := 0. Taking a telescopic sum gives Moreover, combining the orthogonality relationship (A2) with the three term recurrence relation leads to which can also be expressed in terms of the Hankel determinant D n in (A3) through since h n = D n+1 D n .

Ladder Operators, Compatibility Conditions, and Difference Equations
In the theory of Hermitian random matrices, orthogonal polynomials plays an important role, As we now show, there is a recursive algorithm that facilitates the determination of the recurrence coefficients α n and β n . This is implemented through the use of so-called "ladder operators" as well as their associated compatibility conditions. This approach can be traced back to Laguerre From the weight function w AF (x, T,t), one constructs the associated potential v(x) through and therefore, As shown in Ref. 15, a pair of ladder operators, to be satisfied by our orthogonal polynomials of interest, are expressed in terms of v(x) and are given by where x − y P n (y)P n−1 (y)w AF (y)dy, and where, for the sake of brevity, we have dropped the t, T dependence in w AF .
Moreover, there are associated fundamental compatibility conditions to be satisfied by A n (x) and B n (x), which are given by 16 These were initially derived for any polynomial v(x) (see Refs. [54][55][56], and then were shown to hold for all x ∈ C ∪ {∞} in greater generality. 16 We now combine (S 1 ) and (S 2 ) as follows. First, multiplying (S 2 ) by A n (x), it can be seen that the RHS is a first order difference, while (x − α n )A n (x) on the LHS can be replaced by B n+1 (x) + B n (x) + v ′ (x) from (S 1 ). Then, taking a telescopic sum with initial conditions leads to the useful identity The condition (S ′ 2 ) is of considerable interest, since it is intimately related to the logarithm of the Hankel determinant. In order to gain further information about the determinant, we need to find a way to reduce the sum to fixed number of quantities; for which, (S ′ 2 ) ultimately provides a way of going forward.

Remark 6. Since our v ′ (x) is a rational function of x, we see that
is also a rational function of x, which in turn implies that A n (x) and B n (x) are rational functions of x. Consequently, equating the residues of the simple and double pole at x = 0, x = −T , x = −t on both sides of the compatibility conditions (S 1 ), (S 2 ) and (S ′ 2 ), we obtain equations containing numerous n, T and t dependant quantities; which we call the "auxiliary variables" (to be introduced below). The resulting non-linear discrete equations are likely very complicated, but the main idea is to express the recurrence coefficients α n and β n in terms of these auxiliary variables, and eventually take advantage of the product representation (A3) to obtain an equation satisfied by the logarithmic derivative of the Hankel determinant. Now substituting (A14) into (A13) which define A n (x) and B n (x), and followed by integration by parts, we find where R * n (T,t) ≡ N s h n ∞ 0 w AF (y)P 2 n (y) y + t dy, r * n (T,t) ≡ N s h n−1 ∞ 0 w AF (y)P n (y)P n−1 (y) y + t dy, w AF (y)P 2 n (y) y + T dy, r n (T,t) ≡ N s h n−1 ∞ 0 w AF (y)P n (y)P n−1 (y) y + T dy, are the auxiliary variables.

Difference Equations from Compatibility Conditions
Inserting A n (x) and B n (x) into the compatibility conditions (S 1 ), (S 2 ) and (S ′ 2 ), and equating the residues as described above yields a system of 12 equations. Note that there are actually 14 equations, the compatibility conditions (S 1 ) and (S 2 ) also have O(x 0 ) terms which yield equations that lead to 0 = 0.
The compatibility conditions (S 1 ) give the following set of equations, Similarly, the compatibility condition (S 2 ) gives the following set of equations, To obtain a system of difference equations from the compatibility condition (S ′ 2 ) is somewhat more complicated, but carrying out the process, equating all respective residues in (S ′ 2 ) yields six equations. The first three are obtained by equating the residues of the double pole at x = 0, x = −t and x = −T respectively: − 2r n r * n − (2n − N s + α)r n +(2n + N s + α)r * n + n(n + α) = β n − β n (R n R * n−1 + R * n R n−1 ) while the last three are given by equating the residues of the simple pole at x = 0, x = −t and x = −T respectively: n−1 ∑ j=0 R j + (N s − 2r n )(r n − n − r * n ) + αr n T + (N s − 2r n )r * n + N s r n T − t + r n = − β n (R n + R n−1 + 2R n R n−1 )

Remark 7. It can be seen that (A27) is a combination of (A28) and (A29), and serves as a consis-
tency check. We also see that (A28) and (A29) respectively gives us the value of ∑ j R * j and ∑ j R j automatically in closed form. Hence we can also obtain any linear combination of these sums in closed form, which is a crucial step in obtaining a link between a linear combination of the logarithmic partial derivatives of the Hankel determinant, and the quantities β n , r n and r * n . This step would not be possible without (S ′ 2 ), also known as the bilinear identity.

Analysis of Non-linear System
Whilst some of the above difference equations (A18)-(A29) look rather complicated, our aim is to manipulate these equations in such a way as to give us insight into the recurrence coefficients α n and β n . Our dominant strategy is to always try to describe the recurrence coefficients α n and β n in terms of the auxiliary variables R n , r n , R * n and r * n . The first thing we can do is to sum (A18) through (A20) to get a simple expression for the recurrence coefficient α n in terms of the auxiliary variables R n and R * n : α n = T R n − tR * n + α + 2n + 1.
We can immediately see that by taking a telescopic sum of the above from j = 0 to j = n − 1, and recalling equation (A7), gives rise to Looking at (S 2 ) next, we sum equations (A21) through (A23) to get Taking a telescopic sum of the above from j = 0 to j = n − 1 and rearranging yields β n = Tr n − tr * n − p 1 (n).
We are now in a position to derive the following important lemma, which describes the recurrence coefficients α n and β n in terms of the set of auxiliary variables: Lemma 2. The quantities α n and β n are given in terms of the auxiliary variables R n , r n , R * n and r * n as α n = T R n − tR * n + α + 2n + 1, (A34) Proof. Equation (A34) is a restatement of equation (A30).
To obtain (A35), we use the equations obtained from(S ′ 2 ). We eliminate R * n−1 and R n−1 from (A24) using (A25) and (A26) respectively, and then rearrange to express β n in terms of the auxiliary quantities.

Toda Evolution
In this section, n is kept fixed while we vary two parameters in the weight function (22) ∂ t β n = β n (R * n − R * n−1 ).
Applying ∂ T and ∂ t to the orthogonality relation ∞ 0 w AF (x, T,t)P n (x)P n−1 (x)dx = 0 results in the following two relations: ∂ t p 1 (n) = −r * n .
Proof. To obtain (A45), we substitute in for α n in (A43) using (A34) to yield r n − r n+1 = ∂ T α n , Eliminating r n+1 in the above formula using (A20) gives 2r n = 1 − α n − T R n + T ∂ T R n − t∂ T R * n + N s .
We obtain (A46) in a similar method to (A45). This time, we first substitute in for α n in (A44) using (A34). We then proceed to eliminate r * n+1 using (A19), and finally eliminate α n again using (A34) to obtain (A46).
To obtain (A47), we differentiate equation (A33) with respect to T . Following this, we substitute in for ∂ T β n and ∂ T p 1 (n) using (A38) and (A40) respectively to yield T ∂ T r n − t∂ T r * n = β n (R n−1 − R n ).
We then proceed to replace β n R n−1 by r n (r n − N s )/R n using (A26), and finally we eliminate β n in favor of R n , R * n , r n and r * n using (A35) to obtain equation (A47). Similarly, we obtain (A48) by first differentiating (A33) with respect to t. Following this, we substitute in for ∂ t β n and ∂ t p 1 (n) using (A39) and (A41) respectively to yield T ∂ t r n − t∂ t r * n = β n (R * n − R * n−1 ).
We then proceed to replace β n R * n−1 by r * n (r * n − N s )/R * n using (A25), and finally we eliminate β n in favor of R n , R * n , r n and r * n using (A35) to obtain equation (A48), completing our proof.
From Lemma 4, the auxiliary variables r n and r * n appear linearly in equations (A45) and (A46) respectively. By eliminating r n and r * n from (A47) and (A48) using (A45) and (A46), we obtain a second order PDE system for the auxiliary variables R n and R * n . This is stated in the following lemma: Lemma 5. The auxiliary variables R n and R * n satisfy the following second-order PDE system: +2T (tR * n + 1) ∂ T R * n ) + 2t(T R n + 1) ∂ t R * n − 2T 2 R * n ∂ T R n + R n ∂ t R n − 2T (T − t) ∂ t R n +2R * n (T R n − tR * n + t) T R n − tR * n + 2n + α + 1 + (A54)

Toda Evolution of Hankel Determinant
Recall that the moment generating function is related to the Hankel determinant through equation (24): while the Hankel determinant is related to h n (T,t) by the relation D n (T,t) = ∏ n−1 j=0 h j (T,t). We compute the partial derivatives of the Hankel determinant by taking a telescopic sum from j = 0 to j = n − 1 of the partial derivatives of log h n , i.e. equations (A36) and (A37). Using (A3), To continue, we apply t∂ T + t∂ t to equation (A59), and find that (T ∂ T + t∂ t )H n = p 1 (n) + β n , where we have used (A42) to replace (T ∂ T + t∂ t )p 1 (n) by p 1 (n) + β n .
We see that equations (A59) and (A61) can be regarded as a set of simultaneous equations for p 1 (n) and β n . Solving for p 1 (n) and β n , we find that p 1 (n) = H n − n(n + α), and β n = (T ∂ T + t∂ t )H n − H n + n(n + α).
Before completing the proof of Theorem 1, we show here that our Hankel determinant D n (T,t) is related to the τ-function of a Toda PDE.
Substituting the definition (44) into (A62) and together with (A9), we find The above equation is a two parameter generalization of the Toda molecule equation, 57 and hence we identifyD n (T,t) as the corresponding τ-function of the two-parameter Toda equations.
A reduction may be obtained by writingD n (T,t) as Proof. We rewrite log(x + t 2 ) using (B11) and then use (B2) to obtain .
For our problem, where a and b are given by (83) and β ≥ 0, we see that the discriminant of the quadratic form is given by while the constant term is t 2 2 t 2 1 + 2t 1 (2 + β ) + β 2 > 0.
Hence we may invoke (B9) and then take the limit ε → 0 to get Substituting this back into (B13) and after some algebra, we get (B12).


Concept of the abstract program

The aim of this paper is to alter the abstract definition of the program of the theoretical programming model which has been developed at Eotvos Lorand University for many years in order to investigate methods that support designing correct programs. The motivation of this modi?cation was that the dynamic properties of programs appear in the model. This new definition of the program gives a hand to extend the model with the concept of subprograms while the earlier results of the original programming model are preserved.

Introduction
It is a well-known aspiration to exclusively use methods which guarantee the correctness of the developed program with respect to the problem [1,5,6] posed, that is, what makes it essential to find abstract mathematical definition of programming notions. At Eötvös Loránd University, a relational model of programming which treats the most important, fundamental notions of programming-problem, program, state space, variable, data type, etc.-in a uniform and consistent way [2,3,7] has been built for thirty years.
The base of the abstract definitions of programming notions is the state space and their tool is the mathematical relation. The problem, for example, is a homogeneous binary relation over the state space that maps from the possible initial states to the appropriate goal states. The program is defined as all of its executions, so it can be described by the relation that maps from any state to the executions starting from one, where an execution is the sequence of the states concerned. The effect of a program is also a homogeneous binary relation where its domain contains the states from which the program surely terminates, i.e. the executions starting from these states are finite, and maps to states where these executions stop. The definition of the solution (when a program is said to solve a problem) establishes the main connection between the relations of the problem and the effect of the program. This theoretical approach makes it possible to investigate the concept and correctness of several kinds of program-designer methods such as program synthesis, analogous programming, etc.
The main characteristic of this relational model is the special static point of view in which the concept of the program is a mapping and not a working, nonetheless it is originally dynamic. The advantage of this static point of view is the simple definition of the solution. Although it is not self-evident if the number of the components of the state space of a problem and that of a program are different. To prove the correctness, in this case, the program and the problem must be transformed into a new common state space [4].
In the present model however, the state space is persistent, i.e. all variables of the program are global and static (the scope and life time of the variables involve the whole program), the so-called local variables cannot be created and destroyed during an execution. In absence of local variables, the concept of real subprograms cannot be introduced into our programming model. (Only macros can be defined.) Namely, the power of a subprogram is what permits creating new local auxiliary variables as well as parameter variables, and it is not supported by the present concept of the program. It is not a serious problem if simply the correctness is wanted to be proved but it is a great difficulty if our aim is to design a program. Unfortunately, the programming method, which is based on our programming model to synthetize correct programs, cannot produce subprograms with parameter variables, so their creation remains in the sphere of the implementation instead of the designing. It is an additional disadvantage that, without the concept of subprograms, recursive calling cannot be used in our programming model. The aim of this paper is to remedy this shortcoming of the definition of the program, and to preserve the earlier results of the programming model.

State space
All definitions of the programming model are based on the state space. The concept of the state space has already been interpreted in several ways. For many people, the state space is a model of a von Neumann type computer; others, e.g. Dijkstra [1], associate this notion with the problem to be solved, where a state is a compound of the values of the main data types. So, the program is "outside" of the state space operating on it. In our programming model, this second meaning is used. In [2], the notion of the state space is a Cartesian product of the value set of data types. The only mistake of this obvious definition is that it imposes an order on the components. In [3], where the state space is a direct product, this mistake is repaired.
Here it is the definition of the direct product and the notions related to it. Let I be a finite set named as index set. Let A i 's (i P I) be arbitrary sets. The set Ś iPI A i ::" tx : I Ñ Ť iPI A i |@i P I : xpiq P A i u is the direct product of the sets A i (i P I).
The If J Ď I and @j P J : B j " A j then B is the subspace of A, i.e. B ď A. The function pr B : A Ñ B is a projection if B ď A and @a P A : pr B paq " a| J . Not only one state but a set or sequence of states can also be projected if it is done one by one. If J " tju then the function pr B is named as j variable and j is the variable name.

The new concept of the program
The most important notion of the programming model is the program. In our explanation, the program is not a collection of some statements that can be executed on a computer. The statements could only describe a program but the program is the complex of its executions. An execution is a sequence of states. A program, by definition, can always begin, i.e. at least one execution has to start from each state. If several executions start from the same start state, it means that the program is non-deterministic: nobody knows which execution will happen.
In the original programming model, every state of the progam belongs to the same state space. Now we are going to permit that the state space can be permanently changed; the inner states of the executions may have got new components.
The other novelty of our new definition is the idea of the base state space of the program. The first state (start state) of all executions and the last, if the execution is finite, are in this base state space.
Before giving the formal definition of the program, some notions must be introduced. Let H˚˚denote the set of all finite and infinite sequences of the elements of set H. H 8 includes the infinite sequences; H˚contains the finite ones. So, H˚˚" H˚Y H 8 and H˚X H 8 " H. The length of the sequence α P H˚˚is |α|, the value of which is 8 if the sequence is infinite.

Definition 1 Let A be the so-called base state space andĀ be the set of all states which belong to the state spaces B whose subspace is
2. @α P R S XĀ˚: α |α| P A (the last state of the finite executions of the range of S), 3. @a P A and @α P Spaq : |α| ě 1 and α 1 " a.
The variables of the base state space are the base variables; the other variables are the auxiliary variables of the program.
The concept of the program allows to create and destroy new components (variables) during an execution, so the state space changes dynamically. All new components have to be destroyed at the termination, at the very latest, but the base variables should never be removed. The current state always contains the components of the base state space.
The following three definitions will outline that the base state space of any program is denoted arbitrarily. The base variables of a program can be extended, restricted or renamed without changing the essence of the program.
The base variables of any program can be renamed easily without the execution of the program is changed. However, if the new name of a base variable would be identical to an auxiliary variable's name then this auxiliary variable is given a new, unique name. of S onto C is the program (denoted by S 1 ) whose base state space is C and for all c P C : S 1 pcq " tγ P C˚˚| Dα P Sppr A pcqq : |α| " |γ|^@i P r2..|γ|s : if k R D α i then γ i | I " α i and γ i pkq " γ i´1 pkq if k P D α i then γ i " α i u.
It is easy to generalize this definition in case the state space A is a subspace of the state space C.
The following definition shows how a base variable can be degraded to an auxiliary variable: the first step creates it and the last step destroys it.

Definition 4 Let S be a program over the state space A "
Ś iPI A i and let C be the subspace of A. The restriction of S onto C is the program (denoted by S) whose base state space is C and for all c P C: The base state space of a program can be renamed, extended or restricted several times, thus the base state space can be totally transformed. But these transformations can change only the base state space and not the program. The essence of the executions of the program remains the same. The renamed, extended and restricted program hardly differs from its original version.
Definition 5 Two programs are identical if they can transform into the same program, using extensions, restrictions and renaming.

Conclusions
We have introduced a new concept of the program. Now let us observe what kinds of effects this modification has on the original programming model.

Concept of dynamic program
The new concept of the program smuggles the dynamic property of the program back to the static programming model [2,3] developed at Eötvös Loránd University. Certainly, the fact that a program can create and destroy auxiliary variables is not a brilliant discovery. The difficulty of our investigation was to find out how this dynamic property can be embedded into the programming model described in the static point of view.
The most important element of our concept is the base state space. It can be seen as the interface of the program. It determines the components through which the program can communicate, and can make contact with its environment. Only the base variables can get value from outside before the program starts, and their value can be asked after the termination.
The program becomes very flexible because its base state space can be changed easily. This property is reflected in the fact that the same program can possess different interfaces depending on the problem to be solved. The aim of a program can be changed by altering its base state space whereas its operation cannot be changed.

Concept of solution
Despite the above modification, the definition of the solution can be preserved [2,3,7]. The only thing that must be done is to redefine the concept of the effect of the program because the definition of the solution relies on it.
Now we will repeat all definitions that are important to describe the concept of solution.
Let the problem generally be a relation F Ď AˆA where A is a state space. Let S be a program over A.
The program S is said to solve the problem F if The definition of the solution supposes that the state space of the problem is identical to the base state space of the program. As the base state space can be chosen arbitrarily, the program can be extended, restricted or renamed in order that its base state space is identical to the state space of the problem. Since programs of this kind are identical, if one of them can solve the problem, so can the others.
Moreover, if the effects of the programs S 1 and S 2 are identical over all common base state spaces (ppS 1 q " ppS 2 q) and if one of them can solve a problem, so can the others. In this case these programs are called equivalent. This relation is reflexive, symmetrical and transitive, thus it is an equivalence relation. Consequently, if one program belonging to an equivalence class can solve a problem, then every program derived from this class can also solve it.
In the original programming model, the definition of the solution must be generalized [4] when the state spaces of the problem and the program are different. Now, our new concept is used to avoid this. It is enough to fit the base state space of the program to the state space of the problem and transform the program onto this common state space.

Earlier results
The original programming model contains many important results that can be used in the verification of programs. It is apparent that these results are not based directly on the concept of the program; except for the effect of the program. The definition of the effect of the program has not been changed because the new definition of the program fixes that the start state and the end state of the finite execution are in the same state space (that is, the base state space). Accordingly, all earlier results, namely Dijkstra's weakest precondition [1], the theorem of the specification or the derivation rules [1,3] of the program constructions [3] are used in unalterable forms.
Certainly, in the definitions of program constructions, the components (programs and conditions) which are the parts of the construction have to be on the same base state space. Accordingly, before making one of the constructions, this common base state space has to be agreed on.

Subprograms
The concept of real subprograms, which have got parameter variables, could not be introduced into the original programming model because the program has got an unvarying state space, thus local variables cannot be created during the execution, so parameter variables, which are local variables, cannot be used. Now, the situation is changed. In our new programming model, the concept of the subprograms may be defined.
First and foremost the subprogram is a program. It can be executed independently; it has got own base state space. Its only speciality is that it can be built into another program. (This is the host program.) Each program, including subprograms, is equivalent to an assignment. This assignment is appropriate to identify the subprogram that is equivalent to it. Accordingly, in an arbitrary program description language, a subprogram can be denoted with this assignment. The variables of this assignment, which are the formal parameters, form the base state space of the subprogram. The variables at the left-hand side of the assignment are the output parameter variables; the ones at the right-hand side are the input parameter variables. (For simplicity's sake, we shall restrict our consideration to assignments where a parameter variable does not occur more than once.) This assignment can be also used as a calling statement in the description of the host program. Certainly, in the calling statement, the parameter variables can be substituted by actual parameters (arguments), which are the expressions (often, they are only variables) of the host program. An output variable cannot be changed by an expression more general than a variable. The number and the type of the formal and current parameters must be the same. At this calling statement, the execution of the host program is interrupted, the control switches over the subprogram, the values of current parameters are given to the formal parameter variables, and the subprogram is executed. After the termination of the subprogram, the values of the output parameter variables are recopied to the appropriate current parameter variables. (Here, two kinds of parameter passing have been defined: passing by value and passing by value-result but other passing methods may be introduced.) The base state space of the subprogram is the interface between the subprogram and the host program. At the beginning of the subprogram, the variables of this base state space will get their initial values from the current state of the host program. At the end of it these variables give their values to the variables of the host program. The parameter variables and other local variables of the subprogram are created when the subprogram is called and destroyed at the termination of the subprogram.
In addition, the usage of subprograms permits recursive callings because a subprogram can call itself. After each calling, the parameter and the local variables of the subprogram are created again and again without the error message "out of memory" because the memory of the abstract program is infinite.
During planning, it is convenient that all the current variables of the host program are treated as global with respect to the subprogram that is called by the host program. The applied programming style and the facilities of the selected programming language determine if these global variables can be used directly in the subprogram or not. Obviously, their use has to be forbidden if we want to guarantee the independency of the subprogram.
We can make difference between the calling statement and the calling expression. In the latter case, only the right-hand side of the subprogram's head appears inside an expression of the host program, which contains actual parameters instead of formal parameter variables. After the termination, the result value of the left-hand side of the head is given back to the place of the calling expression in the host program.
Certainly, the abstract program description language which is used for planning has to be extended with the notation of the subprogram and of its calling including the connection between the actual parameters and the parameter variables.
To sum up the introduction of the concept of the subprograms makes it possible to design subprograms during planning and not only during implementation.


A Novel Feature-based Bayesian Model for Query Focused Multi-document Summarization

Both supervised learning methods and LDA based topic model have been successfully applied in the field of query focused multi-document summarization. In this paper, we propose a novel supervised approach that can incorporate rich sentence features into Bayesian topic models in a principled way, thus taking advantages of both topic model and feature based supervised learning methods. Experiments on TAC2008 and TAC2009 demonstrate the effectiveness of our approach.

Introduction
Query-focused multi-document summarization (Nenkova et al., 2006;Wan et al., 2007;Ouyang et al., 2010) can facilitate users to grasp the main idea of documents. In query-focused summarization, a specific topic description, such as a query, which expresses the most important topic information is proposed before the document collection, and a summary would be generated according to the given topic.
Supervised models have been widely used in summarization (Li, et al., 2009, Shen et al., 2007, Ouyang et al., 2010. Supervised models usually regard summarization as a classification or regression problem and use various sentence features to build a classifier based on labeled negative or positive samples. However, existing supervised approaches seldom exploit the intrinsic structure among sentences. This disadvantage usually gives rise to serious problems such as unbalance and low recall in summaries. Recently, LDA-based (Blei et al., 2003) Bayesian topic models have widely been applied in multidocument summarization in that Bayesian approaches can offer clear and rigorous probabilistic interpretations for summaries (Daume and Marcu, 2006;Haghighi and Vanderwende, 2009;Jin et al., 2010;Mason and Charniak, 2011;Delort and Alfonseca, 2012). Exiting Bayesian approaches label sentences or words with topics and sentences which are closely related with query or can highly generalize documents are selected into summaries. However, LDA topic model suffers from the intrinsic disadvantages that it only uses word frequency for topic modeling and can not use useful text features such as position, word order etc (Zhu and Xing, 2010). For example, the first sentence in a document may be more important for summary since it is more likely to give a global generalization about the document. It is hard for LDA model to consider such information, making useful information lost.
It naturally comes to our minds that we can improve summarization performance by making full use of both useful text features and the latent semantic structures from by LDA topic model. One related work is from Celikyilmaz and Hakkani-Tur (2010). They built a hierarchical topic model called Hybhsum based on LDA for topic discovery and assumed this model can produce appropriate scores for sentence evaluation. Then the scores are used for tuning the weights of various features that helpful for summary generation. Their work made a good step of combining topic model with feature based supervised learning. However, what their approach confuses us is that whether a topic model only based on word frequency is good enough to generate an appropriate sentence score for regression. Actually, how to incorporate features into LDA topic model has been a open problem. Supervised topic models such as sLDA (Blei and MacAuliffe 2007) give us some inspiration. In sLDA, each document is associated with a labeled feature and sLDA can integrate such feature into LDA for topic modeling in a prin-cipled way.
With reference to the work of supervised LDA models, in this paper, we propose a novel sentence feature based Bayesian model S-sLDA for multidocument summarization. Our approach can naturally combine feature based supervised methods and topic models. The most important and challenging problem in our model is the tuning of feature weights. To solve this problem, we transform the problem of finding optimum feature weights into an optimization algorithm and learn these weights in a supervised way. A set of experiments are conducted based on the benchmark data of DUC2007, TAC2008 and TAC2009, and experimental results show the effectiveness of our model.
The rest of the paper is organized as follows. Section 2 describes some background and related works. Section 3 describes our details of S-sLDA model. Section 4 demonstrates details of our approaches, including learning, inference and summary generation. Section 5 provides experiments results and Section 6 concludes the paper.

Related Work
A variety of approaches have been proposed for query-focused multi-document summarizations such as unsupervised (semi-supervised) approaches, supervised approaches, and Bayesian approaches.
Unsupervised (semi-supervised) approaches such as Lexrank (Erkan and Radex, 2004), manifold (Wan et al., 2007) treat summarization as a graphbased ranking problem. The relatedness between the query and each sentence is achieved by imposing querys influence on each sentence along with the propagation of graph. Most supervised approaches regard summarization task as a sentence level two class classification problem. Supervised machine learning methods such as Support Vector Machine(SVM) (Li, et al., 2009), Maximum Entropy (Osborne, 2002) , Conditional Random Field (Shen et al., 2007 and regression models (Ouyang et al., 2010) have been adopted to leverage the rich sentence features for summarization.
Recently, Bayesian topic models have shown their power in summarization for its clear probabilistic interpretation. Daume and Marcu (2006) proposed Bayesum model for sentence extraction based on query expansion concept in information retrieval. Haghighi and Vanderwende (2009) proposed topicsum and hiersum which use a LDA-like topic model and assign each sentence a distribution over background topic, doc-specific topic and content topics. Celikyilmaz and Hakkani-Tur (2010) made a good step in combining topic model with supervised feature based regression for sentence scoring in summarization. In their model, the score of training sentences are firstly got through a novel hierarchical topic model. Then a featured based support vector regression (SVR) is used for sentence score prediction. The problem of Celikyilmaz and Hakkani-Turs model is that topic model and feature based regression are two separate processes and the score of training sentences may be biased because their topic model only consider word frequency and fail to consider other important features. Supervised feature based topic models have been proposed in recent years to incorporate different kinds of features into LDA model. Blei (2007) proposed sLDA for document response pairs and Daniel et al. (2009) proposed Labeled LDA by defining a one to one correspondence between latent topic and user tags. Zhu and Xing (2010) proposed conditional topic random field (CTRF) which addresses feature and independent limitation in LDA. The hierarchical Bayesian LDA (Blei et al., 2003) models the probability of a corpus on hidden topics as shown in Figure 1(a). Let K be the number of topics , M be the number of documents in the corpus and V be vocabulary size. The topic distribution of each document θ m is drawn from a prior Dirichlet distribution Dir(α), and each document word w mn is sampled from a topic-word distribution φ z specified by a drawn from the topic-document distribution θ m . β is a K × M dimensional matrix and each β k is a distribution over the V terms. The generating procedure of LDA is illustrated in Figure 2. θ m is a mixture proportion over topics of document m and z mn is a K dimensional variable that presents the topic assignment distribution of different words.
Supervised LDA (sLDA) (Blei and McAuliffe 2007) is a document feature based model and intro- duces a response variable to each document for topic discovering, as shown in Figure 1(b). In the generative procedure of sLDA, the document pairwise la-

Problem Formulation
Here we firstly give a standard formulation of the task. Let K be the number of topics, V be the vocabulary size and M be the number of documents.
where N m denotes the number of sentences in m th document. Each sentence is represented with a collection of words {w msn } n=Nms n=1 where N ms denotes the number of words in current sentence.
− − → Y ms denotes the feature vector of current sentence and we assume that these features are independent.

S-sLDA
z ms is the hidden variable indicating the topic of current sentence. In S-sLDA, we make an assumption that words in the same sentence are generated from the same topic which was proposed by Gruber (2007). z msn denotes the topic assignment of current word. According to our assumption, z msn =  Figure 4: generation process for S-sLDA z ms for any n ∈ [1, N ms ]. The generative approach of S-sLDA is shown in Figure 3 and Figure 4. We can see that the generative process involves not only the words within current sentence, but also a series of sentence features. The mixture weights over features in S-sLDA are defined with a generalized linear model (GLM).
Here we assume that each sentence has T features and − − → Y ms is a T × 1 dimensional vector. η is a K × T weight matrix of each feature upon topics, which largely controls the feature generation procedure. Unlike s-LDA where η is a latent variable estimated from the maximum likelihood estimation algorithm, in S-sLDA the value of η is trained through a supervised algorithm which will be illustrated in detail in Section 3.

Posterior Inference and Estimation
Given a document and labels for each sentence, the posterior distribution of the latent variables is: Eqn.
(2) cannot be efficiently computed. By applying the Jensens inequality, we obtain a lower bound of the log likelihood of document where H(q) = −E[logq] and it is the entropy of variational distribution q is defined as here γ a K-dimensional Dirichlet parameter vector and multinomial parameters. The first, third and forth terms of Eqn.
(3) are identical to the corresponding terms for unsupervised LDA (Blei et al., 2003). The second term is the expectation of log probability of features given the latent topic assignments.
The Bayes estimation for S-sLDA model can be got via a variational EM algorithm. In EM procedure, the lower bound is firstly minimized with respect to γ and φ, and then minimized with α and β by fixing γ and φ. E-step: The updating of Dirichlet parameter γ is identical to that of unsupervised LDA, and does not involve where Ψ(·) denotes the log Γ function. m s denotes the document that current sentence comes from and Y st denotes the t th feature of sentence s.

M-step:
The M-step for updating β is the same as the procedure in unsupervised LDA, where the probability of a word generated from a topic is proportional to the number of times this word assigned to the topic. In this subsection, we describe how we learn the feature weight η in a supervised way. The learning process of η is a supervised algorithm combined with variational inference of S-sLDA. Given a topic description Q 1 and a collection of training sentences S from related documents, human assessors assign a score v(v = −2, −1, 0, 1, 1) to each sentence in S. The score is an integer between −2 (the least desired summary sentences) and +2 (the most desired summary sentences), and score 0 denotes neutral atti- is the set containing sentences with score v. Let φ Qk denote the probability that query is generated from topic k. Since query does not belong to any document, we use the following strategy to leverage φ Qk (9) In Equ.(9), w∈Q β kw denotes the probability that all terms in query are generated from topic k and 1 can be seen as the average probability that all documents in the corpus are talking about topic k. Eqn. (9) is based on the assumption that query topic is relevant to the main topic discussed by the document corpus. This is a reasonable assumption and most previous LDA summarization models are based on similar assumptions.
Next, we define φ Ov,k for sentence set O v , which can be interpreted as the probability that all sentences in collection O v are generated from topic k.
|O v | denotes the number of sentences in set O v . Inspired by the idea that desired summary sentences would be more semantically related with the query, we transform problem of finding optimum η to the following optimization problem: T t=1 η kt = 1 (11) 1 We select multiple queries and their related sentences for training where KL(O v ||Q) is the Kullback-Leibler divergence between the topic and sentence set O v as shown in Eqn. (12).
In Eqn. (11), we can see that O 2 , which contain desirable sentences, would be given the largest penalty for its KL divergence from Query. The case is just opposite for undesired set.
Our idea is to incorporate the minimization process of Eqn.(11) into variational inference process of S-sLDA model. Here we perform gradient based optimization method to minimize Eqn.(11). Firstly, we derive the gradient of L(η) with respect to η.
For simplification, we regard β and γ as constant during updating process of η, so ∂φ Qk ∂ηxy = 0. 2 We can further get first derivative for each labeled sentence.

Feature Space
Lots of features have been proven to be useful for summarization (Louis et al., 2010  In addition, we also use the commonly used features including sentence position, paragraph position, sentence length and sentence bigram frequency.

E-step
initialize φ 0 sk := 1/K for all i and s. initialize γ mi := α mi + N )m/K for all i. initialize η kt = 0 for all k and t.

Sentence Selection Strategy
Next we explain our sentence selection strategy. According to our intuition that the desired summary should have a small KL divergence with query, we propose a function to score a set of sentences Sum. We use a decreasing logistic function ζ(x) = 1/(1+ e x ) to refine the score to the range of (0,1).
Let Sum denote the optimum update summary. We can get Sum by maximizing the scoring function.

Score(Sum)
1. Learning: Given labeled set O v , learn the feature weight vector η using algorithm in Figure 5. 2. Given new data set and η, use algorithm in section 3.3 for inference. (The only difference between this step and step (1) is that in this step we do not need minimize L(η). 3. Select sentences for summarization from algorithm in Figure 6. A greedy algorithm is applied by adding sentence one by one to obtain Sum . We use G to denote the sentence set containing selected sentences. The algorithm first initializes G to Φ and X to SU . During each iteration, we select one sentence from X which maximize Score(s m ∪ G). To avoid topic redundancy in the summary, we also revise the MMR strategy (Goldstein et al., 1999;Ouyang et al., 2007) in the process of sentence selection. For each s m , we compute the semantic similarity between s m and each sentence s t in set Y in Eqn. (18).
We need to assure that the value of semantic similarity between two sentences is less than T h sem . The whole procedure for summarization using S-sLDA model is illustrated in Figure 6. T h sem is set to 0.5 in the experiments.

Experiments Set-up
The query-focused multi-document summarization task defined in DUC 3 (Document Understanding Conference) and TAC 4 (Text Analysis Conference) evaluations requires generating a concise and well organized summary for a collection of related news documents according to a given query which describes the users information need. The query usually consists of a title and one or more narrative/question sentences. The system-generated summaries for DUC and TAC are respectively limited to 250 words and 100 words. Stop-words in both documents and queries are removed using a stop-word list of 598 words, and the remaining words are stemmed by Porter Stemmer 6 . As for the automatic evaluation of summarization, ROUGE (Recall-Oriented Understudy for Gisting Evaluation) measures, including ROUGE-1, ROUGE-2, and ROUGE-SU4 7 and their corresponding 95% confidence intervals, are used to evaluate the performance of the summaries. In order to obtain a more comprehensive measure of summary quality, we also conduct manual evaluation on TAC data with reference to (Haghighi and Vanderwende, 2009;Celikyilmaz and Hakkani-Tur, 2011;Delort and Alfonseca, 2011).

Comparison with other Bayesian models
In this subsection, we compare our model with the following Bayesian baselines: KL-sum: It is developed by Haghighi and Vanderwende (Lin et al., 2006) by using a KLdivergence based sentence selection strategy.
where P s is the unigram distribution of candidate summary and Q d denotes the unigram distribution of document collection. Sentences with higher ranking score is selected into the summary. HierSum: A LDA based approach proposed by Haghighi and Vanderwende (2009), where unigram distribution is calculated from LDA topic model in Equ.(14).
For fair comparison, baselines use the same proprecessing methods with our model and all sum-maries are truncated to the same length of 100 words.
From Table 1 and  see that among all the Bayesian baselines, Hybhsum achieves the best result. This further illustrates the advantages of combining topic model with supervised method. In Table 1, we can see that our S-sLDA model performs better than Hybhsum and the improvements are 3.4% and 3.7% with respect to ROUGE-2 and ROUGE-SU4 on TAC2008 data. The comparison can be extended to TAC2009 data as shown in Table 2: the performance of S-sLDA is above Hybhsum by 4.3% in ROUGE-2 and 5.1% in ROUGE-SU4. It is worth explaining that these achievements are significant, because in the TAC2008 evaluation, the performance of the top ranking systems are very close, i.e. the best system is only 4.2% above the 4th best system on ROUGE-2 and 1.2% on ROUGE-SU4.

Comparison with other baselines.
In this subsection, we compare our model with some widely used models in summarization. Manifold: It is the one-layer graph based semisupervised summarization approach developed by Wan et al.(2008). The graph is constructed only considering sentence relations using tf-idf and neglects topic information.
LexRank: Graph based summarization approach (Erkan and Radev, 2004), which is a revised version of famous web ranking algorithm PageRank. It is an unsupervised ranking algorithms compared with Manifold.
SVM: A supervised method -Support Vector Machine (SVM) (Vapnik 1995) which uses the same features as our approach.
MEAD: A centroid based summary algorithm by Radev et al. (2004). Cluster centroids in MEAD consists of words which are central not only to one article in a cluster, but to all the articles. Similarity is measure using tf-idf.
At the same time, we also present the top three participating systems with regard to ROUGE-2 on TAC2008 and TAC2009 for comparison, denoted as (denoted as SysRank 1st, 2nd and 3rd) (Gillick et al., 2008;Zhang et al., 2008;Gillick et al., 2009;Varma et al., 2009). The ROUGE scores of the top TAC system are directly provided by the TAC evaluation.
From Table 3 and Table 4, we can see that our approach outperforms the baselines in terms of ROUGE metrics consistently. When compared with the standard supervised method SVM, the relative improvements over the ROUGE-1, ROUGE-2 and ROUGE-SU4 scores are 4.3%, 13.1%, 8.3% respectively on TAC2008 and 7.2%, 14.9%, 14.3% on TAC2009. Our model is not as good as top participating systems on TAC2008 and TAC2009. But considering the fact that our model neither uses sentence compression algorithm nor leverage domain knowledge bases like Wikipedia or training data, such small difference in ROUGE scores is reasonable.

Manual Evaluations
In order to obtain a more accurate measure of summary quality for our S-sLDA model and Hybhsum, we performed a simple user study concerning the following aspects: (1) Overall quality: Which summary is better overall? (2) Focus: Which summary contains less irrelevant content? (3)Responsiveness: Which summary is more responsive to the query. (4) Non-Redundancy: Which summary is less redundant? 8 judges who specialize in NLP participated in the blind evaluation task. Evaluators are presented with two summaries generated by S-sLDA   and Hybhsum, as well as the four questions above. Then they need to answer which summary is better (tie). We randomly select 20 document collections from TAC 2008 data and randomly assign two summaries for each collection to three different evaluators to judge which model is better in each aspect. As we can see from  Responsiveness, S-sLDA model outputs Hybhsum based on t-test on 95% confidence level. Table 6 shows the example summaries generated respectively by two models for document collection D0803A-A in TAC2008, whose query is "Describe the coal mine accidents in China and actions taken". From table 6, we can see that each sentence in these two summaries is somewhat related to topics of coal mines in China. We also observe that the summary in Table 6(a) is better than that in Table 6(b), tending to select shorter sentences and provide more information. This is because, in S-sLDA model, topic modeling is determined simultaneously by various features including terms and other ones such as sentence length, sentence position and so on, which can contribute to summary quality. As we can see, in Table 6(b), sentences (3) and (5) provide some unimportant information such as "somebody said", though they contain some words which are related to topics about coal mines.
(1)China to close at least 4,000 coal mines this year: official (2)By Oct. 10 this year there had been 43 coal mine accidents that killed 10 or more people, (3)Officials had stakes in coal mines. (4)All the coal mines will be closed down this year. (5) In the first eight months, the death toll of coal mine accidents rose 8.5 percent last year. (6) The government has issued a series of regulations and measures to improve the coun.try's coal mine safety situation. (7)The mining safety technology and equipments have been sold to countries. (8)More than 6,000 miners died in accidents in China (1) In the first eight months, the death toll of coal mine accidents across China rose 8.5 percent from the same period last year. (2)China will close down a number of ill-operated coal mines at the end of this month, said a work safety official here Monday.
(3) Li Yizhong, director of the National Bureau of Production Safety Supervision and Administration, has said the collusion between mine owners and officials is to be condemned.
(4)from January to September this year, 4,228 people were killed in 2,337 coal mine accidents. (5) Chen said officials who refused to register their stakes in coal mines within the required time

Conclusion
In this paper, we propose a novel supervised approach based on revised supervised topic model for query-focused multi document summarization. Our approach naturally combines Bayesian topic model with supervised method and enjoy the advantages of both models. Experiments on benchmark demonstrate good performance of our model.