A hyperbolic tangent identity and the geometry of Padé sign function iterations

The rational iterations obtained from certain Padé approximations associated with computing the matrix sign function are shown to be equivalent to iterations involving the hyperbolic tangent and its inverse. Using this equivalent formulation many results about these Padé iterations, such as global convergence, the semi-group property under composition, and explicit partial fraction decompositions can be obtained easily. In the second part of the paper it is shown that the behavior of points under the Padé iterations can be expressed, via the Cayley transform, as the combined result of a completely regular iteration and a chaotic iteration. These two iterations are decoupled, with the chaotic iteration taking the form of a multiplicative linear congruential random number generator where the multiplier is equal to the order of the Padé approximation.This research was supported in part by the National Science Foundation under Grant No. ECS-9120643, the Air Force Office of Scientific Research under Grant no. F49620-94-1-0104DEF, and the Office of Naval Research under Grant No. N00014-92-J-1706.     

Error analysis of Padé iterations for computing matrix invariant subspaces

The method of Padé matrix iteration is commonly used for computing matrix sign function and invariant subspaces of a real or complex matrix. In this paper, a detailed rounding error analysis is given for two classical schemes of the Pad’e matrix iteration, using basic matrix floating point arithmetics. Error estimations of computing invariant subspaces by the Padé sign iteration are also provided. Numerical experiments are given to show the numerical behaviors of the Padé iterations and the corresponding subspace computation.      

An algorithm for computing a Padé approximant with minimal degree denominator

In this paper, a new definition of a reduced Padé approximant and an algorithm for its computation are proposed. Our approach is based on the investigation of the kernel structure of the Toeplitz matrix. It is shown that the reduced Padé approximant always has nice properties which the classical Padé approximant possesses only in the normal case. The new algorithm allows us to avoid the appearance of Froissart doublets induced by computer roundoff in the non-normal Padé table.     

Computing the Matrix Cosine

An algorithm is developed for computing the matrix cosine, building on a proposal of Serbin and Blalock. The algorithm scales the matrix by a power of 2 to make the -norm less than or equal to 1, evaluates a Padé approximant, and then uses the double angle formula cos(2A)=2cos(A)²–I to recover the cosine of the original matrix. In addition, argument reduction and balancing is used initially to decrease the norm. We give truncation and rounding error analyses to show that an [8,8] Padé approximant produces the cosine of the scaled matrix correct to machine accuracy in IEEE double precision arithmetic, and we show that this Padé approximant can be more efficiently evaluated than a corresponding Taylor series approximation. We also provide error analysis to bound the propagation of errors in the double angle recurrence. Numerical experiments show that our algorithm is competitive in accuracy with the Schur–Parlett method of Davies and Higham, which is designed for general matrix functions, and it is substantially less expensive than that method for matrices of -norm of order 1. The dominant computational kernels in the algorithm are matrix multiplication and solution of a linear system with multiple right-hand sides, so the algorithm is well suited to modern computer architectures.     

一类变形的牛顿法求解矩阵平方根

提出一些改进的方法来计算矩阵A的平方根,也就是应用一些牛顿法的变形来解决二次矩阵方程．研究表明,改进的方法比牛顿算法和一些已有的牛顿算法的变形效果要好．通过迭代方法,举出一些数值例子说明改进的方法的性能．     

Optimized look-ahead recurrences for adjacent rows in the Padé table

