Using the Sherman-Morrison-Woodbury inversion formula for a fast solution of tridiagonal block Toeplitz systems

A fast numerical algorithm for solving systems of linear equations with tridiagonal block Toeplitz matrices is presented. The algorithm is based on a preliminary factorization of the generating quadratic matrix polynomial associated with the Toeplitz matrix, followed by the Sherman-Morrison-Woodbury inversion formula and solution of two bidiagonal and one diagonal block Toeplitz systems. Tight estimates of the condition numbers are provided for the matrix system and the main matrix systems generated during the preliminary factorization. The emphasis is put on rigorous stability analysis to rounding errors of the Sherman-Morrison-Woodbury inversion. Numerical experiments are provided to illustrate the theory.     

Using FFT-based techniques in polynomial and matrix computations: recent advances and applicatons

Computations with univariate polynomials, like the evaluation of product, quotient, remainder, greatest common divisor, etc, are closely related to linear algebra computations performed with structured matrices having the Toeplitz-like or the Hankel-like structures.

The discrete Fourier transform, and the FFT algorithms for its computation, constitute a powerful tool for the design and analysis of fast algorithms for solving problems involving polynomials and structured matrices.

We recall the main correlations between polynomial and matrix computations and present two recent results in this field: in particular, we show how Fourier methods can speed up the solution of a wide class of problems arising in queueing theory where certain Markov chains, defined by infinite block Toeplitz matrices in generalized Hessenberg form, must be solved. Moreover, we present a new method for the numerical factorization of polynomials based on a matrix generalization of Koenig's theorem. This method, that relies on the evaluation/interpolation technique at the Fourier points, reduces the problem of polynomial factorization to the computation of the LU decomposition of a banded Toeplitz matrix with its rows and columns suitably permuted. Numerical experiments that show the effectiveness of our algorithms are presented     

The convergence of padé approximants and the size of the power series coefficients

Let f be a power series ∑a_izⁱ with complex coefficients. The (n. n) Pade approximant to f is a rational function P/Q where P and Q are polynomials, Q(z) ? 0, of degree ≦ n such that f(z)Q(z)-P(z) = Az²ⁿ⁺¹ + higher degree terms. It is proved that if the coefficients a_i satisfy a certain growth condition, then a corresponding subsequence of the sequence of (n, n) Pade approximants converges to f in the region where the power series f converges, except on an exceptional set E having a certain Hausdorff measure 0. It is also proved that the result is best possible in the sense that we may have divergence on E. In particular,there exists an entire function f such that the sequence of (ny n) Pade approximants diverges everywhere (except at 0)     

Characterizations of the polynomial numerical hull of degree k

Six characterizations of the polynomial numerical hull of degree k are established for bounded linear operators on a Hilbert space. It is shown how these characterizations provide a natural distinction between interior and boundary points. One of the characterizations is used to prove that the polynomial numerical hull of any fixed degree k for a Toeplitz matrix whose symbol is piecewise continuous approaches all or most of that of the infinite-dimensional Toeplitz operator, as the matrix size goes to infinity.     

Explicit inverse of a tridiagonal k−Toeplitz matrix

Summary We obtain explicit formulas for the entries of the inverse of a nonsingular and irreducible tridiagonal k–Toeplitz matrix A. The proof is based on results from the theory of orthogonal polynomials and it is shown that the entries of the inverse of such a matrix are given in terms of Chebyshev polynomials of the second kind. We also compute the characteristic polynomial of A which enables us to state some conditions for the existence of A^–1. Our results also extend known results for the case when the residue mod k of the order of A is equal to 0 or k–1 (Numer. Math., 10 (1967), pp. 153–161.).The work was supported by CMUC (Centro de Matemática da Universidade de Coimbra) and by Acção Integrada Luso-Espanhola E-6/03     

Discrete versus continuous Newton's method: A case study

We consider the damped Newton's method N _h (z) = z – hp(z)/p(z), 0<h<1 for polynomialsp(z) with complex coefficients. For the usual Newton's method (h=1) and polynomialsp(z), it is known that the method may fail to converge to a root ofp and rather leads to an attractive periodic cycle.N _h(z) may be interpreted as an Euler step for the differential equation =–p(z)/p(z) with step sizeh. In contrast to the possible failure of Newton's method, we have that for almost all initial conditions to the differential equation that the solutions converge to a root ofp. We show that this property generally carries over to Newton's methodN _h(z) only for certain nondegenerate polynomials and for sufficiently small step sizesh>0. Further we discuss the damped Newton's method applied to the family of polynomials of degree 3.     

Point estimation of simultaneous methods for solving polynomial equations: A survey (II)

