A note on integrating products of linear forms over the unit simplex
aa r X i v : . [ c s . PF ] M a y A note on integrating products of linear forms over the unitsimplex
Giuliano CasaleDepartment of ComputingImperial College London, UK [email protected]
July 4, 2018
Abstract
Integrating a product of linear forms over the unit simplex can be done in polynomial time ifthe number of variables n is fixed [1]. In this note, we highlight that this problem is equivalent toobtaining the normalizing constant of state probabilities for a popular class of Markov processesused in queueing network theory. In light of this equivalence, we survey existing computationalalgorithms developed in queueing theory that can be used for exact integration. For example,under some regularity conditions, queueing theory algorithms can exactly integrate a product oflinear forms of total degree N by solving N systems of linear equations. Let ∆ = { x ∈ R n | x i ≥ , P ni =1 x i = 1 } be the unit simplex and denote by dm the integral Lebesguemeasure on ∆ . Denote by θ , θ , . . . , θ d a collection of linear forms on R n , where θ ij is the i th elementof θ j and define the coefficient matrix θ = ( θ ij ) ∈ R nd . Also let N , . . . , N d be a set of nonnegativeintegers, with N = P di =1 N i and N = ( N , . . . , N d ) .Recently, [1] considers the problem of integrating a homogeneous polynomial function f of degree N over ∆ , observing that this can be reduced to computing a finite collection of integrals of the type J ( θ , N ) = Z ∆ d Y j =1 n X i =1 θ ij x i ! N j dm (1)A polynomial-time algorithm to compute (1) is then proposed in [1], which obtains this integral bydetermining the coefficient of z N · · · z N d d in the Taylor expansion of T ( z , . . . , z d ) = 1 Q ni =1 (1 − P dj =1 z j θ ij ) (2)in O ( N d ) time. Our goal is to examine alternative methods to compute (1) that arise from queueingnetwork theory [4].We first observe that (2) corresponds to the generating function of the normalizing constant G ( θ , N ) of state probabilities for a class of finite Markov processes known as product-form closedqueueing networks , which have been extensively used in performance evaluation of computer and1ommunication systems [2,21]. We provide in Appendix A a brief introduction to this class of Markovmodels; additional background may be found in [4, 28].The connection between normalizing constants and (1) is explicitly noted in [11], which showsthat J ( θ , N ) = N ! N ! · · · N d !( N + n − G ( θ , N ) (3)Therefore, algorithms developed in queueing theory to compute G ( θ , N ) may be readily used tointegrate products of linear forms over the unit simplex.The rest of this note reviews exact computational methods for G ( θ , N ) . Note that queueing theoryimplicitly assumes θ ij ≥ , since the coefficients θ ij represent processing times of jobs, however thiscondition is not used in the algorithms reviewed later.In the following sections, we denote by N − j the vector obtained from N by decreasing its j thelement N j by one unit. We also indicate with θ − θ j the matrix obtained from θ by removing a rowwith elements θ j , . . . , θ jD , and similarly denote by θ + θ j the matrix obtained by appending to θ anew row with elements θ j , . . . , θ jD . Recurrence relations are the standard method used in queueing theory to compute G ( θ , N ) . Existingmethods differ for computational requirements and ease of implementation. We here review twoclassic methods, convolution and RECAL, and the method of moments, a class of algorithms thatscales efficiently under large degree N . We point to [12, 20, 21, 26, 27] for other recursive algorithms. The convolution algorithm is a recurrence relation for normalizing constants presented in its mostgeneral form in [29], extending the case d = 1 first developed in [5]. The method relies on therecurrence relation G ( θ , N ) = G ( θ − θ n , N ) + d X j =1 θ nj G ( θ , N − j ) (4)with termination conditions G ( ∅ , · ) = 0 , G ( · , ) = 1 , where = (0 , . . . , , and G ( · , N − j ) = 0 if N j = 0 . It can be verified that the algorithm requires O ( N d ) time and space for fixed n and d . Variants of this algorithm exist [12, 27], for example a specialized version has been proposed forproblems with sparse θ [26]. Compared to convolution, the RECAL algorithm is more efficient on models with large d [14]. Themethod is based on the recurrence relation G ( θ , N ) = N − d n X i =1 θ id G ( θ + θ i , N − d ) (5)with similar termination conditions as the convolution algorithm. The computational complexity is O ( N n ) time and space for fixed n and d . The method is well-suited for parallel implementation andcan be optimized for sparse θ ij coefficients [17]. 2 .1.3 Method of moments The method of moments simultaneously applies (4) and (5) to a set of normalizing constants that differfor the elements and composition of θ and N . Different rules exist to define this set of normalizingconstants, called the basis , leading to multiple variants of the method [8–10]. Assume that the nor-malizing constants in the basis are grouped in a vector g ( N ) , then the method of moments defines amatrix recurrence relation A ( θ , N ) g ( N ) = B ( θ , N ) g ( N − j ) (6)for any ≤ j ≤ D and where g ( ) can be computed explicitly from the termination conditionsof (4)-(5). Here A ( θ , N ) and B ( θ , N ) are square matrices, with constant or decreasing sizes asthe recurrence progresses, depending on the implementation. The coefficients within these matricesdepend on j . Thus, to avoid to recompute these matrices at each step, the method of moments firstperforms a recursion along dimension j = d , then along j = d − , and so forth up to reaching N = . An explicit example of the basic method may be found in [7]. Provided that A ( θ , N ) isinvertible at all steps of the recursion, it is possible to solve (6) and obtain G ( θ , N ) from the elementsof g ( N ) . A necessary condition for invertibility is that θ ir = θ jr , ∀ i, j, r .Compared to existing algorithms, (6) finds G ( θ , N ) after solving N systems of linear equations.Since the order of A ( θ , N ) does not depend on N , the theoretical complexity of the method is O ( N ) time and O (1) space for fixed n and d . However, implementation complexity is larger due to the costof exact algebra, which is normally required to solve (6) without error propagation [6]. d = 1 In this section we consider the case d = 1 , where we can simplify notation by indicating N with N and θ i with θ i . It has been known since long time that in the case d = 1 , and provided that all θ i coefficients are distinct, the normalizing constant can be written as [23] G ( θ , N ) = n X i =1 θ N + n − i Q k = i ( θ k − θ i ) (7)that is simply the divided difference [ θ , . . . , θ n ] x N + n − . It is useful to note that in the case wheresome coefficients θ ij coincide, (7) generalizes to [3] G ( θ , N ) = [ θ i , . . . , θ n ] x N + n − = n ′ X j =1 ( − m j − θ N + n − m j j X r ≥ r = m j − ( − r j (cid:18) N + r j r j (cid:19) n ′ Y k =1 k = j (cid:18) m k + r k − r k (cid:19) θ r k k ( θ j − θ k ) m k + r k (8)where n ′ ≤ n is the number of distinct coefficients θ i , and m j denotes the number of coefficientsidentical to θ j . 3 .2.2 Case d > For arbitrary values of d , [11] derives some generalizations of (7). The first one is G ( θ , N ) = X ≤ t ≤ N ( − N − t N ! · · · N d ! d Y j =1 (cid:18) N j t j (cid:19) n X i =1 ( P dj =1 t j θ ij ) N + n − Q k = i ( P dj =1 t j ( θ kj − θ ij )) (9)where t = ( t , . . . , t d ) , t = P dj =1 t j , which has the same O ( N d ) time complexity of the convolutionalgorithm, but O (1) space requirements. This expression holds assuming that θ ij values for given j are all distinct, otherwise the following generalized expression based on (8) should be used [11] G ( θ , N ) = X ≤ t ≤ N ( − N − t N ! · · · N R ! R Y r =1 (cid:18) N r t r (cid:19) n ′ X j =1 ( − m j − ( P Rs =1 t s θ js ) N + n − m j × X r ≥ r = m j − ( − r j (cid:18) N + r j r j (cid:19) n ′ Y k =1 k = j (cid:18) m k + r k − r k (cid:19) ( P Rs =1 t s θ ks ) r k ( P Rs =1 t s ( θ js − θ ks )) m k + r k (10)Another explicit expression is given by G ( θ , N ) = X h ≥ : h ≤ N ( − N − h N ! · · · N d ! (cid:18) N + n − N − h (cid:19) d Y j =1 n X i =1 h i θ ij ! N j (11)where h = ( h , . . . , h n ) , h = P ni =1 h i . This expression has O ( N n +1 ) time and O (1) space require-ments.It is also worth noting that since the integrand of J ( θ , N ) is a polynomial, the cubature rulesproposed in [18] provide alternative explicit expressions for G ( θ , N ) , which become exact for asufficiently high interpolation degree, corresponding to a O (( N/ d ) time complexity. We pointto [3, 15, 16] for other works on explicit expressions for the normalizing constant. In this note, we have highlighted a connection between queueing network theory and the integrationof products of linear forms over the unit simplex as in (1). In order to solve queueing network models,a normalizing constant is required, which can be computed with the recurrence relations and theexplicit expressions that we have reviewed. This normalizing constant readily provides the value ofthe integral (1).While the scope of the present note is restricted to exact methods, numerical approximations andasymptotic expansions are also available for G ( θ , N ) , such as [13, 22, 24, 25, 30]. References [1] V. Baldoni, N. Berline, J. A. d. Loera, M. K ¨oppe, and M. Vergne. How to integrate a polynomialover a simplex.
Mathematics of Computation , 80:297–325, 2011.[2] F. Baskett, K. M. Chandy, R. R. Muntz, and F. G. Palacios. Open, closed, and mixed networksof queues with different classes of customers.
JACM , 22:248–260, 1975.43] A. Bertozzi and J. McKenna. Multidimensional residues, generating functions, and their appli-cation to queueing networks.
SIAM Review , 35(2):239–268, 1993.[4] G. Bolch, S. Greiner, H. de Meer, and K. S. Trivedi.
Queueing Networks and Markov Chains .Wiley, 2006.[5] J. P. Buzen. Computational algorithms for closed queueing networks with exponential servers.
Comm. of the ACM , 16(9):527–531, 1973.[6] G. Casale. An efficient algorithm for the exact analysis of multiclass queueing networks withlarge population sizes.
Proc. of ACM SIGMETRICS , 169–180, 2006.[7] G. Casale. An application of exact linear algebra to capacity planning models.
ACM Communi-cations in Computer Algebra , 42(4):202–205, 2008.[8] G. Casale. CoMoM: Efficient class-oriented evaluation of multiclass performance models.
IEEETrans. on Software Engineering , 35(2):162–177, 2009.[9] G. Casale. A generalized method of moments for closed queueing networks.
PerformanceEvaluation , 68(2):180–200, 2011.[10] G. Casale. Exact analysis of performance models by the method of moments.
PerformanceEvaluation , 68(6):487–506, 2011.[11] G. Casale. Accelerating performance inference over closed systems by asymptotic methods.
Proc. ACM Meas. Anal. Comput. Syst (POMACS) , 1(1), 2017. To be presented at ACM SIG-METRICS 2017, preprint available at https://spiral.imperial.ac.uk/handle/10044/1/43431 .[12] K. M. Chandy and C. H. Sauer. Computational algorithms for product-form queueing networksmodels of computing systems.
Comm. of the ACM , 23(10):573–583, 1980.[13] G L. Choudhury, K. K. Leung, and W. Whitt. Calculating normalization constants of closedqueuing networks by numerically inverting their generating functions.
JACM , 42(5):935–970,1995.[14] A. E. Conway and N. d. Georganas. RECAL - A new efficient algorithm for the exact analysisof multiple-chain closed queueing networks.
Journal of the ACM , 33(4):768–791, 1986.[15] A. I. Gerasimov. On normalizing constants in multiclass queueing networks.
Oper. Res. ,43(4):704–711, 1995.[16] J. J. Gordon. The evaluation of normalizing constants in closed queueing networks.
Oper. Res. ,38(5):863–869, 1990.[17] A. G. Greenberg and J. McKenna. Solution of closed, product form, queueing networks via theRECAL and TREE-RECAL methods on a shared memory multiprocessor.
ACM SIGMETRICSPerformance Evaluation Review (PER) , 17(1):127–135, 1989.[18] A. Grundmann, H.M. M¨oller. Invariant integration formulas for the n-simplex by combinatorialmethods.
SIAM J. on Numerical Analysis , 15(2):282–290, 1978.[19] P. G. Harrison. On normalizing constants in queueing networks.
Operations Research ,33(2):464–468, 1985. 520] P. G. Harrison and Ting Ting Lee A New Recursive Algorithm for Computing GeneratingFunctions in Closed Multi-Class Queueing Networks.
Proc. of IEEE MASCOTS , 231–238, 2004.[21] P. G. Harrison and S. Coury. On the asymptotic behaviour of closed multiclass queueing net-works.
Performance Evaluation , 47(2):131–138, 2002.[22] C. Knessl, C. Tier. Asymptotic expansions for large closed queueing networks with multiple jobclasses.
IEEE Trans. Computers , 41(4):480–488, 1992.[23] E. Koenigsberg. Cyclic queues.
Operational Research Quarterly , 9, 1:22–35, 1958.[24] Y. Kogan. Asymptotic expansions for probability distributions in large loss and closed queueingnetworks.
Perform. Eval. Rev. , 29(3):25–27, Dec. 2001.[25] Y. Kogan, A. Yakovlev. Asymptotic analysis for closed multichain queueing networks withbottlenecks.
QUESTA , 23:235–258, 1996.[26] S. Lam. Dynamic scaling and growth behavior of queueing network normalization constants.
Journal of the ACM , 29(2):492–513, 1982.[27] S. Lam. A simple derivation of the
MVA and
LBANC algorithms from the convolution algorithm.
IEEE Trans. on Computers , 32:1062–1064, 1983.[28] S. S. Lavenberg. A perspective on queueing models of computer performance.
PerformanceEvaluation , 10(1):53–76, 1989.[29] M. Reiser and H. Kobayashi. Queueing networks with multiple closed chains: Theory andcomputational algorithms.
IBM J. Res. Dev. , 19(3):283–294, 1975.[30] W. Wang, G. Casale, C. Sutton. A Bayesian Approach to Parameter Inference in QueueingNetworks.
ACM TOMACS , 27(1), 2016.
A Product-form closed queueing networks
In closed queueing networks, a finite set of N jobs circulate among a network of n nodes, calledqueueing stations, according to a set of transition probabilities. Jobs are partitioned into d types,called classes, and the network topology is closed, meaning that the number of jobs inside the networkremains constant over time and equal to N j in class j . Upon visiting a station i , a job of class j isprocessed at rate µ ij = θ − ij k − i when a total of k i jobs are simultaneously present at the station.Therefore, the larger the number of jobs at a station, the smaller the rate.A common question that arises in these models is to determine the joint stationary distribution ofthe number k ij of class- j jobs residing at station i when the Markov process reaches equilibrium. Theanalytic form of this distribution is known to be [2] P ( k ) = 1 G θ ( N ) n Y i =1 k i ! Q dj =1 k ij ! d Y l =1 θ k il il (12)where k ∈ S is a state vector, S is a state space given by S = ( k ∈ R nd (cid:12)(cid:12)(cid:12) k ij ≥ , n X i =1 k ij = N j ) P k ∈S P ( k ) = 1 , i.e. , G ( θ , N ) = X k ∈S n Y i =1 k i ! Q dj =1 k ij ! d Y l =1 θ k il ilil