A new look-ahead algorithm for recursively computing Padé approximants is introduced. It generates a subsequence of the Padé approximants on two adjacent rows (defined by fixed numerator degree) of the Padé table. Its two basic versions reduce to the classical Levinson and Schur algorithms if no look-ahead is required. The new algorithm can be viewed as a combination of the look-ahead sawtooth and the look-ahead Levinson and Schur algorithms that we proposed before, but now the look-ahead step size is minimal (as in the sawtooth version) and the computational costs are as low as in the least expensive competing algorithms (including our look-ahead Levinson and Schur algorithms). The underlying recurrences link well-conditioned basic pairs,i.e., pairs of sufficiently different neighboring Padé forms.The algorithm can be used to solve Toeplitz systems of equationsTx = b. In this application it comes in several versions: anO(N ²) Levinson-type form, anO(N ²) Schur-type form, and a superfastO(N log² N) Schur-type version. As an option of the first two versions, the corresponding block LDU decompositions ofT ^–1 orT, respectively, can be found.     

A Padé family of iterations for the matrix sign function and related problems

In this paper we consider the Pad'e family of iterations for computing the matrix sign function and the Padé family of iterations for computing the matrix p‐sector function. We prove that all the iterations of the Padé family for the matrix sign function have a common convergence region. It completes a similar result of Kenney and Laub for half of the Padé family. We show that the iterations of the Padé family for the matrix p‐sector function are well defined in an analogous common region, depending on p. For this purpose we proved that the Padé approximants to the function (1?z)^?σ, 0<σ<1, are a quotient of hypergeometric functions whose poles we have localized. Furthermore we proved that the coefficients of the power expansion of a certain analytic function form a positive sequence and in a special case this sequence has the log‐concavity property. Copyright © 2011 John Wiley & Sons, Ltd.     

A note on Radau and Lobatto formulae for O.D.E.: s

There exist Runge-Kutta methods based on Radau and Lobatto quadrature formulae. One class gives the set of all first and second above diagonal Padé approximations and another class gives the set of all first and second subdiagonal Padé approximations to the exponential function. A new short proof of the strongA-stability of the latter class of methods and a connection between these two classes are presented.     

Convergence Rates of Padé and Padé-Type Approximants

A comparison is made between Padé and Padé-type approximants. LetQ_nbe thenth orthonormal polynomial with respect to a positive measureμwith compact support inC. We show that for functions of the form[formula]wherewis an analytic function on the support ofμ, Padé-type approximants with denominatorQ_ngive a successful and, in general, better approximation procedure than Padé approximation.     

Defects and the Convergence of Padé Approximants

A selective survey is given of convergence results for sequences of Padé approximants. Various approaches for dealing with the convergence problems due to `defects" are discussed. Attention is drawn to the close relationship between analyticity properties of a function and the `smoothness" of its Taylor series coefficients. A new theorem on the convergence of horizontal sequences of Padé approximants to functions in the Baker–Gammel–Wills conjecture function class is presented.     

Algorithms for the matrix <Emphasis Type="Italic">p</Emphasis>th root

New theoretical results are presented about the principal matrix pth root. In particular, we show that the pth root is related to the matrix sign function and to the Wiener–Hopf factorization, and that it can be expressed as an integral over the unit circle. These results are used in the design and analysis of several new algorithms for the numerical computation of the pth root. We also analyze the convergence and numerical stability properties of Newtons method for the inverse pth root. Preliminary computational experiments are presented to compare the methods. AMS subject classification 15A24, 65H10, 65F30Numerical Analysis Report 454, Manchester Centre for Computational Mathematics, July 2004.Dario A. Bini: This work was supported by MIUR, grant number 2002014121.Nicholas J. Higham: This work was supported by Engineering and Physical Sciences Research Council grant GR/R22612 and by a Royal Society – Wolfson Research Merit Award.     

Generalized Padé approximations to the exponential function

