Simplify $ \sum_{\{x_{ij}|\forall_i \sum_j x_{ij} = \xi_i, \forall_j \sum_i x_{ij} = \eta_j \}} \prod_{ij} \frac{a_{ij}^{x_{ij}}}{x_{ij}!} $ Let $\{x_{ij}\}$ be a finite set of nonnegative integer variables, with $i = 1..m$ and $j = 1..n$. Let $a_{ij}, \xi_i, \eta_j \geq 0$. Here $\xi_i , \eta_j$ are integers, but $a_{ij}$ 
can be a real number.
Simplify the following sum:


*

*$$ \sum_{\{x_{ij}|\forall_i \sum_j x_{ij} = \xi_i, \forall_j \sum_i
    x_{ij} = \eta_j \}} \prod_{ij} \frac{a_{ij}^{x_{ij}}}{x_{ij}!} $$


That is, the sum goes over all values of the variables $\{x_{ij}\}$ constrained by $\sum_j x_{ij} = \xi_i$ for all $i = 1..m$, and $\sum_i x_{ij} = \eta_j$ for all $j = 1..n$. In other words, the matrix with entries $x_{ij}$ has constrained sums of rows and columns. (See here and here)
By simplifying, I mean find an equivalent expression that that is easier to compute. For example, by the multinomial theorem,
$$
\sum_{\{x_{ij}|\forall_i \sum_j x_{ij} = \xi_i\}} \prod_{ij} \frac{a_{ij}^{x_{ij}}}{x_{ij}!} 
= \sum_i \frac{1}{\xi_i!} \left(\sum_j a_{ij}\right)^{\xi_i}
$$
However in (1.) the sum has two constraints, $\sum_j x_{ij} = \xi_i$ and $\sum_i x_{ij} = \eta_j$, so the multinomial theorem cannot be directly applied. I am looking for a similar simplification of (1.).
Any suggestions are appreciated.
 A: We can use generating functions to find an equivalent expression.
Define:
$$f(\alpha_1,\ldots,\alpha_m;\beta_1,\ldots,\beta_n) = \sum_{\{x_{ij}\}} \prod_{ij} \frac{(a_{ij}\alpha_i \beta_j)^{x_{ij}}}{x_{ij}!},$$
where the sum over $x_{ij}$ is unconstrained, except that the $x_{ij}$ must be non-negative integers. It is easy to see that:
$$f(\alpha_1,\ldots,\alpha_m;\beta_1,\ldots,\beta_n) = \exp\left(\sum_{ij} a_{ij} \alpha_i \beta_j\right)$$
The summation (1.) in the question equals the coefficient of
$$\prod_{ij}\alpha_i^{\xi_i}\beta_j^{\eta_j}$$
in the Taylor series of $f(\alpha_1,\ldots,\alpha_m;\beta_1,\ldots,\beta_n)$. This is:
$$
\frac{1}{\prod_i \xi_i!}\frac{1}{\prod_j \eta_j!}
\left.\frac{\partial^{\xi_1 + \cdots + \xi_m}}{\partial\alpha_1^{\xi_1}\ldots\partial\alpha_m^{\xi_m}}
\frac{\partial^{\eta_1 + \cdots + \eta_n}}{\partial\beta_1^{\eta_1}\ldots\partial\beta_n^{\eta_n}}
\exp\left(\sum_{ij} a_{ij} \alpha_i \beta_j\right)
\right|_{\alpha_1 = \cdots = \alpha_m = \beta_1 = \cdots = \beta_n = 0}$$
I am not marking this as the answer. Hopefully someone comes up with something better.
(See also: https://scicomp.stackexchange.com/q/7833/988)
