The growing interest in computing structured matrix functions stems from the fact that preserving and exploiting the structure of matrices can help us gain physically meaningful solutions with less computational cost and memory requirement. The work presented here is divided into two parts. The first part deals with the computation of functions of structured matrices. The second part is concerned with the structured error analysis in the computation of matrix functions. We present algorithms applying the inverse scaling and squaring method and using the Schur-like form of the symplectic matrices as an alternative to the algorithms using the Schur decomposition to compute the logarithm of symplectic matrices. There are two main calculations in the inverse scaling and squaring method: taking a square root and evaluating the Padé approximants. Numerical experiments suggest that using the Schur-like form with the structure preserving iterations for the square root helps us to exploit the Hamiltonian structure of the logarithm of symplectic matrices. Some type of matrices are nearly structured. We discuss the conditions for using the nearest structured matrix to the nearly structured one by analysing the forward error bounds. Since the structure preserving algorithms for computing the functions of matrices provide advantages in terms of accuracy and data storage we suggest to compute the function of the nearest structured matrix. The analysis is applied to the nearly unitary, nearly Hermitian and nearly positive semi-definite matrices for the matrix logarithm, square root, exponential, cosine and sine functions. It is significant to investigate the effect of the structured perturbations in the sensitivity analysis of matrix functions. We study the structured condition number of matrix functions defined between smooth square matrix manifolds. We develop algorithms computing and estimating the structured condition number. We also present the lower and upper bounds on the structured condition number, which are cheaper to compute than the "exact" structured condition number. We observe that the lower bounds give a good estimation for the structured condition numbers. Comparing the structured and unstructured condition number reveals that they can differ by several orders of magnitude. Having discussed how to compute the structured condition number of matrix functions defined between smooth square matrix manifolds we apply the theory of structured condition numbers to the structured matrix factorizations. We measure the sensitivity of matrix factors to the structured perturbations for the structured polar decomposition, structured sign factorization and the generalized polar decomposition. Finally, we consider the unstructured perturbation analysis for the canonical generalized polar decomposition by using three different methods. Apart from theoretical aspect of the perturbation analysis, perturbation bounds obtained from these methods are compared numerically and our findings show an improvement on the sharpness of the perturbation bounds in the literature.

In this thesis, we consider real matrix functions that depend on two parameters and study the problem of how to detect and approximate parameters' values where the singular values coalesce. We prove several results connecting the existence of coalescing points to the periodic structure of the smooth singular values decomposition computed around the boundary of a domain enclosing the points. This is further used to develop algorithms for the detection and approximation of coalescing points in planar regions. Finally, we present techniques for continuing curves of coalescing singular values of matrices depending on three parameters, and illustrate how these techniques can be used to locate coalescing singular values of complex-valued matrices depending on three parameters.

In Markov chain models in finance and healthcare a transition matrix over a certain time interval is needed but only a transition matrix over a longer time interval may be available. The problem arises of determining a stochastic $p$th root of astochastic matrix (the given transition matrix). By exploiting the theory of functions of matrices, we develop results on the existence and characterization of stochastic $p$th roots. Our contributions include characterization of when a real matrix hasa real $p$th root, a classification of $p$th roots of a possibly singular matrix,a sufficient condition for a $p$th root of a stochastic matrix to have unit row sums,and the identification of two classes of stochastic matrices that have stochastic $p$th roots for all $p$. We also delineate a wide variety of possible configurationsas regards existence, nature (primary or nonprimary), and number of stochastic roots,and develop a necessary condition for existence of a stochastic root in terms of the spectrum of the given matrix. On the computational side, we emphasize finding an approximate stochastic root: perturb the principal root $A^{1/p}$ or the principal logarithm $\log(A)$ to the nearest stochastic matrix or the nearest intensity matrix, respectively, if they are not valid ones;minimize the residual $\normF{X^p-A}$ over all stochastic matrices $X$ and also over stochastic matrices that are primary functions of $A$. For the first two nearness problems, the global minimizers are found in the Frobenius norm. For the last two nonlinear programming problems, we derive explicit formulae for the gradient and Hessian of the objective function $\normF{X^p-A}^2$ and investigate Newton's method, a spectral projected gradient method (SPGM) and the sequential quadratic programming method to solve the problem as well as various matrices to start the iteration. Numerical experiments show that SPGM starting with the perturbed $A^{1/p}$to minimize $\normF{X^p-A}$ over all stochastic matrices is method of choice.Finally, a new algorithm is developed for computing arbitrary real powers $A^\a$ of a matrix $A\in\mathbb{C}^{n\times n}$. The algorithm starts with a Schur decomposition,takes $k$ square roots of the triangular factor $T$, evaluates an $[m/m]$ Pad\'e approximant of $(1-x)^\a$ at $I - T^$, and squares the result $k$ times. The parameters $k$ and $m$ are chosen to minimize the cost subject to achieving double precision accuracy in the evaluation of the Pad\'e approximant, making use of a result that bounds the error in the matrix Pad\'e approximant by the error in the scalar Pad\'e approximant with argument the norm of the matrix. The Pad\'e approximant is evaluated from the continued fraction representation in bottom-up fashion, which is shown to be numerically stable. In the squaring phase the diagonal and first superdiagonal are computed from explicit formulae for $T^$, yielding increased accuracy. Since the basic algorithm is designed for $\a\in(-1,1)$, a criterion for reducing an arbitrary real $\a$ to this range is developed, making use of bounds for the condition number of the $A^\a$ problem. How best to compute $A^k$ for a negative integer $k$ is also investigated. In numerical experiments the new algorithm is found to be superior in accuracy and stability to several alternatives,including the use of an eigendecomposition, a method based on the Schur--Parlett\alg\ with our new algorithm applied to the diagonal blocks and approaches based on the formula $A^\a = \exp(\a\log(A))$.