The construction of computationally verifiable initial conditions which provide both the guaranteed and fast convergence of the numerical root-finding algorithm is one of the most important problems in solving nonlinear equations. Smale's “point estimation theory” from 1981 was a great advance in this topic; it treats convergence conditions and the domain of convergence in solving an equation f(z)=0$f(z)=0$ using only the information of f   at the initial point _z0$z_0$. The study of a general problem of the construction of initial conditions of practical interest providing guaranteed convergence is very difficult, even in the case of algebraic polynomials. In the light of Smale's point estimation theory, an efficient approach based on some results concerning localization of polynomial zeros and convergent sequences is applied in this paper to iterative methods for the simultaneous determination of simple zeros of polynomials. We state new, improved initial conditions which provide the guaranteed convergence of frequently used simultaneous methods for solving algebraic equations: Ehrlich–Aberth's method, Ehrlich–Aberth's method with Newton's correction, Börsch-Supan's method with Weierstrass’ correction and Halley-like (or Wang–Zheng) method. The introduced concept offers not only a clear insight into the convergence analysis of sequences generated by the considered methods, but also explicitly gives their order of convergence. The stated initial conditions are of significant practical importance since they are computationally verifiable; they depend only on the coefficients of a given polynomial, its degree n and initial approximations to polynomial zeros.     

QR factorization of Toeplitz matrices

Summary This paper presents a new algorithm for computing theQR factorization of anm×n Toeplitz matrix inO(mn) operations. The algorithm exploits the procedure for the rank-1 modification and the fact that both principal (m–1)×(n–1) submatrices of the Toeplitz matrix are identical. An efficient parallel implementation of the algorithm is possible.     

A direct method to solve block banded block Toeplitz systems with non-banded Toeplitz blocks

A fast solution algorithm is proposed for solving block banded block Toeplitz systems with non-banded Toeplitz blocks. The algorithm constructs the circulant transformation of a given Toeplitz system and then by means of the Sherman-Morrison-Woodbury formula transforms its inverse to an inverse of the original matrix. The block circulant matrix with Toeplitz blocks is converted to a block diagonal matrix with Toeplitz blocks, and the resulting Toeplitz systems are solved by means of a fast Toeplitz solver.The computational complexity in the case one uses fast Toeplitz solvers is equal to ξ(m,n,k)=O(mn³)+O(k³n³) flops, there are m block rows and m block columns in the matrix, n is the order of blocks, 2k+1 is the bandwidth. The validity of the approach is illustrated by numerical experiments.     

Complexity reduction of least squares problems involving special vandermonde matrices

This paper presents a new QRD factorization of a rectangular Vandermonde matrix for a special point distribution, including the symmetric case, based on ak-dimensional block decomposition of the matrix and some properties of the Kronecker product. The computational reduction factor with respect to any QR method isk ², in the general case, and 4 in the symmetric one. By the resulting matrix factorization, new formulas are devised for the least squares system solution, whose implementation produces an algorithm of reduced computational cost and computer storage. Finally the perturbation bounds of this new factorization are devised.     

A communication-less parallel algorithm for tridiagonal Toeplitz systems

Diagonally dominant tridiagonal Toeplitz systems of linear equations arise in many application areas and have been well studied in the past. Modern interest in numerical linear algebra is often focusing on solving classic problems in parallel. In McNally [Fast parallel algorithms for tri-diagonal symmetric Toeplitz systems, MCS Thesis, University of New Brunswick, Saint John, 1999], an m processor Split & Correct algorithm was presented for approximating the solution to a symmetric tridiagonal Toeplitz linear system of equations. Nemani [Perturbation methods for circulant-banded systems and their parallel implementation, Ph.D. Thesis, University of New Brunswick, Saint John, 2001] and McNally (2003) adapted the works of Rojo [A new method for solving symmetric circulant tri-diagonal system of linear equations, Comput. Math. Appl. 20 (1990) 61–67], Yan and Chung [A fast algorithm for solving special tri-diagonal systems, Computing 52 (1994) 203–211] and McNally et al. [A split-correct parallel algorithm for solving tri-diagonal symmetric Toeplitz systems, Internat. J. Comput. Math. 75 (2000) 303–313] to the non-symmetric case. In this paper we present relevant background from these methods and then introduce an m processor scalable communication-less approximation algorithm for solving a diagonally dominant tridiagonal Toeplitz system of linear equations.     

Block band Toeplitz preconditioners derived from generating function approximations: analysis and applications

We are concerned with the study and the design of optimal preconditioners for ill-conditioned Toeplitz systems that arise from a priori known real-valued nonnegative generating functions f(x,y) having roots of even multiplicities. Our preconditioned matrix is constructed by using a trigonometric polynomial θ(x,y) obtained from Fourier/kernel approximations or from the use of a proper interpolation scheme. Both of the above techniques produce a trigonometric polynomial θ(x,y) which approximates the generating function f(x,y), and hence the preconditioned matrix is forced to have clustered spectrum. As θ(x,y) is chosen to be a trigonometric polynomial, the preconditioner is a block band Toeplitz matrix with Toeplitz blocks, and therefore its inversion does not increase the total complexity of the PCG method. Preconditioning by block Toeplitz matrices has been treated in the literature in several papers. We compare our method with their results and we show the efficiency of our proposal through various numerical experiments.This research was co-funded by the European Union in the framework of the program “Pythagoras I” of the “Operational Program for Education and Initial Vocational Training” of the 3rd Community Support Framework of the Hellenic Ministry of Education, funded by national sources (25%) and by the European Social Fund - ESF (75%). The work of the second and of the third author was partially supported by MIUR (Italian Ministry of University and Research), grant number 2004015437.     

