5.5 Numerical stability of the method

The algorithm proposed in Section 5.1 results in a finite system of linear equations that can be solved with numerical methods. Because the linear system is a result of matrix additions, subtractions, and multiplications, its numerical stability should be examined. However, because of the presence of a linear system, and because our matrices are not M-matrices, an analytic examination of the numerical stability is not easily feasible. In this section we argue via experimentation that ETAQA-M/G/1 is numerically stable and compare its stability behavior with Ramaswami's recursive formula. Ramaswami's recursive formula for the computation of the steady state probability vector of an M/G/1-type process consists of matrix additions of non-negative matrices and these computations are known to be numerically stable ^5.3.

In the following, we focus on the stability of the method used to solve the original problem, rather than the stability of the problem itself. The latter is measured by a condition number (or conditioning), which depends on a specific instance of the problem, but not on the method used.

The stability of a method $A: \Re^N \longrightarrow \Re^M$, given an input $x\in \Re^N$, is determined as follows:

$$\Vert A(x) - A(x + \delta) \Vert < \kappa(x) \cdot \Vert\delta\Vert$$

where $\delta \in \Re^N$ is a small perturbation of the input, and $\kappa(x)$ is the conditioning of the problem with input $x$. Any norms of the corresponding vector spaces can be used, but in the following we limit our discussion to the infinity (maximum) norm. [109] states that a stable algorithm gives nearly the right answer to nearly the right question. In other words, if we change the input of a stable algorithm by a small $\delta$ we should obtain an output that is perturbed, within the constant factor $\kappa(x)$, by a corresponding amount.

We follow the above definition to examine experimentally the stability of ETAQA-M/G/1 versus that of Ramaswami's formula. The output of the aggregate scheme is a probability vector of $m + 2n$ elements and is denoted as $A(x)$, where $x$ belongs to the domain of the method, i.e., it is a choice of all the elements of the input matrices. The output of Ramaswami's is again a probability vector of $m + 2n$ elements and is denoted as $R(x)$. Note that Ramaswami's original output is post-processed to produce the same aggregate probabilities that $A(x)$ produces. We run two sets of experiments, one for a well conditioned instance of the problem, and one for an ill-conditioned instance. This is performed via the following steps:

1. Select a CTMC and solve it using both ETAQA-M/G/1 and Ramaswami's formula and check the closeness of the results.

2. Perturb the input of the selected CTMC by small random numbers. We select three different perturbation magnitudes: $10^{\rm -12}$, $10^{\rm -10}$ and $10^{\rm -6}$, and solve the CTMC with the perturbed input.

3. Repeat the perturbation experiment $50$ times with different sets of random perturbation values within the same magnitude range to achieve statistical accuracy.

4. Calculate the perturbation of output for each randomly perturbed input for ETAQA-M/G/1 solutions considering as base case the output obtained by using ETAQA-M/G/1 to solve the CTMC without any perturbation of input, i.e., $\Vert A(x)-A(x+\delta_i)\Vert$, for each experiment $i$. Calculate the same for the set of perturbed solutions using Ramaswami's formula, where the base case is the solution obtained using Ramaswami's formula on the un-perturbed input, i.e., $\Vert R(x)-R(x+\delta_i)\Vert$, for each experiment $i$..

5. Calculate the absolute difference between the solution obtained by using ETAQA-M/G/1 and Ramaswami's formula, i.e., $\Vert A(x+\delta_i) - R(x+\delta_i)\Vert$, for each experiment $i$.

Figure: Numerical behavior of algorithms under well-conditioned input. Graphics (a), (c), and (e) plot $\Vert A(x)-A(x+\delta_i)\Vert$ for ETAQA-M/G/1 and $\Vert R(x)-R(x+\delta_i)\Vert$ for Ramaswami's recursive formula for different perturbation magnitudes and for $1 \le i \le 50$ distinct experiments. Graphics (b), (d), and (f) illustrate $\Vert A(x+\delta_i) - R(x+\delta_i)\Vert$ for the same perturbation magnitudes and the same $1 \le i \le 50$ distinct experiments.

For our experiments, we selected a CTMC that models a bursty hyper-exponential server with burst sizes ranging from 1 to $p=64$. The dimension of matrices ${\mathbf{B}}$, $\L$ and ${\mathbf{F}}^{(i)}$ for $1\leq i \leq p$ is $16 \times 16$ and matrices $\widehat{{\mathbf{L}}}$, $\widehat{{\mathbf{B}}}$ a $\widehat{{\mathbf{F}}}^{(i)}$ for $1\leq i \leq p$ are of dimensions $1 \times 1$, $16 \times 1$ and $1 \times 16$ respectively. Since the corresponding ${\mathbf{G}}$ matrix for the process as well as matrices $\widehat{{\mathbf{S}}}^{(i)}$ and $\S^{(i)}$ for $1\leq i \leq p$ are full, we consider the case to a representative one.^5.4 All experiments are conducted on a Pentium III with 64-bit double precision arithmetic, and $10^{-16}$ machine precision.

