3.7.3 Fast FFT Ramaswami's formula

Meini in [58] gives an improved version of Ramaswami's formula. Once $\mbox{\boldmath {$\pi$}}^{(0)}$ is known using Eq.(3.28), the stationary probability vector is computed using matrix-generating functions associated with block triangular Toeplitz matrices^3.7. These matrix-generating functions are computed efficiently by using fast Fourier transforms (FFT).

The algorithm of the Fast FFT Ramaswami's formula is based on the fact that in practice it is not possible to store an infinite number of matrices to express the M/G/1-type process. Assuming that only $p$ matrices can be stored then the infinitesimal generator ${\mathbf{Q}}_{M/G/1}$ has the following structure

$${\mathbf{Q}}_{M/G/1} = \left[ \begin{array}{c c c c c c c c ... ...\vdots & \vdots & \vdots & \vdots & \ddots \end{array}\right].$$ (3.30)

Because there are only $p$ matrices of type $\widehat{{\mathbf{F}}}^{(i)}$ and ${\mathbf{F}}^{(i)}$, there are only $p$ sums of type $\widehat{{\mathbf{S}}}^{(i)}$ and $\S^{(i)}$ to be computed. Therefore, the computation of $\mbox{\boldmath {$\pi$}}^{(i)}$ for $i > 0$ using Ramaswami's formula, i.e., Eq.(3.26), depends only on $p$ vectors $\mbox{\boldmath {$\pi$}}^{(j)}$ for $max(0, i - p) \le j < i$. Define

$$\tilde{\mbox{\boldmath {$\pi$}}}^{(1)} = [\mbox{\boldmath {$... ...,...,\mbox{\boldmath {$\pi$}}^{(pi)}] ~~~\mbox{for}~~~i\geq 2.$$ (3.31)

The above definition simplifies the formalization of Ramaswami's formula since $\tilde{\mbox{\boldmath {$\pi$}}}^{(i)}$ depends only on the values of $\tilde{\mbox{\boldmath {$\pi$}}}^{(i-1)}$ for $i > 1$. If we apply Ramaswami's formula for vectors $\mbox{\boldmath {$\pi$}}^{(1)}$ to $\mbox{\boldmath {$\pi$}}^{(p)}$, we obtain the following equations

$$\begin{array}{c @{=} l} \mbox{\boldmath {$\pi$}}^{(1)} & - \... ...math {$\pi$}}^{(p-1)}\S^{(1)})(\S^{(0)})^{-1} \\ \end{array}.$$ (3.32)

We rewrite the above equations in the following form:

$$\begin{array}{l @{=} c} \mbox{\boldmath {$\pi$}}^{(1)}\S^{(0... ...ath {$\pi$}}^{(0)}\widehat{{\mathbf{S}}}^{(p)}\\ \end{array}.$$ (3.33)

Define

$${\mathbf{Y}}= \left[ \begin{array}{c c c c c} \S^{(0)} & \S^... ...athbf{S}}}^{(3)},\cdots,~\widehat{{\mathbf{S}}}^{(p)}\right] .$$ (3.34)

The set of equations in Eq.(3.33) can be written in a compact way by using the definitions in Eq.(3.31) and Eq.(3.34).

$$\tilde{\mbox{\boldmath {$\pi$}}}^{(1)} = - \mbox{\boldmath {$\pi$}}^{(0)} \cdot {\mathbf{b}}\cdot {\mathbf{Y}}^{-1}.$$ (3.35)

We apply Ramswami's formula for all vectors $\mbox{\boldmath {$\pi$}}^{(j)}$,  $p(i-1)+1 \leq j \leq pi$ in $\tilde{\mbox{\boldmath {$\pi$}}}^{(i)}$ for $i > 1$.

$$\begin{array}{c @{=} c @{+...+} c c} \mbox{\boldmath {$\pi$}... ...i$}}^{(p(i-1)+p-1)}\S^{(1)}) &(\S^{(0)})^{-1} \\ \end{array}.$$

These equations can be rewritten in the following form

$$\begin{array}{l@{=}l } \mbox{\boldmath {$\pi$}}^{(p(i-1)+1)}... ... -\mbox{\boldmath {$\pi$}}^{(p(i-1))}\S^{(p)} \\ \end{array}.$$

The above set of equations can be written in a matrix form as

$$\tilde{\mbox{\boldmath {$\pi$}}}^{(i)} = - \tilde{\mbox{\bol... ...}}}^{(i-1)} \cdot {\mathbf{Z}}{\mathbf{Y}}^{-1} ~~~~~ i\geq 2,$$ (3.36)

where matrix ${\mathbf{Y}}$ is defined in Eq.(3.34) and the definition of matrix ${\mathbf{Z}}$ is given by the following

$${\mathbf{Z}}= \left[ \begin{array}{c c c c c} \S^{(p)} & {\m... ...^{(2)} & \cdots & \S^{(p-1)}& \S^{(p)}\\ \end{array} \right].$$ (3.37)

The Fast Ramaswami's Formula consists of the set of equations defined in Eq.(3.35) and Eq.(3.36). The effectiveness of the representation of the Ramaswami's formula in the form of Eq.(3.35) and Eq.(3.36) comes from the special structure of matrices ${\mathbf{Y}}$ and ${\mathbf{Z}}$. The matrix ${\mathbf{Y}}$ is an upper block triangular Toeplitz matrix and the matrix ${\mathbf{Z}}$ is a lower block triangular Toeplitz matrix. Using fast Fourier transforms one can compute efficiently the inverse of a Toeplitz matrix or the multiplication of a vector with a Toeplitz matrix [58]. Although the use of Fourier transforms for matrix operations may result in numerical instabilities, in numerous test cases the above algorithm has not experienced instability [58,59].