A note on integrating products of linear forms over the unit simplex

TL;DR

This paper reveals the equivalence between integrating products of linear forms over the unit simplex and computing normalizing constants in queueing networks, highlighting existing algorithms for exact integration.

Contribution

It establishes a novel connection between polynomial integration and queueing theory, enabling the use of queueing algorithms for exact integration tasks.

Findings

Integration can be performed in polynomial time for fixed variables

Queueing algorithms can compute exact integrals under certain conditions

N systems of linear equations suffice for degree N polynomial integration

Abstract

Integrating a product of linear forms over the unit simplex can be done in polynomial time if the number of variables n is fixed (V. Baldoni et al., 2011). In this note, we highlight that this problem is equivalent to obtaining the normalizing constant of state probabilities for a popular class of Markov processes used in queueing network theory. In light of this equivalence, we survey existing computational algorithms 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 of linear forms of total degree N by solving N systems of linear equations.

Equations

$$J(\mathbf{\theta },\mathbf{N})={\text{bigintsss}}_{\Delta }\prod _{j=1}^d{\left(\sum _{i=1}^n{\theta }_{ij}{x}_i\right)}^{N_j}dm$$

$$T(z_1,\dots ,z_d)=\frac{1}{{\prod }_{i=1}^n(1-{\sum }_{j=1}^d{z}_j{\theta }_{ij})}$$

$$J(\mathbf{\theta },\mathbf{N})=\frac{N_1!N_2!\cdots N_d!}{(N+n-1)!}G(\mathbf{\theta },\mathbf{N})$$

$$G(\mathbf{\theta },\mathbf{N})=G(\mathbf{\theta }-{\theta }_n,\mathbf{N})+\sum _{j=1}^d{\theta }_{nj}G(\mathbf{\theta },\mathbf{N}-1_j)$$

$$G(\mathbf{\theta },\mathbf{N})=N_d^{-1}\sum _{i=1}^n{\theta }_{id}G(\mathbf{\theta }+{\theta }_i,\mathbf{N}-1_d)$$

$$\mathbf{A}(\mathbf{\theta },\mathbf{N})\mathbf{g}(\mathbf{N})=\mathbf{B}(\mathbf{\theta },\mathbf{N})\mathbf{g}(\mathbf{N}-1_j)$$

$$G(\mathbf{\theta },N)=\sum _{i=1}^n\frac{{\theta }_i^{N+n-1}}{{\prod }_{k\ne i}({\theta }_i-{\theta }_k)}$$

$$G(\mathbf{\theta },N)$$

$$=\sum _{j=1}^{n^{\prime }}(-1{)}^{m_j-1}{\theta }_j^{N+n-m_j}\sum _{\begin{array}{c}\mathbf{r}\ge 0\\ r=m_j-1\end{array}}(-1{)}^{r_j}\left(\genfrac{}{}{0ex}{}{N+r_j}{r_j}\right)\prod _{\genfrac{}{}{0ex}{}{k=1}{k\ne j}}^{n^{\prime }}\left(\genfrac{}{}{0ex}{}{m_k+r_k-1}{r_k}\right)\frac{{\theta }_k^{r_k}}{({\theta }_j-{\theta }_k{)}^{m_k+r_k}}$$

$$G(\mathbf{\theta },\mathbf{N})=\sum _{0\le \mathbf{t}\le \mathbf{N}}\frac{(-1{)}^{N-t}}{N_1!\cdots N_d!}\prod _{j=1}^d\left(\genfrac{}{}{0ex}{}{N_j}{t_j}\right)\sum _{i=1}^n\frac{({\sum }_{j=1}^d{t}_j{\theta }_{ij}{)}^{N+n-1}}{{\prod }_{k\ne i}({\sum }_{j=1}^d{t}_j({\theta }_{ij}-{\theta }_{kj}))}$$

