Search

© 2014 SIAM.The Fréchet derivative Lf of a matrix function f: ℂn×n → ℂn×n controls the sensitivity of the function to small perturbations in the matrix. While much is known about the properties of Lf and how to compute it, little attention has been given to higher order Fréchet derivatives. We derive sufficient conditions for the kth Fréchet derivative to exist and be continuous in its arguments and we develop algorithms for computing the kth derivative and its Kronecker form. We analyze the level-2 absolute condition number of a matrix function ("the condition number of the condition number") and bound it in terms of the second Fréchet derivative. For normal matrices and the exponential we show that in the 2-norm the level-1 and level-2 absolute condition numbers are equal and that the relative condition numbers are within a small constant factor of each other. We also obtain an exact relationship between the level-1 and level-2 absolute condition numbers for the matrix inverse and arbitrary nonsingular matrices, as well as a weaker connection for Hermitian matrices for a class of functions that includes the logarithm and square root. Finally, the relation between the level-1 and level-2 condition numbers is investigated more generally through numerical experiments.

To approximate $f(A){b}$---the action of a matrix function on a vector---by a Krylov subspace method, restarts may become mandatory due to storage requirements for the Arnoldi basis or due to the growing computational complexity of evaluating $f$ on a Hessenberg matrix of growing size. A number of restarting methods have been proposed in the literature in recent years and there has been substantial algorithmic advancement concerning their stability and computational efficiency. However, the question under which circumstances convergence of these methods can be guaranteed has remained largely unanswered. In this paper we consider the class of Stieltjes functions and a related class, which contain important functions like the (inverse) square root and the matrix logarithm. For these classes of functions we present new theoretical results which prove convergence for Hermitian positive definite matrices $A$ and arbitrary restart lengths. We also propose a modification of the Arnoldi approximation which guaran...

The Schur--Pade algorithm [N. J. Higham and L. Lin, SIAM J. Matrix Anal. Appl., 32 (2011), pp. 1056--1078] computes arbitrary real powers $A^t$ of a matrix $A\in\mathbb{C}^{n\times n}$ using the building blocks of Schur decomposition, matrix square roots, and Pade approximants. We improve the algorithm by basing the underlying error analysis on the quantities $\|(I-A)^k\|^{1/k}$, for several small $k$, instead of $\|I-A\|$. We extend the algorithm so that it computes along with $A^t$ one or more Frechet derivatives, with reuse of information when more than one Frechet derivative is required, as is the case in condition number estimation. We also derive a version of the extended algorithm that works entirely in real arithmetic when the data is real. Our numerical experiments show the new algorithms to be superior in accuracy to, and often faster than, the original Schur--Pade algorithm for computing matrix powers and more accurate than several alternative methods for computing the Frechet derivative. They ...

Newton’s method for the inverse matrix $p$th root, $A^{-1/p}$, has the attraction that it involves only matrix multiplication. We show that if the starting matrix is $c^{-1}I$ for $c\in\mathbb{R}^+$ then the iteration converges quadratically to $A^{-1/p}$ if the eigenvalues of $A$ lie in a wedge-shaped convex set containing the disc $\{z: |z-c^p| < c^p\}$. We derive an optimal choice of $c$ for the case where $A$ has real, positive eigenvalues. An application is described to roots of transition matrices from Markov models, in which for certain problems the convergence condition is satisfied with $c=1$. Although the basic Newton iteration is numerically unstable, a coupled version is stable and a simple modification of it provides a new coupled iteration for the matrix $p$th root. For general matrices we develop a hybrid algorithm that computes a Schur decomposition, takes square roots of the upper (quasi-) triangular factor, and applies the coupled Newton iteration to a matrix for which fast convergence is guaranteed. The new algorithm can be used to compute either $A^{1/p}$ or $A^{-1/p}$, and for large $p$ that are not highly composite it is more efficient than the method of Smith based entirely on the Schur decomposition.

In this paper we explore the computation of the matrix exponential in a manner that is consistent with Lie-group structure. Our point of departure is the method of generalized polar decompositions, which we modify and combine with similarity transformations that bring the underlying matrix to a form more amenable to efficient computation. We develop techniques valid for a range of Lie groups: the orthogonal group, the symplectic group, Lorentz, isotropy, and scaling groups. However, the GPD approach is equally promising in a more general context. Even when Lie-group structure is not at issue, our algorithm is more efficient in many settings than classical methods for the computation of the matrix exponential.

