No abstract

Given a symmetric matrix what is the nearest correlation matrix, that is, the nearest symmetric positive semidefinite matrix with unit diagonal? This problem arises in the finance industry, where the correlations are between stocks. For distance measured in two weighted Frobenius norms we characterize the solution using convex analysis. We show how the modified alternating projections method can be used to compute the solution for the more commonly used of the weighted Frobenius norms. In the finance application the original matrix has many zero or negative eigenvalues; we show that for a certain class of weights the nearest correlation matrix has correspondingly many zero eigenvalues and that this fact can be exploited in the computation.

The scalar sign function is defined for z ∈ C lying off the imaginary axis by sign(z) = 1, Re z > 0, −1, Re z < 0. The matrix sign function can be obtained from any of the definitions in Chapter 1. Note that in the case of the Jordan canonical form and interpolating polynomial definitions, the derivatives sign (k) (z) are zero for k ≥ 1. Throughout this chapter, A ∈ C n×n is assumed to have no eigenvalues on the imaginary axis, so that sign(A) is defined. Note that this assumption implies that A is nonsingular. As we noted in Section 2.4, if A = ZJZ −1 is a Jordan canonical form arranged so that J = diag(J 1 , J 2), where the eigenvalues of J 1 ∈ C p×p lie in the open left half-plane and those of J 2 ∈ C q×q lie in the open right half-plane, then

Abstract. A new algorithm is developed for computing e tA B, where A is an n × n matrix and B is n × n 0 with n 0 ≪ n. The algorithm works for any A, its computational cost is dominated by the formation of products of A with n × n 0 matrices, and the only input parameter is a backward error tolerance. The algorithm can return a single matrix e tA B or a sequence e t k A B on an equally spaced grid of points t k . 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 recent analysis of Al-Mohy and Higham [SIAM J. Matrix Anal. Appl. 31 (2009), pp. 970-989], which provides sharp truncation error bounds expressed in terms of the quantities A k 1/k for a few values of k, where the norms are estimated using a matrix norm estimator. Shifting and balancing are used as preprocessing steps to reduce the cost of the algorithm. Numerical experiments show that the algorithm performs in a numerically stable fashion across a wide range of problems, and analysis of rounding errors and of the conditioning of the problem provides theoretical support. Experimental comparisons with two Krylov-based MATLAB codes show the new algorithm to be sometimes much superior in terms of computational cost and accuracy. An important application of the algorithm is to exponential integrators for ordinary differential equations. It is shown that the sums of the form p k=0 ϕ k (A)u k that arise in exponential integrators, where the ϕ 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, and the algorithm of this paper can therefore be employed.

No abstract

No abstract

Abstract. The scaling and squaring method for the matrix exponential is based on the approximation e A ≈ (rm(2 −s A)) 2 s , where rm(x) is the [m/m] Padé approximant to e x and the integers m and s are to be chosen. Several authors have identified a weakness of existing scaling and squaring algorithms termed overscaling, in which a value of s much larger than necessary is chosen, causing a loss of accuracy in floating point arithmetic. Building on the scaling and squaring algorithm of Higham [SIAM J. Matrix Anal. Appl., 26 (2005), pp. 1179-1193], which is used by the MATLAB function expm, we derive a new algorithm that alleviates the overscaling problem. Two key ideas are employed. The first, specific to triangular matrices, is to compute the diagonal elements in the squaring phase as exponentials instead of from powers of rm. The second idea is to base the backward error analysis that underlies the algorithm on members of the sequence { A k 1/k } instead of A , since for nonnormal matrices it is possible that A k 1/k is much smaller than A , and indeed this is likely when overscaling occurs in existing algorithms. The terms A k 1/k are estimated without computing powers of A by using a matrix 1-norm estimator in conjunction with a bound of the form A k 1/k ≤ max A p 1/p , A q 1/q that holds for certain fixed p and q less than k. The improvements to the truncation error bounds have to be balanced by the potential for a large A to cause inaccurate evaluation of rm in floating point arithmetic. We employ rigorous error bounds along with some heuristics to ensure that rounding errors are kept under control. Our numerical experiments show that the new algorithm generally provides accuracy at least as good as the existing algorithm of Higham at no higher cost, while for matrices that are triangular or cause overscaling it usually yields significant improvements in accuracy, cost, or both.

Abstract. A quadratically convergent Newton method for computing the polar decomposition of a full-rank matrix is presented and analysed. Acceleration parameters are introduced so as to enhance the initial rate of convergence and it is shown how reliable estimates of the optimal parameters may be computed in practice.To add to the known best approximation property of the unitary polar factor, the Hermitian polar factor H of a nonsingular Hermitian matrix A is shown to be a good positive definite approximation to A and 1/2(A / H) is shown to be a best Hermitian positive semi-definite approximation to A. Perturbation bounds for the polar factors are derived.Applications of the polar decomposition to factor analysis, aerospace computations and optimisation are outlined; and a new method is derived for computing the square root of a symmetric positive definite matrix.