$$\begin{gathered} G(\mathbf{\theta },\mathbf{N})=\sum _{0\le \mathbf{t}\le \mathbf{N}}\frac{(-1{)}^{N-t}}{N_1!\cdots N_R!}\prod _{r=1}^R\left(\genfrac{}{}{0ex}{}{N_r}{t_r}\right)\sum _{j=1}^{n^{\prime }}(-1{)}^{m_j-1}({\sum }_{s=1}^R{t}_s{\theta }_{js}{)}^{N+n-m_j} \\ \times {\sum }_{\begin{array}{c}\mathbf{r}\ge 0\\ r=m_j-1\end{array}}(-1{)}^{r_j}(\genfrac{}{}{0ex}{}{N+r_j}{r_j}\left){\prod }_{\genfrac{}{}{0ex}{}{k=1}{k\ne j}}^{n^{\prime }}\right(\genfrac{}{}{0ex}{}{m_k+r_k-1}{r_k})\frac{({\sum }_{s=1}^R{t}_s{\theta }_{ks}{)}^{r_k}}{({\sum }_{s=1}^R{t}_s({\theta }_{js}-{\theta }_{ks}){)}^{m_k+r_k}} \end{gathered}$$

$$G(\mathbf{\theta },\mathbf{N})=\sum _{\genfrac{}{}{0ex}{}{\mathbf{h}\ge 0:}{h\le N}}\frac{(-1{)}^{N-h}}{N_1!\cdots N_d!}\left(\genfrac{}{}{0ex}{}{N+n-1}{N-h}\right)\prod _{j=1}^d{\left(\sum _{i=1}^n{h}_i{\theta }_{ij}\right)}^{N_j}$$

$$\mathbb{P}(\mathbf{k})=\frac{1}{G_{\mathbf{\theta }}(\mathbf{N})}\prod _{i=1}^n\frac{k_i!}{{\prod }_{j=1}^d{k}_{ij}!}\prod _{l=1}^d{\theta }_{il}^{k_{il}}$$

$${\cal S}=\Biggl{\{}\boldsymbol{k}\in\mathbb{R}^{nd}\;\Big{|}\Bigr{.}\;k_{ij}\geq 0,\;\sum_{i=1}^{n}k_{ij}=N_{j}\Biggr{\}}$$

$$G(\mathbf{\theta },\mathbf{N})=\sum _{\mathbf{k}\in \mathcal{S}}\prod _{i=1}^n\frac{k_i!}{{\prod }_{j=1}^d{k}_{ij}!}\prod _{l=1}^d{\theta }_{il}^{k_{il}}$$

Full text

1 Introduction

Let $\Delta=\{\boldsymbol{x}\in\mathbb{R}^{n}|x_{i}\geq 0,\sum_{i=1}^{n}x_{i}=1\}$ be the unit simplex and denote by $dm$ the integral Lebesgue measure on $\Delta$. Denote by $\theta_{1},\theta_{2},\ldots,\theta_{d}$ a collection of linear forms on $\mathbb{R}^{n}$, where $\theta_{ij}$ is the $i$th element of $\theta_{j}$ and define the coefficient matrix $\boldsymbol{\theta}=(\theta_{ij})\in\mathbb{R}^{nd}$. Also let $N_{1}$,$\ldots$,$N_{d}$ be a set of nonnegative integers, with $N=\sum_{i=1}^{d}N_{i}$ and $\boldsymbol{N}=(N_{1},\ldots,N_{d})$.

Recently, [1] considers the problem of integrating a homogeneous polynomial function $f$ of degree $N$ over $\Delta$, observing that this can be reduced to computing a finite collection of integrals of the type

[TABLE]

A polynomial-time algorithm to compute (1) is then proposed in [1], which obtains this integral by determining the coefficient of $z_{1}^{N_{1}}\cdots z_{d}^{N_{d}}$ in the Taylor expansion of

[TABLE]

in $O(N^{d})$ time. Our goal is to examine alternative methods to compute (1) that arise from queueing network theory [4].

We first observe that (2) corresponds to the generating function of the normalizing constant $G({\boldsymbol{\theta}},\boldsymbol{N})$ of state probabilities for a class of finite Markov processes known as product-form closed queueing networks, which have been extensively used in performance evaluation of computer and communication systems [2, 21]. We provide in Appendix A a brief introduction to this class of Markov models; additional background may be found in [28, 4].

The connection between normalizing constants and (1) is explicitly noted in [11], which shows that