An algorithm for computing matrix functions is presented. It employs a Schur decomposition with reordering and blocking followed by the block form of a recurrence of Parlett, with functions of the nontrivial diagonal blocks evaluated via a Taylor series. A parameter is used to balance the conflicting requirements of producing small diagonal blocks and keeping the separations of the blocks large. The algorithm is intended primarily for functions having a Taylor series with an infinite radius of convergence, but it can be adapted for certain other functions, such as the logarithm. Novel features introduced here include a convergence test that avoids premature termination of the Taylor series evaluation and an algorithm for reordering and blocking the Schur form. Numerical experiments show that the algorithm is competitive with existing special-purpose algorithms for the matrix exponential, logarithm, and cosine. Nevertheless, the algorithm can be numerically unstable with the default choice of its blocking parameter (or in certain cases for all choices), and we explain why determining the optimal parameter appears to be a very difficult problem. A MATLAB implementation is available that is much more reliable than the function funm in MATLAB 6.5 (R13).

The inverse scaling and squaring method for evaluating the logarithm of a matrix takes repeated square roots to bring the matrix close to the identity, computes a Pade approximant, and then scales back. We analyze several methods for evaluating the Pade approximant, including Horner's method (used in some existing codes), suitably customized versions of the Paterson--Stockmeyer method and Van Loan's variant, and methods based on continued fraction and partial fraction expansions. The computational cost, storage, and numerical accuracy of the methods are compared. We find the partial fraction method to be the best method overall and illustrate the benefits it brings to a transformation-free form of the inverse scaling and squaring method recently proposed by Cheng, Higham, Kenney, and Laub [ SIAM J. Matrix Anal. Appl., 22 (2001), pp. 1112--1125]. We comment briefly on how the analysis carries over to the matrix exponential.

The standard inverse scaling and squaring algorithm for computing the matrix logarithm begins by transforming the matrix to Schur triangular form in order to facilitate subsequent matrix square root and Pade approximation computations. A transformation-free form of this method that exploits incomplete Denman--Beavers square root iterations and aims for a specified accuracy (ignoring roundoff) is presented. The error introduced by using approximate square roots is accounted for by a novel splitting lemma for logarithms of matrix products. The number of square root stages and the degree of the final Pade approximation are chosen to minimize the computational work. This new method is attractive for high-performance computation since it uses only the basic building blocks of matrix multiplication, LU factorization and matrix inversion.

In this work we (i) use the theory of piecewise analytic functions to represent the Frechet derivative of any primary matrix function, in particular of primary logarithms; (ii) propose an indicator to assess inherent difficulties to compute a logarithm; and (iii) revisit Pade approximation techniques for the principal logarithm of block triangular matrices.

The Schur--Frechet method of evaluating matrix functions consists of putting the matrix in upper triangular form, computing the scalar function values along the main diagonal, and then using the Frechet derivative of the function to evaluate the upper diagonals. This approach requires a reliable method of computing the Frechet derivative. For the logarithm this can be done by using repeated square roots and a hyperbolic tangent form of the logarithmic Frechet derivative. Pade approximations of the hyperbolic tangent lead to a Schur--Frechet algorithm for the logarithm that avoids problems associated with the standard "inverse scaling and squaring" method. Inverting the order of evaluation in the logarithmic Frechet derivative gives a method of evaluating the derivative of the exponential. The resulting Schur--Frechet algorithm for the exponential gives superior results compared to standard methods on a set of test problems from the literature.

In principle, the exponential of a matrix could be computed in many ways. Methods involving approximation theory, differential equations, the matrix eigenvalues, and the matrix characteristic polyn...

A sensitivity theory based on Frechet derivatives is presented that has both theoretical and computational advantages. Theoretical results such as a generalization of Van Loan’s work on the matrix exponential are easily obtained: matrix functions are least sensitive at normal matrices. Computationally, the central problem is to estimate the norm of the Frechet derivative, since this is equal to the function’s condition number. Two norm-estimation procedures are given; the first is based on a finite-difference approximation of the Frechet derivative and costs only two extra function evaluations. The second method was developed specifically for the exponential and logarithmic functions; it is based on a trapezoidal approximation scheme suggested by the chain rule for the identity $e^X = ( e^{X/2^n } )^{2^n } $. This results in an infinite sequence of coupled Sylvester equations that, when truncated, is uniquely suited to the “scaling and squaring” procedure for $e^X $ or the “inverse scaling and squaring” p...

Matrix power series may slowly converge or even diverge if some eigenvalues of the matrix are near the boundary or outside the disk of convergence. In this case it is proposed to apply suitably chosen summability methods to accelerate or generate convergence; special attention is paid to Euler methods. The matrix logarithm appearing in connection with stationary Markov chains is considered as an example.