No description available.

New algorithms for the matrix exponential and its Fréchet derivative are presented. First, we derive a new scaling and squaring algorithm (denoted expm[new]) for computing eA, where A is any square matrix, that mitigates the overscaling problem. The algorithm is built on the algorithm of Higham [SIAM J.Matrix Anal. Appl., 26(4): 1179-1193, 2005] but improves on it by two key features. The first, specific to triangular matrices, is to compute the diagonal elements in the squaring phase as exponentials instead of powering them. The second is to base the backward error analysis that underlies the algorithm on members of the sequence {||Ak||1/k} instead of ||A||. The terms ||Ak||1/k are estimated without computing powers of A by using a matrix 1-norm estimator. Second, a new algorithm is developed for computing the action of the matrix exponential on a matrix, etAB, where A is an n x n matrix and B is n x n₀ with n₀ << n. The algorithm works for any A, its computational cost is dominated by the formation of products of A with n x n₀ matrices, and the only input parameter is a backward error tolerance. The algorithm can return a single matrix etAB or a sequence etkAB on an equally spaced grid of points tk. It uses the scaling part of the scaling and squaring method together with a truncated Taylor series approximation to the exponential. It determines the amount of scaling and the Taylor degree using the strategy of expm[new].Preprocessing steps are used to reduce the cost of the algorithm. An important application of the algorithm is to exponential integrators for ordinary differential equations. It is shown that the sums of the form $\sum_{k=0}^p\varphi_k(A)u_k$ that arise in exponential integrators, where the $\varphi_k$ are related to the exponential function, can be expressed in terms of a single exponential of a matrix of dimension $n+p$ built by augmenting $A$ with additional rows and columns. Third, a general framework for simultaneously computing a matrix function, $f(A)$, and its Fréchet derivative in the direction $E$, $L_f(A,E)$, is established for a wide range of matrix functions. In particular, we extend the algorithm of Higham and $\mathrm{expm_{new}}$ to two algorithms that intertwine the evaluation of both $e^A$ and $L(A,E)$ at a cost about three times that for computing $e^A$ alone. These two extended algorithms are then adapted to algorithms that simultaneously calculate $e^A$ together with an estimate of its condition number. Finally, we show that $L_f(A,E)$, where $f$ is a real-valued matrix function and $A$ and $E$ are real matrices, can be approximated by $\Im f(A+ihE)/h$ for some suitably small $h$. This approximation generalizes the complex step approximation known in the scalar case, and is proved to be of second order in $h$ for analytic functions $f$ and also for the matrix sign function. It is shown that it does not suffer the inherent cancellation that limits the accuracy of finite difference approximations in floating point arithmetic. However, cancellation does nevertheless vitiate the approximation when the underlying method for evaluating $f$ employs complex arithmetic. The complex step approximation is attractive when specialized methods for evaluating the Fréchet derivative are not available.

Many problems in applied mathematics and physics are formulated most naturally in terms of matrices, and can be solved by computing functions of these matrices. For example, in quantum mechanics, the coherent dynamics of physical systems is described by the matrix exponential of their Hamiltonian. In state of the art experiments, one can now observe such unitary evolution of many-body systems, which is of fundamental interest in the study of many-body quantum phenomena. On the other hand the theoretical simulation of such non-equilibrium many-body dynamics is very challenging. In this thesis, we develop a symbolic approach to matrix functions and quantum dynamics based on a novel algebraic structure we identify for sets of walks on graphs. We begin by establishing the graph theoretic equivalent to the fundamental theorem of arithmetic: all the walks on any finite digraph uniquely factorise into products of prime elements. These are the simple paths and simple cycles, walks forbidden from visiting any vertex more than once. We give an algorithm that efficiently factorises individual walks and obtain a recursive formula to factorise sets of walks. This yields a universal continued fraction representation for the formal series of all walks on digraphs. It only involves simple paths and simple cycles and is thus called a path-sum. In the second part, we recast matrix functions into path-sums. We present explicit results for a matrix raised to a complex power, the matrix exponential, matrix inverse, and matrix logarithm. We introduce generalised matrix powers which extend desirable properties of the Drazin inverse to all powers of a matrix. In the third part, we derive an intermediary form of path-sum, called walk-sum, relying solely on physical considerations. Walk-sum describes the dynamics of a quantum system as resulting from the coherent superposition of its histories, a discrete analogue to the Feynman path-integrals. Using walk-sum we simulate the dynamics of quantum random walks and of Rydberg-excited Mott insulators. Using path-sum, we demonstrate many-body Anderson localisation in an interacting disordered spin system. We give two observable signatures of this phenomenon: localisation of the system magnetisation and of the linear magnetic response function. Lastly we return to the study of sets of walks. We show that one can construct as many representations of series of walks as there are ways to define a walk product such that the factorisation of a walk always exist and is unique. Illustrating this result we briefly present three further methods to evaluate functions of matrices. Regardless of the method used, we show that graphs are uniquely characterised, up to an isomorphism, by the prime walks they sustain.