[TABLE]

Therefore, algorithms developed in queueing theory to compute $G({\boldsymbol{\theta}},\boldsymbol{N})$ may be readily used to integrate products of linear forms over the unit simplex.

The rest of this note reviews exact computational methods for $G({\boldsymbol{\theta}},\boldsymbol{N})$. Although in queueing theory it is assumed that $\theta_{ij}\geq 0$, computational algorithms developed therein for $G({\boldsymbol{\theta}},\boldsymbol{N})$ do not rely on this assumption and can therefore be applied also to linear forms with arbitrary values of $\theta_{ij}$.

In the following sections, we denote by $\boldsymbol{N}-\boldsymbol{1}_{j}$ the vector obtained from $\boldsymbol{N}$ by decreasing its $j$th element $N_{j}$ by one unit. We also indicate with $\boldsymbol{\theta}-\theta_{j}$ the matrix obtained from $\boldsymbol{\theta}$ by removing a row with elements $\theta_{j1},\ldots,\theta_{jD}$, and similarly denote by $\boldsymbol{\theta}+\theta_{j}$ the matrix obtained by appending to $\boldsymbol{\theta}$ a new row with elements $\theta_{j1},\ldots,\theta_{jD}$.

2 Exact computational algorithms

2.1 Recurrence relations

Recurrence relations are the standard method used in queueing theory to compute $G({\boldsymbol{\theta}},\boldsymbol{N})$. Existing methods differ for computational requirements and ease of implementation. We here review two classic methods, convolution and RECAL, and the method of moments, a class of algorithms that scales efficiently under large degree $N$. We point to [21, 20, 12, 26, 27] for other recursive algorithms.

2.1.1 Convolution algorithm

The convolution algorithm is a recurrence relation for normalizing constants presented in its most general form in [29], extending the case $d=1$ first developed in [5]. The method relies on the recurrence relation

[TABLE]