The stability properties of the Padé rational approximations to the exponential function are of importance in determining the linear stability properties of several classes of Runge-Kutta methods. It is well known that the Padé approximationR _n,m (z) =N _n,m (z)/M _n,m (z), whereN _n,m (z) is of degreen andM _n,m (z) is of degreem, is A-stable if and only if 0 m – n 2, a result first conjectured by Ehle. In the study of the linear stability properties of the broader class of general linear methods one must generalize these rational approximations. In this paper we introduce a generalization of the Padé approximations to the exponential function and present a method of constructing these approximations for arbitrary order and degree. A generalization of the Ehle inequality is considered and, in the case of the quadratic Padé approximations, evidence is presented that suggests the inequality is both necessary and sufficient for A-stability. However, in the case of the cubic Padé approximations, the inequality is shown to be insufficient for A-stability. A generalization of the restricted Padé approximation, in which the denominator has a singlem-fold zero, is also introduced. A procedure for the construction of these restricted approximations is described, and results are presented on the A-stability of the restricted quadratic Padé approximations. Finally, to demonstrate the connection between a generalized Padé approximation and a general linear method, a specific general linear method is constructed with a stability region corresponding to a given quadratic Padé approximation.     

Efficient sixth order methods for nonlinear oscillation problems

A class of two-step (hybrid) methods is considered for solving pure oscillation second order initial value problems. The nonlinear system, which results on applying methods of this type to a nonlinear differential system, may be solved using a modified Newton iteration scheme. From this class the author has derived methods which are fourth order accurate,P-stable, require only two (new) function evaluations per iteration and have a true real perfect square iteration matrix. Now, we propose an extension to sixth order,P-stable methods which require only three (new) function evaluations per iteration and for which the iteration matrix is a true realperfect cube. This implies that at most one real matrix must be factorised at each step. These methods have been implemented in a new variable step, local error controlling code.     

A superfast algorithm for multi-dimensional Padé systems

For a vector ofk+1 matrix power series, a superfast algorithm is given for the computation of multi-dimensional Padé systems. The algorithm provides a method for obtaining matrix Padé, matrix Hermite Padé and matrix simultaneous Padé approximants. When the matrix power series is normal or perfect, the algorithm is shown to calculate multi-dimensional matrix Padé systems of type (n ₀,...,n _k ) inO(n · log²n) block-matrix operations, where n=n ₀+...+n _k . Whenk=1 and the power series is scalar, this is the same complexity as that of other superfast algorithms for computing Padé systems. Whenk>1, the fastest methods presently compute these matrix Padé approximants with a complexity ofO(n²). The algorithm succeeds also in the non-normal and non-perfect case, but with a possibility of an increase in the cost complexity.Supported in part by NSERC grant No. A8035.Partially supported by NSERC operating grant No. 6194.     

A primal—dual potential reduction method for problems involving matrix inequalities

We describe a potential reduction method for convex optimization problems involving matrix inequalities. The method is based on the theory developed by Nesterov and Nemirovsky and generalizes Gonzaga and Todd's method for linear programming. A worst-case analysis shows that the number of iterations grows as the square root of the problem size, but in practice it appears to grow more slowly. As in other interior-point methods the overall computational effort is therefore dominated by the least-squares system that must be solved in each iteration. A type of conjugate-gradient algorithm can be used for this purpose, which results in important savings for two reasons. First, it allows us to take advantage of the special structure the problems often have (e.g., Lyapunov or algebraic Riccati inequalities). Second, we show that the polynomial bound on the number of iterations remains valid even if the conjugate-gradient algorithm is not run until completion, which in practice can greatly reduce the computational effort per iteration.We describe in detail how the algorithm works for optimization problems withL Lyapunov inequalities, each of sizem. We prove an overallworst-case operation count of O(m ^5.5L^1.5). Theaverage-case complexity appears to be closer to O(m ⁴L^1.5). This estimate is justified by extensive numerical experimentation, and is consistent with other researchers' experience with the practical performance of interior-point algorithms for linear programming.This result means that the computational cost of extending current control theory based on the solution of Lyapunov or Riccatiequations to a theory that is based on the solution of (multiple, coupled) Lyapunov or Riccatiinequalities is modest.Supported by the Belgian National Fund for Scientific Research (NFWO). Research supported in part by the Belgian program on Interuniversity Attraction Poles (IUAP 17 and 50) initiated by the Belgian State, Prime Minister's Office, Science Policy Programming.Research supported in part by AFOSR (under F49620-92-J-0013), NSF (under ECS-9222391) and ARPA (under F49620-93-1-0085).     