Our first set of experiments considers well-conditioned input matrices, where the values of their elements differ at most by two orders of magnitude. Figure 5.4 illustrates the behavior of ETAQA-M/G/1 and Ramaswami's formula under well-conditioned input for 50 distinct experiments. Each experiment corresponds to a different $\delta_i$ but within the same magnitude range. Figure 5.4(a) shows the perturbation of solution for each of $50$ experiments for ETAQA-M/G/1 and Ramaswami's formula, is within the same magnitude range of $10^{-9}$. Observe that Figure 5.4(a) does present two lines, one for ETAQA-M/G/1 and one for Ramaswami's formula but the lines are almost indistinguishable at this level. The proximity of the two solutions is better illustrated in Figure 5.4(b) where the difference between the solutions obtained by the two different methods is plotted and is in the range of $[0.0,10^{-16}]$. The two methods are equal for all numerical purposes. Figures 5.4(c) and 5.4(e) illustrate the perturbation of solution for both methods with $\delta_i$'s in the range of $10^{-10}$ and $10^{-12}$, respectively. Across the three experiments, the degree of perturbation in the solution (i.e., the conditioning of the problem) is within three orders of magnitude less than $\delta_i$. Consistently with Figure 5.4(b), Figures 5.4(d) and 5.4(f) illustrate that the two methods agree to machine precision. Regardless of the magnitude of the input perturbation, the differences between the solutions are consistently within the same range, i.e., $10^{-16}$.

Next, we turn to a worse conditioned problem, where the elements within the various input matrices vary significantly. We use the same CTMC as the one in the previous set of experiments but the entries in all input matrices vary in magnitude up to $10^{11}$ with the largest element in the range of $10^{2}$ and the smallest in the range of $10^{-9}$. Therefore, by increasing the stiffness of the problem the possibility of numerical errors increases. Again, we perturb the input with random values within ranges of $10^{\rm -12}$, $10^{\rm -10}$, and $10^{\rm -6}$. Results are presented in Figure 5.5.

Figure: Numerical behavior of the two algorithms with ill-conditioned input. Graphics (a), (c), and (e) plot $\Vert A(x)-A(x+\delta_i)\Vert$ for ETAQA-M/G/1 and $\Vert R(x)-R(x+\delta_i)\Vert$ for Ramaswami's recursive formula for different perturbation magnitudes and for $1 \le i \le 50$ distinct experiments. Graphics (b), (d), and (f) illustrate $\Vert A(x+\delta_i) - R(x+\delta_i)\Vert$ for the same perturbation magnitudes and the same $1 \le i \le 50$ distinct experiments.

The perturbation of input matrices with values at the level of $10^{-6}$ introduces a perturbation of the solution in the range of $10^{-7}$, higher than the perturbation of solution in the well-conditioned case (compare Figures 5.5(a) and 5.4(a)). We point out that there are two lines on top of each-other in Figure 5.5(a) corresponding to ETAQA-M/G/1 and Ramaswami's output respectively. The differences between the solutions obtained by both methods for each experiments are presented in Figure 5.5(b) and are in the range of $[0.0,1.8 \times 10^{-14})$. Comparing to the results of the well-conditioned case we note an increase on the difference among the two solutions, but still very small and clearly less than the perturbation value. Figures 5.5(c) and 5.5(e) illustrate the perturbation of solutions for perturbation of inputs in the ranges of $10^{-10}$ and $10^{-12}$ respectively. Comparing to the results of Figure 5.4, we observe that the conditioning of the problem increases. The degree of perturbation remains constant for all three experiments and is one order of magnitude less than $\delta_i$, consistently across experiments. The difference of solutions between the two methods in the case of input perturbation ranges of $10^{-10}$ and $10^{-12}$ are presented in Figures 5.5(d) and 5.5(f). The differences are within the same range as for the experiment depicted in Figure 5.5(b).

The results presented in Figures 5.4 and 5.5 show that both methods, ETAQA-M/G/1 and Ramaswami's formula behave very similarly under different numerical scenarios. Since for nearly the same input we obtain, in both cases, nearly the same output, we argue that the stability of Ramaswami's recursive formula is re-confirmed. Our experiments also illustrate that ETAQA-M/G/1 and Ramaswami's recursive formula are in good agreement.