Newton's iteration for structured matrices and linear systems of equations

[abstract]

Preconditioned conjugate gradient (PCG) methods are quite effective for the solution of a large class of Toeplitz and Toeplitz-like linear systems of equations. In a few iterations they provide an approximate solution and then improve the approximations with linear convergence rate. While this performance is sufficient in several applications, there are many cases where it is still desirable to compute a highly accurate solution at a faster convergence rate, assuming that a rough initial approximation has been computed by a PCG algorithm. Alternatively, any Toeplitz and Toeplitz-like linear system of equations can be effectively solved by the known direct methods, but in many cases the computed solution must be refined to counter the effect of rounding errors.

We describe and analyze four quadratically convergent algorithms for the refinement of rough initial approximations to the inverses of n x n nonsingular Toeplitz and Toeplitz-like matrices and to the solutions of nonsingular Toeplitz and Toeplitz-like linear systems of n equations. The algorithms are variations of Newton's iteration, which we adjust to structured matrix computations based on the known inversion formulas for Toeplitz matrices and the displacement representations of Toeplitz-like matrices. Due to using matrix structure, we perform each iteration step in O(n log n) time. Numerical experiments confirm the results of our analysis and the efficiency of the algorithms. The algorithms can be extended to the case of structured matrices of other classes.

We rely on the variations of our earlier techniques to compute the desired refinement by means of Newton's iteration, which represents the class of the residual correction methods and which we simplify in the case of Toeplitz and Toeplitz-like matrices. Our first two modifications of Newton's iteration (Algs. 1 and 2) exploit the displacement structure of the input matrices to simplify the computations. They work in the more general case of a Toeplitz-like input. Our third and fourth modifications (Algs.3 and 4) specialize Newton's iteration as the residual correction of the solution of a Toeplitz linear system of equations. They rely on the inversion formulae known for Toeplitz matrices and simplify Algs. 1 and 2, respectively.

Alg. 1 is a little more costly to perform but implies the convergence under milder assumptions about the initial approximation than Alg. 2. Alg. 4 runs roughly twice as fast as Alg. 3 and does not seem to require any stronger initial assumptions. Algs. 1 and 3 are more convenient to apply where the triangular Toeplitz representation is used for the input and output matrices, whereas Algs. 2 and 4 are more convenient to apply where a factor-circulant (/-circulant) representation is used. Some of our algorithms can be extended to other classes of structured matrices.

All presented algorithms require an initial approximation lying sufficiently close to the solution, for otherwise the iterations may converge too slowly or diverge. A partial remedy can be obtained by means of the homotopy techniques.