Computing exponentials of essentially non-negative matrices entrywise to high relative accuracy

Abstract A real square matrix is said to be essentially non-negative if all of its off-diagonal entries are non-negative. It has recently been shown that the exponential of an essentially non-negative matrix is determined entrywise to high relative accuracy by its entries up to a condition number intrinsic to the exponential function (Numer. Math., 110 (2008), pp. 393-403). Thus the smaller entries of the exponential may be computed to the same relative accuracy as the bigger entries. This paper develops algorithms to compute exponentials of essentially non-negative matrices entrywise to high relative accuracy. Keywords: matrix exponential; essentially non-negative matrix; high relative accuracy algorithms

Computing the Exponential of Large Block-Triangular Block-Toeplitz Matrices Encountered in Fluid Queues

The Erlangian approximation of Markovian fluid queues leads to the problem of computing the matrix exponential of a subgenerator having a block-triangular, block-Toeplitz structure. To this end, we propose some algorithms which exploit the Toeplitz structure and the properties of generators. Such algorithms allow to compute the exponential of very large matrices, which would otherwise be untreatable with standard methods. We also prove interesting decay properties of the exponential of a generator having a block-triangular, block-Toeplitz structure

Componentwise accurate fluid queue computations using doubling algorithms

Markov-modulated fluid queues are popular stochastic processes frequently used for modelling real-life applications. An important performance measure to evaluate in these applications is their steady-state behaviour, which is determined by the stationary density. Computing it requires solving a (nonsymmetric) M-matrix algebraic Riccati equation, and indeed computing the stationary density is the most important application of this class of equations. Xue et al. (Numer Math 120:671–700, 2012) provided a componentwise first-order perturbation analysis of this equation, proving that the solution can be computed to high relative accuracy even in the smallest entries, and suggested several algorithms for computing it. An important step in all proposed algorithms is using so-called triplet representations, which are special representations for M-matrices that allow for a high-accuracy variant of Gaussian elimination, the GTH-like algorithm. However, triplet representations for all the M-matrices needed in the algorithm were not found explicitly. This can lead to an accuracy loss that prevents the algorithms from converging in the componentwise sense. In this paper, we focus on the structured doubling algorithm, the most efficient among the proposed methods in Xue et al., and build upon their results, providing (i) explicit and cancellation-free expressions for the needed triplet representations, allowing the algorithm to be performed in a really cancellation-free fashion; (ii) an algorithm to evaluate the final part of the computation to obtain the stationary density; and (iii) a componentwise error analysis for the resulting algorithm, the first explicit one for this class of algorithms. We also present numerical results to illustrate the accuracy advantage of our method over standard (normwise-accurate) algorithms. © 2014, Springer-Verlag Berlin Heidelberg

The Krylov Subspace Methods for the Computation of Matrix Exponentials

The problem of computing the matrix exponential etA arises in many theoretical and practical problems. Many methods have been developed to accurately and efficiently compute this matrix function or its product with a vector, i.e., etAv. In the past few decades, with the increasing need of the computation for large sparse matrices, iterative methods such as the Krylov subspace methods have proved to be a powerful class of methods in dealing with many linear algebra problems. The Krylov subspace methods have been introduced for computing matrix exponentials by Gallopoulos and Saad, and the corresponding error bounds that aim at explaining the convergence properties have been extensively studied. Many of those bounds show that the speed of convergence depends on the norm of the matrix, while some others emphasize the important role played by the spectral distribution. For example, it is shown in a recent work by Ye that the speed of convergence is also determined by the condition number for a symmetric negative definite matrix. Namely the convergence is fast for a well-conditioned matrix no matter how large the norm is. In this dissertation, we derive new error bounds for computing etAv for non-symmetric A, using the spectral information of A. Our result is based on the assumption that A is negative definite, i.e., the field of values of A lies entirely in the left half of the complex plane, such that the underlying dynamic system is stable. The new bounds show that the speed of convergence is related to the size and shape of the rectangle containing the field of values, and they agree with the existing results when A is symmetric. Furthermore, we also derive a simpler error bound for the special case when A is skew-Hermitian. This bound explains an observed convergence behavior where the approximation error initially stagnates for certain number of iterations before it starts to converge. In deriving our new error bounds, we use sharper estimates of the decay property of exponentials of Hessenberg matrices, by constructing Faber polynomial approximating exponential function in the region containing the field of values. The Jacobi elliptic functions are used to construct the conformal mappings and generate the Faber polynomials. We also present numerical tests to demonstrate the behavior of the new error bounds

Approximations for time-dependent distributions in Markovian fluid models

In this paper we study the distribution of the level at time $\theta$ of Markovian fluid queues and Markovian continuous time random walks, the maximum (and minimum) level over $[0,\theta]$, and their joint distributions. We approximate $\theta$ by a random variable $T$ with Erlang distribution and we use an alternative way, with respect to the usual Laplace transform approach, to compute the distributions. We present probabilistic interpretation of the equations and provide a numerical illustration