A New sqrtm for Matlab

Байду номын сангаас
A New sqrtm for Matlab
Nicholas J. Higham January 4, 1999
Abstract
information about the stability and conditioning of the problem. number, stability
Matlab's function sqrtm computes a square root of a matrix. We propose a replacement for the sqrtm in Matlab 5.2 that is more accurate and returns useful
2 Matrix Square Roots
We begin by summarizing some salient features of the theory of matrix square roots. Further details can be found in 4], 6, Sec. 6.4]. The matrix X 2 C is a square root of A 2 C if X 2 = A. Any nonsingular matrix has a square root, but whether a singular matrix has a square root depends on the Jordan structure of the zero eigenvalues. If A
1 Introduction
root of a matrix. The function works by reducing the matrix to Schur form and then applying a recurrence of Parlett for computing a general function of a triangular matrix. An error estimate is computed and if it is too large then an attempt is made to improve the accuracy of the computed square root. The function sqrtm in versions of Matlab up to version 5.2 can be much less accurate than is warranted by the condition of the problem. We propose a replacement for sqrtm that is more accurate and returns useful information about the stability and conditioning of the problem. In Sections 2 and 3 we present some background theory on the existence of matrix square roots and their stability and conditioning. In Section 4 we describe the existing function sqrtm and then in Section 5 the proposed replacement. Numerical experiments in Section 6 compare the new and the old routines and conclusions are given in Section 7.
1
is nonsingular and has s distinct eigenvalues then it has precisely 2 square roots that are expressible as polynomials in the matrix A; if some eigenvalue appears in more than one Jordan block then there are in nitely many additional square roots, none of which is expressible as a polynomial in A. To illustrate, the singular matrix 0 1 (2.1) 0 0 has no square root, a diagonal matrix diag(d ) 2 C with distinct diagonal entries has p exactly 2 square roots, diag( d ), and the n n identity matrix has in nitely many square roots for n > 2, including any Householder matrix. Although a square root is never unique when it exists, there is a distinguished square root of particular interest: the one all of whose eigenvalues lie in the right half-plane. To make this square root uniquely de ned, we (arbitrarily) map any eigenvalues on the negative real axis to the positive imaginary axis. This square root is called the principal square root and it is a polynomial in the original matrix. When A is symmetric positive de nite the principal square root is the unique symmetric positive de nite square root, and when A is real and has a real square root that is a polynomial in A, the principal square root is real.
2
F
2
F
F
F
F
F
2
where F 0(X ) : C ! C is the Frechet derivative of F (X ) = X 2 ? A at X . Hence X = F 0(X )?1( A ? X 2 ) and taking norms and solving the resulting quadratic inequality in k X k gives the sharp inequality
Key words. matrix square root, Matlab, Schur decomposition, condition AMS subject classi cations. 65F30
Matlab has included since at least version 3 a function sqrtm for computing a square
s n n i n i
3 Stability and Conditioning of Square Roots
b Two measures of the accuracy of a computed square root X of A are the relative residual b 2 k =kAk and the relative error kX ? X k =kX k , where X is the square root of b kA ? X P interest. Here, we are using the Frobenius norm, kAk = ( ja j2)1 2 . To know how small the relative residual can be expected to be we consider the correctly rounded exact e square root: X = fl(X ) = X + E , where kE k ukX k , with u the unit roundo . We have e kA ? X 2k = kXE + EX + E 2k2 2kX k kE k + kE k2 (2u + u2)kX k2 (2u + u2) (X )kAk ; (3.1) where the stability factor
ISSN 1360-1725
UMIST
A New sqrtm for Matlab N. J. Higham Numerical Analysis Report No. 336 January 1999
Manchester Centre for Computational Mathematics Numerical Analysis Reports
F F F F F i;j ij = F F F F F F F F
k kX (X ) = kXk = kX 2k (3.2) kA k 1: b Thus the best we can expect is that the relative residual of a computed X is of order b (X )u. The stability factor (X ) can be arbitrarily large and it is related to the condition number with respect to inversion, (X ) = kX k kX ?1k , by (X ) (X ) (X ): (A) b The relative error of X depends on the conditioning of X , so we now derive a condition number, following 4]. If A = X 2 and A + A = (X + X )2 then X X + XX = A ? X 2 , which can be written F 0(X ) X = A ? X 2;