Computing a logarithm of a unitary matrix with general spectrum

We analyze an algorithm for computing a skew-Hermitian logarithm of a unitary matrix. This algorithm is very easy to implement using standard software and it works well even for unitary matrices with no spectral conditions assumed. Certain examples, with many eigenvalues near -1, lead to very non-Hermitian output for other basic methods of calculating matrix logarithms. Altering the output of these algorithms to force an Hermitian output creates accuracy issues which are avoided in the considered algorithm. A modification is introduced to deal properly with the $J$-skew symmetric unitary matrices. Applications to numerical studies of topological insulators in two symmetry classes are discussed.Comment: Added discussion of Floquet Hamiltonian

Computing the Wave-Kernel Matrix Functions

We derive an algorithm for computing the wave-kernel functions cosh \surd A and sinhc\surd A for an arbitrary square matrix A, where sinhcz = sinh(z)/z. The algorithm is based on Pad\'e approximation and the use of double angle formulas. We show that the backward error of any approximation to cosh \surd A can be explicitly expressed in terms of a hypergeometric function. To bound the backward error we derive and exploit a new bound for \| Ak\| 1/k that is sharper than one previously obtained by Al-Mohy and Higham [SIAM J. Matrix Anal. Appl., 31 (2009), pp. 970-- 989]. The amount of scaling and the degree of the Pad\'e approximant are chosen to minimize the computational cost subject to achieving backward stability for cosh \surd A in exact arithmetic. Numerical experiments show that the algorithm behaves in a forward stable manner in floating-point arithmetic and is superior in this respect to the general purpose Schur--Parlett algorithm applied to these functions

Computing the square roots of matrices with central symmetry

For computing square roots of a nonsingular matrix A, which are functions of A, two well known fast and stable algorithms, which are based on the Schur decomposition of A, were proposed by Bj¨ork and Hammarling [3], for square roots of general complex matrices, and by Higham [10], for real square roots of real matrices. In this paper we further consider (the computation of) the square roots of matrices with central symmetry. We first investigate the structure of the square roots of these matrices and then develop several algorithms for computing the square roots. We show that our algorithms ensure significant savings in computational costs as compared to the use of standard algorithms for arbitrary matrices.Fundação para a Ciência e a Tecnologia (FCT

Computing the matrix Mittag-Leffler function with applications to fractional calculus

The computation of the Mittag-Leffler (ML) function with matrix arguments, and some applications in fractional calculus, are discussed. In general the evaluation of a scalar function in matrix arguments may require the computation of derivatives of possible high order depending on the matrix spectrum. Regarding the ML function, the numerical computation of its derivatives of arbitrary order is a completely unexplored topic; in this paper we address this issue and three different methods are tailored and investigated. The methods are combined together with an original derivatives balancing technique in order to devise an algorithm capable of providing high accuracy. The conditioning of the evaluation of matrix ML functions is also studied. The numerical experiments presented in the paper show that the proposed algorithm provides high accuracy, very often close to the machine precision

Matrix Scaling and Balancing via Box Constrained Newton's Method and Interior Point Methods

In this paper, we study matrix scaling and balancing, which are fundamental problems in scientific computing, with a long line of work on them that dates back to the 1960s. We provide algorithms for both these problems that, ignoring logarithmic factors involving the dimension of the input matrix and the size of its entries, both run in time $\widetilde{O}\left(m\log \kappa \log^2 (1/\epsilon)\right)$ where $\epsilon$ is the amount of error we are willing to tolerate. Here, $\kappa$ represents the ratio between the largest and the smallest entries of the optimal scalings. This implies that our algorithms run in nearly-linear time whenever $\kappa$ is quasi-polynomial, which includes, in particular, the case of strictly positive matrices. We complement our results by providing a separate algorithm that uses an interior-point method and runs in time $\widetilde{O}(m^{3/2} \log (1/\epsilon))$. In order to establish these results, we develop a new second-order optimization framework that enables us to treat both problems in a unified and principled manner. This framework identifies a certain generalization of linear system solving that we can use to efficiently minimize a broad class of functions, which we call second-order robust. We then show that in the context of the specific functions capturing matrix scaling and balancing, we can leverage and generalize the work on Laplacian system solving to make the algorithms obtained via this framework very efficient.Comment: To appear in FOCS 201

Analytical solutions for compartmental models of contaminant transport in enclosed spaces

Understanding the transport of hazardous airborne materials within buildings and other enclosed spaces is important for predicting and mitigating the impacts of deliberate terrorist releases of chemical and biological materials. Multizone models provide an approach to modelling the contamination levels in enclosed spaces but in certain cases they can be computationally expensive. Alternative methods are being explored at DSTL that involve the direct solution of the contaminant dynamics equation. The Study Group was asked to identify the limits to this alternative approach and to explore its extension