On the new fourth-order methods for the simultaneous approximation of polynomial zeros

A new iterative method of the fourth-order for the simultaneous determination of polynomial zeros is proposed. This method is based on a suitable zero-relation derived from the fourth-order method for a single zero belonging to the Schröder basic sequence. One of the most important problems in solving polynomial equations, the construction of initial conditions that enable both guaranteed and fast convergence, is studied in detail for the proposed method. These conditions are computationally verifiable since they depend only on initial approximations, the polynomial coefficients and the polynomial degree, which is of practical importance. The construction of improved methods in ordinary complex arithmetic and complex circular arithmetic is discussed. Finally, numerical examples and the comparison with existing fourth-order methods are given.     

Pivoting strategy for rank-one modification ofLDM t-like factorization

The rank-one modification algorithm of theLDM ^t factorization was given by Bennett [1]. His method, however, could break down even when the matrix is nonsingular and well-conditioned. We introduce a pivoting strategy for avoiding possible break-down as well as for suppressing error growth in the modification process. The method is based on a symbolic formula of the rank-one modification of the factorization of a possibly singular nonsymmetric matrix. A new symbolic formula is also obtained for the inverses of the factor matrices. Repeated application of our method produces theLDM ^t-like product form factorization of a matrix. A numerical example is given to illustrate our pivoting method. An incomplete factorization algorithm is also introduced for updating positive definite matrix useful in quasi-Newton methods, in which the Fletcher and Powell algorithm [2] and the Gill, Murray and Saunders algorithm [4] are usually used.This paper is presented at the Japan SIAM Annual Meeting held at University of Tokyo, Japan, October 7–9, 1991.     

Matrix computations and polynomial root-finding with preprocessing

Solution of homogeneous linear systems of equations is a basic operation of matrix computations. The customary algorithms rely on pivoting, orthogonalization and SVD, but we employ randomized preprocessing instead. This enables us to accelerate the solution dramatically, both in terms of the estimated arithmetic cost and the observed CPU time. The approach is effective in the cases of both general and structured input matrices and we extend it and its computational advantages to the solution of nonhomogeneous linear systems of equations, matrix eigen-solving, the solution of polynomial and secular equations, and approximation of a matrix by a nearby matrix that has a smaller rank or a fixed structure (e.g., of the Toeplitz or Hankel type). Our analysis and extensive experiments show the power of the presented algorithms.     

A capacitance matrix method for Dirichlet problem on polygon region

Summary An efficient algorithm for the solution of linear equations arising in a finite element method for the Dirichlet problem is given. The cost of the algorithm is proportional toN ²log₂ N (N=1/h) where the cost of solving the capacitance matrix equations isNlog₂ N on regular grids andN ^3/2log₂ N on irregular ones.     

Fast solution of unsymmetric banded Toeplitz systems by means of spectral factorizations and Woodbury's formula

A fast algorithm for solving systems of linear equations with banded Toeplitz matrices is studied. An important step in the algorithm is a novel method for the spectral factorization of the generating function associated with the Toeplitz matrix. The spectral factorization is extracted from the right deflating subspaces corresponding to the eigenvalues inside and outside the open unit disk of a companion matrix pencil constructed from the coefficients of the generating function. The factorization is followed by the Woodbury inversion formula and solution of several banded triangular systems. Stability of the algorithm is discussed and its performance is demonstrated by numerical experiments. Copyright © 2012 John Wiley & Sons, Ltd.     

Vandermonde matrices on integer nodes: the rectangular case

Summary. We consider a rectangular Vandermonde matrix V on integer nodes. Using combinatorial identities, we give an explicit Cholesky factorization of the matrix , a factorization of the pseudo-inverse of V, and asymptotic estimates of the extremal eigenvalues of B. We also discuss the numerical properties of the proposed formulas. Received March 9, 1999 / Revised version received July 12, 1999 / Published online August 24, 2000     

Block LU factorizations of M–matrices

Summary. It is well known that any nonsingular M–matrix admits an LU factorization into M–matrices (with L and U lower and upper triangular respectively) and any singular M–matrix is permutation similar to an M–matrix which admits an LU factorization into M–matrices. Varga and Cai establish necessary and sufficient conditions for a singular M–matrix (without permutation) to allow an LU factorization with L nonsingular. We generalize these results in two directions. First, we find necessary and sufficient conditions for the existence of an LU factorization of a singular M-matrix where L and U are both permitted to be singular. Second, we establish the minimal block structure that a block LU factorization of a singular M–matrix can have when L and U are M–matrices. Received November 21, 1994 / Revised version received August 4, 1997     

On neighbouring matrices with quadratic elementary divisors

Summary Algorithms are presented which compute theQR factorization of an order-n Toeplitz matrix inO(n ²) operations. The first algorithm computes onlyR explicitly, and the second computes bothQ andR. The algorithms are derived from a well-known procedure for performing the rank-1 update ofQR factors, using the shift-invariance property of the Toeplitz matrix. The algorithms can be used to solve the Toeplitz least-squares problem, and can be modified to solve Toeplitz systems inO(n) space.