This paper describes a new O(n log 3 (n)) solver for the positive definite Toeplitz system T x = b. Instead of computing generators for the inverse of T , the new algorithm adjoins b to T and applies a superfast Schur algorithm to the resulting augmented matrix. The generators of this augmented matrix and its Schur complements are used by a divide-and-conquer block back-substitution routine to complete the solution of the system. The goal is to avoid the well-known numerical instability

more »

... in explicit inversion. Experiments suggest that the algorithm is backward stable in most cases. and successively solve the two smaller systems (2) When dividing the system T x = b into two systems, it will simplify the presentation to assume that we have divided them in half. Thus T 11 is m/2 × m/2. We will assume that T so partitioned. A superfast Schur algorithm applied to the augmented matrix formed from T and b can be used to compute the triangular system (1). The block back-substitution (2) is nothing more than the solution of two smaller Toeplitz-like systems that can be combined to solve the full system using a divide-and-conquer procedure. The multiplication T 12 x 2 can be performed using the FFT. The resulting algorithm avoids the suspect step of multiplication by T −1 but at the cost of increasing the complexity from O(n log 2 (n)) to O(n log 3 (n)) floating point operations. The increase in complexity occurs because the second system of (2) involves a modified right-hand-side b 1 − T 12 x 2 and consequently a modified augmented matrix. The new augmented matrix requires its own O(n log 2 (n)) superfast Schur factorization so that the overall procedure is O(n log 3 (n)). The algorithm of [2] is the model for the derivation of the new algorithm as well as the benchmark for evaluating stability and efficiency. In the remainder of this section we will describe both the algorithm of [2] and the generalized Schur algorithm for a matrix with arbitrary displacement rank. In §2 we show how the same divide-and-conquer idea can be applied to an augmented system that incorporates the right hand side vector b. In §3 we show how the information computed by a superfast block triangularization of the augmented matrix can be used to solve the system T x = b without the need for explicit matrix inversion. In §4 we evaluate the computational complexity of the algorithm. In §5 we present the result of numerical experiments that demonstrate the improved stability of the algorithm. Finally, in §6 we make some observations on the possibility of a proof of numerical stability and compare the new method to another stabilized superfast algorithm.