with termination conditions $G(\emptyset,\cdot)=0$, $G(\cdot,\boldsymbol{0})=1$, where $\boldsymbol{0}=(0,\ldots,0)$, and $G(\cdot,\boldsymbol{N}-\boldsymbol{1}_{j})=0$ if $N_{j}=0$. It can be verified that the algorithm requires ${\cal 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 for problems with sparse $\boldsymbol{\theta}$ [26].

2.1.2 Recurrence by chain algorithm (RECAL)

Compared to convolution, the RECAL algorithm is more efficient on models with large $d$ [14]. The method is based on the recurrence relation

[TABLE]

with similar termination conditions as the convolution algorithm. The computational complexity is ${\cal O}(N^{n})$ time and space for fixed $n$ and $d$. The method is well-suited for parallel implementation and can be optimized for sparse $\theta_{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 differ for the elements and composition of $\boldsymbol{\theta}$ and $\boldsymbol{N}$. Different rules exist to define this set of normalizing constants, called the basis, leading to multiple variants of the method [8, 10, 9]. Assume that the normalizing constants in the basis are grouped in a vector $\boldsymbol{g}(\boldsymbol{N})$, then the method of moments defines a matrix recurrence relation

[TABLE]

for any $1\leq j\leq D$ and where $\boldsymbol{g}(\boldsymbol{0})$ can be computed explicitly from the termination conditions of (4)-(5). Here $\boldsymbol{A}(\boldsymbol{\theta},\boldsymbol{N})$ and $\boldsymbol{B}(\boldsymbol{\theta},\boldsymbol{N})$ are square matrices, with constant or decreasing sizes as the recurrence progresses, depending on the implementation. The coefficients within these matrices depend on $j$. Thus, to avoid to recompute these matrices at each step, the method of moments first performs a recursion along dimension $j=d$, then along $j={d-1}$, and so forth up to reaching $\boldsymbol{N}=\boldsymbol{0}$. An explicit example of the basic method may be found in [7]. Provided that $\boldsymbol{A}(\boldsymbol{\theta},\boldsymbol{N})$ is invertible at all steps of the recursion, it is possible to solve (6) and obtain $G(\boldsymbol{\theta},\boldsymbol{N})$ from the elements of $\boldsymbol{g}(\boldsymbol{N})$. A necessary condition for invertibility is that $\theta_{ir}\neq\theta_{jr}$, $\forall i,j,r$.

Compared to existing algorithms, (6) finds $G(\boldsymbol{\theta},\boldsymbol{N})$ after solving $N$ systems of linear equations. Since the order of $\boldsymbol{A}(\boldsymbol{\theta},\boldsymbol{N})$ does not depend on $N$, the theoretical complexity of the method is ${\cal O}(N)$ time and ${\cal O}(1)$ space for fixed $n$ and $d$. However, implementation complexity is larger due to the cost of exact algebra, which is normally required to solve (6) without error propagation [6].

2.2 Explicit solutions

2.2.1 Case $d=1$

In this section we consider the case $d=1$, where we can simplify notation by indicating $\boldsymbol{N}$ with $N$ and $\theta_{i1}$ with $\theta_{i}$. It has been known since long time that in the case $d=1$, and provided that all $\theta_{i}$ coefficients are distinct, the normalizing constant can be written as [23]

[TABLE]

that is simply the divided difference $[\theta_{1},\ldots,\theta_{n}]x^{N+n-1}$. It is useful to note that in the case where some coefficients $\theta_{ij}$ coincide, (7) generalizes to [3]

[TABLE]

where $n^{\prime}\leq n$ is the number of distinct coefficients $\theta_{i}$, and $m_{j}$ denotes the number of coefficients identical to $\theta_{j}$.

2.2.2 Case $d>1$

For arbitrary values of $d$, [11] derives some generalizations of (7). The first one is

[TABLE]

where $\boldsymbol{t}=(t_{1},\ldots,t_{d})$, $t=\sum_{j=1}^{d}t_{j}$, which has the same ${\cal O}(N^{d})$ time complexity of the convolution algorithm, but ${\cal O}(1)$ space requirements. This expression holds assuming that $\theta_{ij}$ values for given $j$ are all distinct, otherwise the following generalized expression based on (8) should be used [11]

[TABLE]

Another explicit expression is given by

[TABLE]

where $\boldsymbol{h}=(h_{1},\ldots,h_{n})$, $h=\sum_{i=1}^{n}h_{i}$. This expression has ${\cal O}(N^{n+1})$ time and ${\cal O}(1)$ space requirements.

It is also worth noting that since the integrand of $\boldsymbol{J}(\boldsymbol{\theta},\boldsymbol{N})$ is a polynomial, the cubature rules proposed in [18] provide alternative explicit expressions for $G({\boldsymbol{\theta}},\boldsymbol{N})$, which become exact for a sufficiently high interpolation degree, corresponding to a ${\cal O}((N/2)^{d})$ time complexity. We point to [16, 15, 3] for other works on explicit expressions for the normalizing constant.

3 Conclusion

In this note, we have highlighted a connection between queueing network theory and the integration of 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 the explicit expressions that we have reviewed. This normalizing constant readily provides the value of the integral (1).

While the scope of the present note is restricted to exact methods, numerical approximations and asymptotic expansions are also available for $G({\boldsymbol{\theta}},\boldsymbol{N})$, such as [13, 30, 22, 25, 24].

4 Acknowledgement

The author thanks Yosef Rinott for helpful comments on an earlier version of this manuscript.

Appendix A Product-form closed queueing networks

In closed queueing networks, a finite set of $N$ jobs circulate among a network of $n$ nodes, called queueing 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 network remains constant over time and equal to $N_{j}$ in class $j$. Upon visiting a station $i$, a job of class $j$ is processed at rate $\mu_{ij}=\theta^{-1}_{ij}k^{-1}_{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 of the number $k_{ij}$ of class-$j$ jobs residing at station $i$ when the Markov process reaches equilibrium. The analytic form of this distribution is known to be [2]

[TABLE]

where $\boldsymbol{k}\in{\cal S}$ is a state vector, ${\cal S}$ is a state space given by

[TABLE]

and the normalizing constant is obtained by requiring that $\sum_{\boldsymbol{k}\in{\cal S}}\mathbb{P}(\boldsymbol{k})=1$, i.e.,

[TABLE]