Stability analysis and fast algorithms for triangularization of Toeplitz matrices

Summary. A new algorithm for triangularizing an Toeplitz matrix is presented. The algorithm is based on the previously developed recursive algorithms that exploit the Toeplitz structure and compute each row of the triangular factor via updating and downdating steps. A forward error analysis for this existing recursive algorithm is presented, which allows us to monitor the conditioning of the problem, and use the method of corrected semi-normal equations to obtain higher accuracy for certain ill-conditioned matrices. Numerical experiments show that the new algorithm improves the accuracy significantly while the computational complexity stays in . Received April 30, 1995 / Revised version received February 12, 1996     

Symmetric Gaussian elimination for cauchy-type matrices with application to positive definite Toeplitz matrices

The problem of solving linear equations with a Toeplitz matrix appears in many applications. Often is positive definite but ill-conditioned with many small eigenvalues. In this case fast and superfast algorithms may show a very poor behavior or even break down. In recent papers the transformation of a Toeplitz matrix into a Cauchy-type matrix is proposed. The resulting new linear equations can be solved in operations using standard pivoting strategies which leads to very stable fast methods also for ill-conditioned systems. The basic tool is the formulation of Gaussian elimination for matrices with low displacement rank. In this paper, we will transform a Hermitian Toeplitz matrix into a Cauchy-type matrix by applying the Fourier transform. We will prove some useful properties of and formulate a symmetric Gaussian elimination algorithm for positive definite . Using the symmetry and persymmetry of we can reduce the total costs of this algorithm compared with unsymmetric Gaussian elimination. For complex Hermitian , the complexity of the new algorithm is then nearly the same as for the Schur algorithm. Furthermore, it is possible to include some strategies for ill-conditioned positive definite matrices that are well-known in optimization. Numerical examples show that this new algorithm is fast and reliable. Received March 24, 1995 / Revised version received December 13, 1995     

An acceleration theorem for the ρ-algorithm

Summary. The ρ-algorithm of Wynn is an excellent device for accelerating the convergence of some logarithmically convergent sequences. Until now a convergence theorem and an acceleration theorem for the ρ-algorithm have not been obtained. The purpose of this paper is to give an acceleration theorem for the ρ-algorithm. Moreover, it is proved that the ρ-algorithm cannot accelerate linear convergence. Numerical examples are given. Received October 20, 1994 / Revised version received July 2, 1995     

Inverse unitary eigenproblems and related orthogonal functions

Summary. This paper explores the relationship between certain inverse unitary eigenvalue problems and orthogonal functions. In particular, the inverse eigenvalue problems for unitary Hessenberg matrices and for Schur parameter pencils are considered. The Szeg? recursion is known to be identical to the Arnoldi process and can be seen as an algorithm for solving an inverse unitary Hessenberg eigenvalue problem. Reformulation of this inverse unitary Hessenberg eigenvalue problem yields an inverse eigenvalue problem for Schur parameter pencils. It is shown that solving this inverse eigenvalue problem is equivalent to computing Laurent polynomials orthogonal on the unit circle. Efficient and reliable algorithms for solving the inverse unitary eigenvalue problems are given which require only O() arithmetic operations as compared with O() operations needed for algorithms that ignore the structure of the problem. Received April 3, 1995 / Revised version received August 29, 1996     

A central limit theorem for stationary random fields

Summary. We prove a central limit theorem for strictly stationary random fields under a projective assumption. Our criterion is similar to projective criteria for stationary sequences derived from Gordin's theorem about approximating martingales. However our approach is completely different, for we establish our result by adapting Lindeberg's method. The criterion that it provides is weaker than martingale-type conditions, and moreover we obtain as a straightforward consequence, central limit theorems for α-mixing or φ-mixing random fields. Received: 19 February 1997 / In revised form: 2 September 1997     

Weierstrass formula and zero-finding methods

Summary. Classical Weierstrass' formula [29] has been often the subject of investigation of many authors. In this paper we give some further applications of this formula for finding the zeros of polynomials and analytic functions. We are concerned with the problems of localization of polynomial zeros and the construction of iterative methods for the simultaneous approximation and inclusion of these zeros. Conditions for the safe convergence of Weierstrass' method, depending only on initial approximations, are given. In particular, we study polynomials with interval coefficients. Using an interval version of Weierstrass' method enclosures in the form of disks for the complex-valued set containing all zeros of a polynomial with varying coefficients are obtained. We also present Weierstrass-like algorithm for approximating, simultaneously, all zeros of a class of analytic functions in a given closed region. To demonstrate the proposed algorithms, three numerical examples are included. Received September 13, 1993     

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     

Pivoting strategies leading to diagonal dominance by rows

Pivoting strategies for Gaussian elimination leading to upper triangular matrices which are diagonally dominant by rows are studied. Forward error analysis of triangular systems whose coefficient matrices are diagonally dominant by rows is performed. We also obtain small bounds of the backward errors for the pivoting strategies mentioned above. Our examples of matrices include H-matrices and some generalizations of diagonally dominant matrices, and scaled partial pivoting for the 1-norm is an example of these pivoting strategies. In the case of an M-matrix, a pivoting strategy of computational complexity is proposed, which satisfies all the results of the paper. Received June 6, 1997 / Revised version received October 27, 1997     

Numerical solution of a secular equation

Summary. A method is proposed for the solution of a secular equation, arising in modified symmetric eigenvalue problems and in several other areas. This equation has singularities which make the application of standard root-finding methods difficult. In order to solve the equation, a class of transformations of variables is considered, which transform the equation into one for which Newton's method converges from any point in a certain given interval. In addition, the form of the transformed equation suggests a convergence accelerating modification of Newton's method. The same ideas are applied to the secant method and numerical results are presented. Received July 1, 1994     