The multipoint Padé table and general recurrences for rational interpolation

We first review briefly the Newton-Padé approximation problem and the analogous problem with additional interpolation conditions at infinity, which we call multipoint Padé approximation problem. General recurrence formulas for the Newton-Padé table combine either two pairs of Newton-Padé forms or one such pair and a pair of multipoint Padé forms. We show that, likewise, certain general recurrences for the multipoint Padé table compose two pairs of multipoint Padé forms to get a new pair of multipoint Padé forms. We also discuss the possibility of superfast, i.e.,O(n log² n) algorithms for certain rational interpolation problems.     

Enriched Methods for Large-Scale Unconstrained Optimization

This paper describes a class of optimization methods that interlace iterations of the limited memory BFGS method (L-BFGS) and a Hessian-free Newton method (HFN) in such a way that the information collected by one type of iteration improves the performance of the other. Curvature information about the objective function is stored in the form of a limited memory matrix, and plays the dual role of preconditioning the inner conjugate gradient iteration in the HFN method and of providing an initial matrix for L-BFGS iterations. The lengths of the L-BFGS and HFN cycles are adjusted dynamically during the course of the optimization. Numerical experiments indicate that the new algorithms are both effective and not sensitive to the choice of parameters.     

A type of matrix Padé approximant inspired by scalar component models

In this paper we define a type of matrix Padé approximant inspired by the identification stage of multivariate time series models considering scalar component models. Of course, the formalization of certain properties in the matrix Padé approximation framework can be applied to time series models and in other fields. Specifically, we want to study matrix Padé approximants as follows: to find rational representations (or rational approximations) of a matrix formal power series, with both matrix polynomials, numerator and denominator, satisfying three conditions: (a) minimum row degrees for the numerator and denominator, (b) an invertible denominator at the origin, and (c) canonical representation (without free parameters).     

Smoothing schemes for reaction-diffusion systems with nonsmooth data

Cox and Matthews [S.M. Cox, P.C. Matthews, Exponential time differencing for stiff systems, J. Comput. Phys. 176 (2002) 430–455] developed a class of Exponential Time Differencing Runge–Kutta schemes (ETDRK) for nonlinear parabolic equations; Kassam and Trefethen [A.K. Kassam, Ll. N. Trefethen, Fourth-order time stepping for stiff pdes, SIAM J. Sci. Comput. 26 (2005) 1214–1233] have shown that these schemes can suffer from numerical instability and they proposed a modified form of the fourth-order (ETDRK4) scheme. They use complex contour integration to implement these schemes in a way that avoids inaccuracies when inverting matrix polynomials, but this approach creates new difficulties in choosing and evaluating the contour for larger problems. Neither treatment addresses problems with nonsmooth data, where spurious oscillations can swamp the numerical approximations if one does not treat the problem carefully. Such problems with irregular initial data or mismatched initial and boundary conditions are important in various applications, including computational chemistry and financial engineering. We introduce a new version of the fourth-order Cox–Matthews, Kassam–Trefethen ETDRK4 scheme designed to eliminate the remaining computational difficulties. This new scheme utilizes an exponential time differencing Runge–Kutta ETDRK scheme using a diagonal Padé approximation of matrix exponential functions, while to deal with the problem of nonsmooth data we use several steps of an ETDRK scheme using a sub-diagonal Padé formula. The new algorithm improves computational efficiency with respect to evaluation of the high degree polynomial functions of matrices, having an advantage of splitting the matrix polynomial inversion problem into a sum of linear problems that can be solved in parallel. In this approach it is only required that several backward Euler linear problems be solved, in serial or parallel. Numerical experiments are described to support the new scheme.