Relative-error bounds for the LU decomposition via the GTH algorithm

Summary. Recently the author showed that the Grassmann-Taksar-Heyman (GTH) algorithm computes the steady-state distribution of a finite-state Markov chain with low relative error. Here it is shown that the LU decomposition computed in the course of the GTH algorithm also has low relative error. The proof requires a refinement of the methods used in the earlier paper. Received September 2, 1994 / Revised version received July 17, 1995     

Analysis of the quasi-Laguerre method

The quasi-Laguerre's iteration formula, using first order logarithmic derivatives at two points, is derived for finding roots of polynomials. Three different derivations are presented, each revealing some different properties of the method. For polynomials with only real roots, the method is shown to be optimal, and the global and monotone convergence, as well as the non-overshooting property, of the method is justified. Different ways of forming quasi-Laguerre's iteration sequence are addressed. Local convergence of the method is proved for general polynomials that may have complex roots and the order of convergence is . Received June 30, 1996 / Revised version received August 12, 1996     

A look-ahead Bareiss algorithm for general Toeplitz matrices

Summary. The Bareiss algorithm is one of the classical fast solvers for systems of linear equations with Toeplitz coefficient matrices. The method takes advantage of the special structure, and it computes the solution of a Toeplitz system of order~ with only~ arithmetic operations, instead of~ operations. However, the original Bareiss algorithm requires that all leading principal submatrices be nonsingular, and the algorithm is numerically unstable if singular or ill-conditioned submatrices occur. In this paper, an extension of the Bareiss algorithm to general Toeplitz systems is presented. Using look-ahead techniques, the proposed algorithm can skip over arbitrary blocks of singular or ill-conditioned submatrices, and at the same time, it still fully exploits the Toeplitz structure. Implementation details and operations counts are given, and numerical experiments are reported. We also discuss special versions of the proposed look-ahead Bareiss algorithm for Hermitian indefinite Toeplitz systems and banded Toeplitz systems. Received August 27, 1993 / Revised version received March 1994     

On the convergence of multipoint iterations

Summary. This note gives a new convergence proof for iterations based on multipoint formulas. It rests on the very general assumption that if the desired fixed point appears as an argument in the formula then the formula returns the fixed point. Received March 24, 1993 / Revised version received January 1994     

Tangent Graeffe iteration

Summary.   Graeffe iteration was the choice algorithm for solving univariate polynomials in the XIX-th and early XX-th century. In this paper, a new variation of Graeffe iteration is given, suitable to IEEE floating-point arithmetics of modern digital computers. We prove that under a certain generic assumption the proposed algorithm converges. We also estimate the error after N iterations and the running cost. The main ideas from which this algorithm is built are: classical Graeffe iteration and Newton Diagrams, changes of scale (renormalization), and replacement of a difference technique by a differentiation one. The algorithm was implemented successfully and a number of numerical experiments are displayed. Received May 29, 1998 / Revised version received September 13, 1999 / Published online April 5, 2001     

A new line-search step based on the Weierstrass \wp-function for minimizing a class of logarithmic barrier functions

Summary. We define the notion of self-concordance of order two for the restriction of a logarithmic barrier function to a given line. Based on this notion we prove an inner approximation of the domain of , as well as a lower bound of the distance from a point $t$ to the minimum of . These results provide the theoretical tools to develop a simple and efficient search step for finding the minimum of the barrier function along a given line. The new bound on the size of the line-search step is better than the optimal bound known for the case of a self-concordant function (of order one). We conclude with some numerical examples that illustrate the promise of the new line-search step. Received May 24, 1993 / Revised version received February 1994     

Fast verification of solutions of matrix equations

Summary. In this paper, we are concerned with a matrix equation where A is an real matrix and x and b are n-vectors. Assume that an approximate solution is given together with an approximate LU decomposition. We will present fast algorithms for proving nonsingularity of A and for calculating rigorous error bounds for . The emphasis is on rigour of the bounds. The purpose of this paper is to propose different algorithms, the fastest with flops computational cost for the verification step, the same as for the LU decomposition. The presented algorithms exclusively use library routines for LU decomposition and for all other matrix and vector operations. Received June 16, 1999 / Revised version received January 25, 2001 / Published online June 20, 2001     

Approximation by translates of refinable functions

Summary. The functions are refinable if they are combinations of the rescaled and translated functions . This is very common in scientific computing on a regular mesh. The space of approximating functions with meshwidth is a subspace of with meshwidth . These refinable spaces have refinable basis functions. The accuracy of the computations depends on , the order of approximation, which is determined by the degree of polynomials that lie in . Most refinable functions (such as scaling functions in the theory of wavelets) have no simple formulas. The functions are known only through the coefficients in the refinement equation – scalars in the traditional case, matrices for multiwavelets. The scalar "sum rules" that determine are well known. We find the conditions on the matrices that yield approximation of order from . These are equivalent to the Strang–Fix conditions on the Fourier transforms , but for refinable functions they can be explicitly verified from the . Received August 31, 1994 / Revised version received May 2, 1995     

Hybrid procedures for solving linear systems

Summary. In this paper, we introduce the notion of hybrid procedures for solving a system of linear equations. A hybrid procedure consists in a combination of two arbitrary approximate solutions with coefficients summing up to one. Thus the combination only depends on one parameter whose value is chosen in order to minimize the Euclidean norm of the residual vector obtained by the hybrid procedure. Properties of such procedures are studied in detail. The two approximate solutions which are combined in a hybrid procedure are usually obtained by two iterative methods. Several strategies for combining these two methods together or with the previous iterate of the hybrid procedure itself are discussed and their properties are analyzed. Numerical experiments illustrate the various procedures. Received October 21, 1992/Revised version received May 28, 1993