Error Analysis of a Partial Pivoting Method for Structured Matrices^{∗}

Douglas R. Sweet Maritime Operations Division

Defence Science and Technology Organisation

Salisbury, SA 5108

           Doug.Sweet@dsto.defence.gov.au

Richard P. Brent Computer Sciences Laboratory Australian National University Canberra, ACT 0200

                                                          rpb@cslab.anu.edu.au

Report TR-CS-95-03 7 June 1995

Abstract

Many matrices that arise in the solution of signal processing problems have a special displacement structure. For example, adaptive filtering and direction-of-arrival estimation yield matrices of Toeplitz type. A recent method of Gohberg, Kailath and Olshevsky (GKO) allows fast Gaussian elimination with partial pivoting for such structured matrices. In this paper, a rounding error analysis is performed on the Cauchy and Toeplitz variants of the GKO method. It is shown the error growth depends on the growth in certain auxiliary vectors, the generators, which are computed by the GKO algorithms. It is also shown that in certain circumstances, the growth in the generators can be large, and so the error growth is much larger than would be encountered with normal Gaussian elimination with partial pivoting. A modification of the algorithm to perform a type of row-column pivoting is proposed which may ameliorate this problem.

Keywords: Structured matrices, fast algorithms, displacement rank, generators, pivoting, error analysis, stability

1 Introduction

Many problems which occur in signal processing, control theory and interpolation lead to a square or rectangular system with a special structure, for which the exact or least-squares solution is required. For example, adaptive filtering requires either the exact solution of a square Toeplitz system or the least-squares solution of a rectangular Toeplitz system. A Toeplitz matrix is one whose entries along the NW to SE diagonals are constant, i.e. element t_{ij} depends only on i − j. Other types of structured matrices which arise are Hankel matrices whose entries along the SW to NE diagonals are constant, Vandermonde matrices whose entries have the form v_{ij} = x_{i}^{j−1}, and Cauchy matrices whose entries have the form c_{ij} = 1/(t_{i} − s_{j}), where the t_{i} and s_{j} are the elements of vectors t and s.

Normally, the exact or least-squares solution of a linear system requires O(n^{3}) operations to solve, where n is the order of the system. However, the structure of the systems mentioned above has been exploited in the past [1, 4, 7] to derive fast solvers, i.e. those that require O(n^{2}) or fewer operations. These fast algorithms are in general numerically unstable for indefinite systems [3, 5, 11]. Recently, methods have been proposed [6, 10, 11] which are numerically stable, but which attempt to retain the O(n^{2}) complexity. However, all of these algorithms will require O(n^{3}) operations in the worst case. The BBH method [2] requires O(n^{2}) operations in the worst case and can be shown to be weakly stable, but not stable in the usual sense of backward error analysis. Thus there is an interest in fast algorithms which require O(n^{2}) operations in the worst case and can be shown to be stable.

∗Copyright \begin{matrix}c \\ ◯ \end{matrix} 1995, the authors. Shorter version to appear in Advanced Signal Processing Algorithms, Proc.

SPIE 40th Annual Meeting, San Diego, July 1995. rpb157tr

Recently, Gohberg, Kailath and Olshevsky [8] have shown how to perform Gaussian elimination in a fast way with matrices with a special displacement structure. Such matrices include Toeplitz, Vandermonde, Hankel and Cauchy matrices, and generalizations thereof, called Toeplitz type, etc. They also show how to incorporate partial pivoting into the Cauchy and Vandermonde solvers. They point out that although pivoting cannot be incorporated directly into the corresponding Toeplitz or Hankel solvers, the Toeplitz and Hankel problems can be transformed by simple orthogonal operations into Cauchy problems. The solution to the original systems can be recovered from those of the transformed systems by the reverse orthogonal operations. Thus fast Gaussian elimination with partial pivoting can be carried out on Toeplitz, Vandermonde, Hankel and Cauchy systems.

It might be assumed that such fast solvers should have the same stability properties as Gaussian elimination with partial pivoting. One of the aims of this paper is analyse the error behaviour of these algorithms by means of a backward error analysis. It is shown that error propagation depends on the magnitude of both the triangular factors L and U (as in Gaussian elimination) and the generators, auxiliary vectors which are computed during the course of the algorithm.

It is shown that in some cases the generators can suffer a large growth and cause a corresponding growth in the backward and forward error. A modification is proposed which may prevent this growth, and so restore the stability of the algorithm in these cases. However, we can not prove that the modification is always successful.

The paper is structured as follows. In §2, the Gohberg-Kailath-Olshevsky (GKO) algorithm for Cauchy and Toeplitz matrices is briefly described. The error analyses of the Cauchy and Toeplitz variants of the GKO algorithm are carried out in §3 and §4 respectively. In §5, examples for both variants are given where a large growth occurs in the generators and hence in the errors in the solutions. The modified version of the GKO algorithm is proposed in §6, and numerical tests of this are carried out there. Some conclusions are drawn and suggestions for future work are given in §7.

Notation. The following notation is used. ϵ is the machine epsilon, and n is the order of the matrix to be factorized. Scalars of the form c_{i} and k_{i} are small constants. e_{j} denotes the jth column of the identity matrix. Elementwise matrix multiplication is denoted by the centred circle ∘. For a matrix A, |A| is the matrix of moduli of the {a_{ij}}, A^{I} denotes elementwise inversion, and A′ denotes augmentation of A to order n by adding zero rows and zero columns respectively above and to the left of A. Other submatrices are indicated in MATLAB style, i.e. for a matrix A, A_{p:q,r:s} selects rows p to q of columns r to s, and a colon without an index range selects all of the rows or columns.

2 The Gohberg-Kailath-Olshevsky (GKO) Algorithm

In this section, we first define the displacement operator, displacement equation and displacement rank for structured matrices; we then give the general Gaussian elimination algorithm for structured matrices, followed by the variants for Cauchy and Toeplitz matrices.

2.1 Displacement structure

Gohberg et al [8] show that structured matrices satisfy a Sylvester equation which has the form

∇_{{A_{f} ,A_{b}}}(R) = A_{f}R − RA_{b} = ΦΨ ,  (1)

where A_{f} and A_{b} have some simple structure (usually banded, with 3 or fewer full diagonals), Φ and Ψ are n × α and α × n respectively, and α is some small integer (usually 4 or less). The pair of matrices Φ, Ψ is called the {A_{f}, A_{b}}-generator of R, and α is called the {A_{f}, A_{b}}-displacement rank of R.

Particular choices of A_{f} and A_{b} lead to definitions of basic classes of matrices. Thus, for a Cauchy matrix

C(t, s) = [\frac{1}{t_{i} − s_{j}} ]_{ij} ,

we have

A_{f} = D_{t} = diag(t_{1}, t_{2}, …, t_{n}), A_{b} = D_{s} = diag(s_{1}, s_{2}, …, s_{n})  (2)

and

Φ^{T} = Ψ = [1, 1, …, 1] .

More general matrices, where A_{f} and A_{b} are as in (2) but Φ and Ψ are general rank-α matrices, are called Cauchy-type.

Similarly, for a Toeplitz matrix T = [t_{ij}] = [a_{i−j}]

A_{f} = Z_{1} = [ \begin{matrix}0 & 0 & ⋯ & 0 & 1 \\ 1 & 0 & & & 0 \\ 0 & 1 & & & ⋮ \\ ⋮ & & ... & & ⋮ \\ 0 & ⋯ & 0 & 1 & 0 \end{matrix} ] , A_{b} = Z_{−1} = [ \begin{matrix}0 & 0 & ⋯ & 0 & −1 \\ 1 & 0 & & & 0 \\ 0 & 1 & & & ⋮ \\ ⋮ & & ... & & ⋮ \\ 0 & ⋯ & 0 & 1 & 0 \end{matrix} ] ,  (3)

Φ = [ \begin{matrix}1 & 0 & ⋯ & ⋯ & 0 \\ a_{0} & a_{1−n} + a_{1} & ⋯ & a_{−2} + a_{n−2} & a_{−1} + a_{n−1} \end{matrix} ]^{T}  (4)

and

Ψ = [ \begin{matrix}a_{n−1} − a_{−1} & a_{n−2} − a_{−2} & ⋯ & a_{1} − a_{1−n} & a_{0} \\ 0 & ⋯ & ⋯ & 0 & 1 \end{matrix} ] .  (5)

2.2 Gaussian elimination for structured matrices

Let the input matrix, R_{1}, have the partitioning R_{1} = [ \begin{matrix}d_{1} & w_{1}^{T} \\ y_{1} & Ṙ_{1} \end{matrix} ]. The first step of normal

Gaussian elimination is to premultiply R_{1} by [ \begin{matrix}1 & 0^{T} \\ −y_{1}/d_{1} & I \end{matrix} ], which reduces R_{1} to [ \begin{matrix}d_{1} & w_{1}^{T} \\ 0 & R_{2} \end{matrix} ], where R_{2} = Ṙ_{1} − y_{1}w_{1}^{T} /d_{1} is the Schur complement of d_{1} in R_{1}. At this stage, R_{1} has the factorization

R_{1} = [ \begin{matrix}1 & 0^{T} \\ y_{1}/d_{1} & I \end{matrix} ] [ \begin{matrix}d_{1} & w_{1}^{T} \\ 0 & R_{2} \end{matrix} ] .

One then proceeds recursively with the Schur complement R_{2} = [ \begin{matrix}d_{2} & w_{2}^{T} \\ y_{2} & Ṙ_{2} \end{matrix} ], eventually yielding

a factorization R_{1} = LU, where column k of L is [0^{T} 1 y_{k}^{T}]^{T}, and row k of U is [0^{T} 1 w_{k}^{T}].

The genesis of structured Gaussian elimination is the fact that the displacement structure is preserved under Schur complementation, and that the generators for the Schur complement R_{k+1} can be computed from the generators of R_{k} in O(n) operations. This is expressed constructively in the following theorem, which is proved in [8].

Theorem 2.1 Let matrix R_{1} = [ \begin{matrix}d_{1} & w_{1}^{T} \\ y_{1} & Ṙ_{1} \end{matrix} ] satisfy the Sylvester equation

∇_{{A_{f,1},A_{b,1}}}(R_{1}) = A_{f,1}R_{1} − R_{1}A_{b,1} = Φ^{(1)}Ψ^{(1)} ,  (6)

where Φ^{(1)} = [φ_{1}^{(1)T} φ_{2}^{(1)T} ⋯ φ_{n}^{(1)T}]^{T}, Ψ^{(1)} = [ψ_{1}^{(1)} ψ_{2}^{(1)} ⋯ ψ_{n}^{(1)}], φ_{i}^{(1)} ∈ C^{1×α} and ψ_{i}^{(1)} ∈ C^{1×α}, (i = 1, 2, …, n). Then R_{2}, the Schur complement of d_{1} in R_{1}, satisfies the Sylvester equation

∇_{{A_{f,2},A_{b,2}}}(R_{2}) = A_{f,2}R_{2} − R_{2}A_{b,2} = Φ^{(2)}Ψ^{(2)} ,

where A_{f,2} and A_{b,2} are respectively A_{f,1} and A_{b,1} with their first rows and first columns deleted, and where Φ^{(2)} = [0, φ_{2}^{(2)T} , φ_{3}^{(2)T} , ⋯, φ_{n}^{(2)T}]^{T} and Ψ^{(2)} = [0, ψ_{2}^{(2)}, ψ_{3}^{(2)}, ⋯, ψ_{n}^{(2)}] are given by

Φ_{2:n,:}^{(2)} = Φ_{2:n,:}^{(1)} − y_{1}φ_{1}^{(1)}/d_{1} ,  (7)
Ψ_{:,2:n}^{(2)} = Ψ_{:,2:n}^{(1)} − ψ_{1}^{(1)}w_{1}^{T} /d_{1} .  (8)

Equations (7) and (8) form the basis of the following general structured Gaussian elimination algorithm.

Algorithm 2.1 (Structured Gaussian elimination)

1. Recover from the generator Φ^{(1)}, Ψ^{(1)} the first row and column of R_{1} = [ \begin{matrix}d_{1} & w_{1}^{T} \\ y_{1} & R_{22}^{(1)} \end{matrix} ] .

2. [1 y_{1}^{T} /d_{1}]^{T} and [d_{1} w_{1}^{T}] are respectively the first column and row of L_{1} and U_{1} in the LU factorization of R_{1}.

3. Compute by equations (7) and (8), the generator Φ^{(2)}, Ψ^{(2)} for the Schur complement R_{2}.

4. Proceed recursively with Φ^{(2)} and Ψ^{(2)}. Each major step yields [1 y_{k}^{T} /d_{k}]^{T} and [d_{k} w_{k}^{T}], which are respectively the first column and row of L_{k} and U_{k} in the LU factorization of R_{k}. Column k of L and row k of U are respectively [0_{k−1}^{T} 1 y_{k}^{T} /d_{k}]^{T} and [0_{k−1}^{T} d_{k} w_{k}^{T}].

Pivoting. Gaussian elimination without pivoting is unstable in general. One normally uses partial pivoting (swapping rows to bring the largest element in the first column to the pivot position) or complete pivoting (swapping rows and columns to bring the largest element in the whole matrix to the pivot position) to improve the accuracy. Row and/or column interchanges can destroy the structure of certain matrices, such as Toeplitz matrices. However, if A_{f,1} in (6) is diagonal (which is the case for Cauchy and Vandermonde type matrices), then the structure is preserved under row permutations.

Partial pivoting can also be incorporated into structured Gaussian elimination. Suppose we wish to swap rows 1 and q of R_{1}. Let P_{1} be the matrix which applies this permutation. Then it is easy to see that P_{1}R_{1} satisfies (6) with the (1, 1) and (q, q) entries of A_{f,1} swapped, and with swapped row vectors φ_{1}^{(1)} and φ_{q}^{(1)}. Thus, pivoting can be incorporated into Algorithm 2.1 by adding the following steps:

0.5 Initialization step. Set permutation matrix P = I.

2.5 After step 2 of Algorithm 2.1. Let (y_{1})_{q} be the largest entry by magnitude in y_{1}. Swap rows 1 and q of P and Φ^{(1)}, and the first and q-th diagonal entries in A_{f,1}. Recover the first row of P_{1}R_{1} from Ψ^{(1)} and the swapped Φ^{(1)}.

Note. The computation of the first row of the original R^{(1)} in step 1 of Algorithm 2.1 may be omitted – we only require the first row of the swapped R^{(1)}.

It may be seen that the pivoted algorithm computes upper and lower triangular matrices L and U which satisfy

R^{(1)} = P^{T} LU .

Note that for Cauchy-type matrices, where both A_{f,1} and A_{b,1} are diagonal, both row and column pivoting may be performed. However, complete pivoting requires the computation of all the entries in the matrix, which would require O(n^{2}) operations at each step and O(n^{3}) operations in all. It will be seen in §6 that a restricted version of row-column pivoting can be used to improve the performance of the GKO algorithm.

2.3 The Cauchy variant of the GKO algorithm (GKO-Cauchy)

Recall that a Cauchy-type matrix satisfies the Sylvester equation (6) with

A_{f,1} = D_{t} = diag(t_{1}, t_{2}, …, t_{n}) and A_{b,1} = D_{s} = diag(s_{1}, s_{2}, …, s_{n}) .

It can be easily verified that if t_{i} ≠ s_{j}, then the (i, j) entry of R^{(1)} = R is given by

r_{ij} = \frac{φ_{i}ψ_{j}}{t_{i} − s_{j}} .

There may be some cases where t_{i} = s_{j} and φ_{i}ψ_{j} = 0 for some (i, j), and r_{ij} cannot be recovered from its generator. We do not consider these cases in this paper.

In general, at major step k, the reduced matrix R^{(k)} has the form

R^{(k)} = [ \begin{matrix}d_{1} & & & w_{1}^{T} \\ 0 & ... & & ⋮ \\ ⋮ & & d_{k−1} & w_{k−1}^{T} \\ 0 & ⋯ & 0 & R_{k} \end{matrix} ] .

The entries of the k-th Schur component R_{k}, may be computed by

r_{ij}^{(k)} = \frac{φ_{i}^{(k)}ψ_{j}^{(k)}}{t_{i} − s_{j}} , k ≤ i, j ≤ n  (9)
= (R_{k})_{i−k+1,j−k+1}

Equation (9) can be used in Algorithm 2.1 with pivoting to yield the Cauchy version of the GKO algorithm.

Algorithm 2.2 (GKO-Cauchy)

Input. Cauchy-type matrix R_{1}, specified by t, s, Φ^{(1)} and Ψ^{(1)}. Output. Factorization R_{1} = P^{T} LU, where P is a permutation, and L and

U are lower and upper-triangular respectively.

% Initialization L ← 0; U ← 0; P ← I for k ← 1 : n % k: Iteration number

for j ← k : n % recover col.1 of R_{k}

r_{jk}^{(k)} ← \frac{φ_{j}^{(k)}ψ_{k}^{(k)}}{t_{j}−s_{k}} end % Carry out row interchanges Find k ≤ q ≤ n such that |r_{qk}^{(k)}| = max_{k≤j≤n} |r_{jk}^{(k)}| t_{k} ↔ t_{q}; φ_{k}^{(k)} ↔ φ_{q}^{(k)}; r_{kk}^{(k)} ↔ r_{qk}^{(k)} swap k-th and q-th rows of L swap k-th and q-th rows of P

for j ← k + 1 : n % Recover row 1 of swapped R_{k}

r_{kj}^{(k)} ← \frac{φ_{k}^{(k)}ψ_{j}^{(k)}}{t_{k}−s_{j}} end u_{kk} ← r_{kk}^{(k)} % Compute row and col k of L and U, and update Φ and Ψ using (7) and (8) for j ← k + 1 : n

l_{jk} ← r_{jk}^{(k)}/r_{kk}^{(k)}

u_{kj} ← r_{kj}^{(k)}

ψ_{j}^{(k+1)} ← ψ_{j}^{(k)} − ψ_{k}^{(k)}u_{kj}/u_{kk} φ_{j}^{(k+1)} ← φ_{j}^{(k)} − φ_{k}^{(k)}l_{jk}

end end

2.4 The Toeplitz variant of the GKO algorithm (GKO-Toeplitz)

Recall that a Toeplitz matrix satisfies the Sylvester equation (1), with A_{f}, A_{b}, Φ and Ψ being given by equations (3) to (5), and a Toeplitz-type matrix is one with A_{f} and A_{b} given by (3), and with general low-rank Φ and Ψ. The first row and column of the Toeplitz-type matrix can be simply generated from Φ and Ψ, and this generating formula can be used in Algorithm 2.1 to yield a structured Gaussian elimination algorithm for Toeplitz-type matrices.

Because neither A_{f} nor A_{b} is diagonal, pivoting cannot be introduced directly into this structured algorithm – pivoting will destroy the Toeplitz-type property. However, the Toeplitztype matrix can be easily converted, by fast orthogonal transformations, into a Cauchy-type matrix which can be factorized as in Algorithm 2.1. The inverse orthogonal transforms yield the factorization of the original matrix. The following result of [8] shows how this conversion may be done.

Theorem 2.2 Let T be a Toeplitz-type matrix, satisfying

∇_{{Z_{1},Z_{−1}}}(T) = ΩΓ ,
Ω = [ω_{1}^{T} ω_{2}^{T} ⋯ ω_{n}^{T}]^{T} , Γ = [γ_{1} γ_{2} ⋯ γ_{n}],

where the {ω_{i}} and the {γ_{i}} are 1 × α and α × 1 respectively. Then

R = F T D^{−1}F^{∗}  (10)

is a Cauchy-type matrix, satisfying

∇_{{D_{F} ,D_{F} }} = ΦΨ ,

where F = \frac{1}{\sqrt{n}}[e^{2πi(k−1)(j−1)/n}]_{1≤k,j≤n} is the Discrete Fourier Transform matrix,

D_{F} = diag(1, e^{2πi/n}, …, e^{2πi(n−1)/n}) , D_{F} = diag(e^{πi/n}, e^{3πi/n}, …, e^{πi(2n−1)/n}) ,  (11)
D = diag(1, e^{πi/n}, …, e^{πi(n−1)/n})

and

Φ = F Ω ,  Ψ^{∗} = F DΓ^{∗} .  (12)

Theorem 2.2 allows the generators of T to be converted to the generators of R in O(2αn log n) operations via FFTs. R can then be factorized as R = P^{T} LU using Algorithm 2.2. Using (10), we then obtain

T = F^{∗}P^{T} LUF D .  (13)

From this factorization, a linear system in T can be solved in n^{2} + 2n log n operations, so the whole procedure of conversion of Cauchy form, factorization and solution requires O(n^{2}) operations.

3 Error Analysis of the GKO-Cauchy Algorithm

In this section, a backward error analysis will be carried out, which yields a bound for the perturbation matrix E, defined by

L̃Ũ = R + E ,  (14)

where R is the matrix to be factorized, and L̃ and Ũ are the computed factors. In the analysis, we first derive some preliminary results which apply to any algorithm for structured Gaussian elimination (SGE), and indicate a general methodology for error analysis of SGE algorithms. We then carry out the analysis for Cauchy-type matrices in general and for the Cauchy-type matrix derived from a Toeplitz matrix by equation (10).

3.1 Preliminary results

The following two lemmas may be used for the error analysis of SGE algorithms in general, and the GKO-Cauchy algorithm in particular. The first lemma shows that if G is the perturbation in the Sylvester equation caused by replacing R by L̃Ũ, then the displacement of E is G.

Lemma 3.1 Let R be a general structured matrix that satisfies (1), let A_{f}, A_{b}, Φ and Ψ be as defined above, and let L̃, Ũ and E be as in (14). Suppose L̃ and Ũ satisfy

A_{f} L̃Ũ − L̃ŨA_{b} = ΨΦ + G ;  (15)

then E satisfies

∇_{{A_{f} ,A_{b}}}(E) ≡ A_{f}E − EA_{b} = G .  (16)

Proof. From (14) and (15),

A_{f}(R + E) − (R + E)A_{b} = ΦΨ + G .

Expanding the above, and using (1) we obtain (16).

Corollary 3.2 If R is a Cauchy-type matrix with A_{f} = D_{t} and A_{b} = D_{s}, then E satisfies

D_{t}E − ED_{s} = G  (17)

and

e_{ij} = \frac{g_{ij}}{t_{i} − s_{j}} , i, j = 1, …, n  (18)

Proof. (17) follows directly from (16), and (18) follows by evaluating each component of (17).

If R is a Toeplitz-type matrix, A_{f} = Z_{1}, A_{b} = Z_{−1}, and E satisfies Z_{1}E − EZ_{−1} = G. Because Z_{1} and Z_{−1} are not diagonal, the recovery formula for E is a little more involved, and will be derived in the next section (Lemma 4.3).

The second lemma of this section shows that G is the sum of the local perturbation matrices incurred in each step of the relevant structured Gaussian elimination (SGE) algorithm.

Lemma 3.3 Let ∇_{{A_{f} ,A_{b}}} be the displacement operator as defined in (1); let L̃, Ũ and G be as defined above; let the {Φ̃^{(k)}, Ψ̃^{(k)}}_{k=1,2,...} be the computed generators of the {R′_{k}}_{k=1,2,...}, the reduced matrices at step k of SGE, and define Φ̃^{(n+1)} = Ψ̃^{(n+1)} = 0. Then

G = ∑_{k=1}^{n} H_{k} ,  (19)

where H_{k}, the local perturbation in each step of SGE, is defined by

∇_{{A_{f} ,A_{b}}}(l̃_{:k}ũ_{k:}) = Φ̃^{(k)}Ψ̃^{(k)} − Φ̃^{(k+1)}Ψ̃^{(k+1)} + H_{k} , k = 1, …, n .  (20)

Proof. Writing (20) explicitly, we get

A_{f}l̃_{:k}ũ_{k:} − l̃_{:k}ũ_{k:}A_{b} = Φ̃^{(k)}Ψ̃^{(k)} − Φ̃^{(k+1)}Ψ̃^{(k+1)} + H_{k} , k = 1, …, n .  (21)

Summing the members of (21), we obtain

A_{f} ∑_{i=1}^{n} l̃_{:k}ũ_{k:} − ∑_{i=1}^{n} l̃_{:k}ũ_{k:}A_{b} = Φ̃^{(1)}Ψ̃^{(1)} − Φ̃^{(n+1)}Ψ̃^{(n+1)} + ∑_{i=1}^{n} H_{k} .  (22)

Now ∑_{i=1}^{n} l̃_{:k}ũ_{k:} = L̃Ũ, Φ̃^{(1)} = Φ, Ψ̃^{(1)} = Ψ and Φ̃^{(n+1)} ≡ Ψ̃^{(n+1)} ≡ 0. Substituting these identities into (22) and comparing the resulting relation with (15), we obtain (19)

3.2 Methodology of error analysis for SGE algorithms

Lemmas 3.1 and 3.3 may be used in a general methodology for the error analysis of SGE algorithms similar to Algorithm 2.1.

In the following methodology and the subsequent analysis of the GKO algorithm, we now let Φ^{(k)} and Ψ^{(k)} be the computed values of these quantities, u_{k:}, r_{k:n,k}^{(k)}, l_{:k}, Φ^{(k+1)} and Ψ^{(k+1)} be the values of these quantities computed in exact arithmetic from Φ^{(k)} and Ψ^{(k)} using steps 1 to 3 of Algorithm 2.1, and ũ_{k:}, r̃_{k:n,k}^{(k)}, l̃_{:k}, Φ̃^{(k+1)} and Ψ̃^{(k+1)} be the actual computed values of

u_{k:}, r_{k:n,k}^{(k)}, l_{:k}, Φ^{(k+1)} and Ψ^{(k+1)} respectively. The methodology is as follows:

1. Using a standard rounding error analysis, derive expressions of the form

ũ_{k:} = u_{k:} + δũ_{k:}  (23)
r̃_{k:n,k}^{(k)} = r_{k:n,k}^{(k)} + δr̃_{k:n,k}^{(k)}  (24)
l̃_{:k} = l_{:k} + δl̃_{:k}  (25)
Φ̃^{(k+1)} = Φ^{(k)} − l̃_{:k}ϕ_{k}^{(k)} + δΦ̃^{(k+1)}  (26)
Ψ̃^{(k+1)} = Ψ^{(k)} − ψ_{k}^{(k)}ũ_{k:}/r̃_{kk}^{(k)} + δΨ̃^{(k+1)}  (27)

where δũ_{k:}, etc. are error terms.

2. Evaluate Φ^{(k)}Ψ^{(k)} − Φ̃^{(k+1)}Ψ̃^{(k+1)} using (23) to (27). This can be expressed in the form

Φ^{(k)}Ψ^{(k)} − Φ̃^{(k+1)}Ψ̃^{(k+1)} = A_{f}l̃_{:k}ũ_{k:} − l̃_{:k}ũ_{k:}A_{b} + F_{k} ,  (28)

where F_{k} is an error term. By (21),

H_{k} = −F_{k} .

3. After some manipulation, express F_{k} as a sum of terms of the form

S(A_{f}, A_{b}) ∘ T (V^{(k)}) ∘ l̃_{:k}ũ_{k:} ∘ Δ̂ or S(A_{f}, A_{b}) ∘ T (V^{(k+1)}) ∘ L_{:,k+1:n}U_{k+1:n,:} ∘ Δ̂ .

Here, the S(A_{f}, A_{b}) are matrices formed from A_{f} and A_{b}, Δ is a matrix whose elements are bounded in magnitude by ϵ, and V^{(k)} is defined by

|Φ^{(k)}||Ψ^{(k)}| ≡ V^{(k)} ∘ Φ^{(k)}Ψ^{(k)} .  (29)

4. Apply (19) to derive an expression for G.

5. Lemma 3.1 shows that G satisfies

∇_{{A_{f} ,A_{b}}}E = G .  (30)

Using the appropriate algorithm to recover a structured matrix from its generators, derive an expression for E from the expression for G. Note that in general, G will be of full rank. However, (30) will still be satisfied by E and G.

6. Derive bounds for ∥E∥ using some norm.

3.3 Error analysis of GKO for Cauchy-type matrices

In this subsection, we use the above methodology to derive the first of our main results — a bound for ∥E∥ when a Cauchy matrix R is factorized by the GKO algorithm. The results are encapsulated in three theorems, which yield expressions for the {H_{k}}, an elementwise bound for G, and a bound for ∥E∥ respectively. We then discuss the size of the bound for ∥E∥.

Theorem 3.4 Let R be a Cauchy matrix to be factorized by the GKO algorithm and let F_{k}, H_{k}, V^{(k)}, l̃_{:k}, ũ_{k:} be as defined above. Then

F_{k} = c_{1}Δ̂^{(1)} ∘ D_{vc}^{(k)}D_{p}l̃_{:k}ũ_{k:} + c_{2}l̃_{:k}ũ_{k:}D_{q}D_{vr}^{(k)} ∘ Δ̂^{(2)} + c_{3}(r_{kk}^{(k)})^{−1}v_{kk}^{(k)}Δ̂^{(3)} ∘ l̃_{:k}ũ_{k:} +
c_{4}Δ̂^{(4)} ∘ B^{I} ∘ V^{(k+1)} ∘ L̃_{:,k+1:n}Ũ_{k+1:n,:}

where c_{1} to c_{4} are small constants, D_{vc}^{(k)} = diag(v_{:k}^{(k)}), D_{vr}^{(k)} = diag(v_{k:}^{(k)}), D_{p} = diag{t_{i} − s_{k}}_{i}, D_{q} = diag{t_{k} − s_{j}}_{j}, B = [1/(t_{i} − s_{j})] is the ordinary Cauchy matrix with displacement operator ∇_{{D_{s},D_{t}}}, the Δ̂^{(·)} are matrices whose elements are less than ϵ in magnitude, and H_{k} = −F_{k}.

Proof. In the following, we simplify our notation and drop the superscript (k); where the superscript is (k + 1) we indicate this by a prime (′); and we drop the subscripts : k, k : and k : n, k. In the following, we do not give all the steps in the derivation of the various expressions, as these are straightforward but very tedious. However, we indicate how key intermediate expressions are derived.

We use the normal properties of floating point operations performed with at least one guard digit, viz. fl(a) = a(1+δ_{1}) and fl(a⋆b) = (a⋆b)(1+δ_{2}), where fl(a) denotes rounding, fl(a⋆b) is the computed result of any of the four basic floating-point operations, and |δ_{1}|, |δ_{2}| < ϵ.

Following step 1 of the above methodology, we evaluate expressions for the computed values of r̃, l̃ and ũ (subscripts and superscripts dropped), yielding after a few steps

ũ = u + 2ũΔ^{(1)} + ϕ_{k}D^{(1)}ΨD_{q}^{−1} ,  (31)
r̃ = r + 2Δ^{(2)}r̃ + D_{p}^{−1}D^{(2)}Φψ_{k} ,  (32)
l̃ = l + 5Δ^{(3)}l̃ + r̃_{kk}^{−1}D_{p}^{−1}D^{(2)}Φψ_{k} − b_{kk}r̃_{kk}^{−1}∂_{k}^{(2)}ϕ_{k}ψ_{k}l̃ .

Here and below the Δ^{(·)} denote diagonal matrices with elements of magnitude less than ϵ; the D^{(·)} are elementwise operators which multiply each element of their matrix operands by a factor less than ϵ, and the ∂_{k}^{(·)} are similar elementwise vector operators.

Similarly, it can be shown that the computed values of Φ′ and Ψ′ satisfy

Φ̃′ = Φ − l̃ϕ_{k} + D^{(3)}Φ′ + D^{(4)}(l̃ϕ_{k}) ,
Ψ̃′ = Ψ − ψ_{k}ũ/r̃_{kk} + D^{(5)}Ψ′ + 2D^{(6)}(ψ_{k}ũ)/r̃_{kk} .

Carrying out step 2 of the above methodology, we obtain

ΦΨ − Φ̃′Ψ̃′ = Φψ_{k}ũ/r̃_{kk} + l̃ϕ_{k}Ψ − l̃ϕ_{k}ψ_{k}ũ/r̃_{kk} − 2ΦD^{(6)}(ψ_{k}ũ)/r̃_{kk} − Φ′D^{(5)}Ψ′ −
D^{(4)}(l̃ϕ_{k})Ψ − D^{(3)}Φ′Ψ′ + D^{(4)}(l̃ϕ_{k})ψ_{k}ũ/r̃_{kk} + 2l̃ϕ_{k}D^{(6)}(ψ_{k}ũ)/r̃_{kk} . (33)

Let T_{3} denote the first three terms in (33). From Algorithm 2.2, we have Φψ_{k} = D_{p}r and ϕ_{k}Ψ = ũD_{q}. Using these relations in T_{3}, and expressing r in terms of (r̃ - error terms) using (32) and u in terms of (ũ - error terms) using (31), we can show that

T_{3} = D_{t}l̃ũ − l̃ũD_{s} − 3D_{p}Δ^{(4)}l̃ũ − 2l̃ũΔ^{(5)}D_{q} + 2r_{kk}^{−1}δl̃ũ −
D^{(2)}Φψ_{k}ũ/r_{kk} − l̃ϕ_{k}D^{(1)}Ψ + r̃_{kk}^{−1}∂_{k}^{(2)}ϕ_{k}ψ_{k}l̃ũ ,  (34)

where |δ| < ϵ. By using (34) for the first three terms of (33), we get an equation of the form (28), where F_{k} is given by the last six terms in (33) plus the last six terms in (34). Terms involving the D^{(·)} may be expressed in terms of l̃ũ or L̃Ũ by using the definition of V, which in the current notation is

v_{ij} = \frac{|ϕ_{i}||ψ_{j}|}{ϕ_{i}ψ_{j}} .

Consider the factor ΦD^{(6)}(ψ_{k}ũ)/r̃_{kk} in the term −2ΦD^{(6)}(ψ_{k}ũ)/r̃_{kk}. We have

(ΦD^{(6)}(ψ_{k}ũ)/r̃_{kk})_{ij} = ϕ_{i}∂_{j}^{(6)}(ψ_{k}ũ_{j})/r̃_{kk} .

Recall that ϕ_{i} = [ϕ_{i1}, ϕ_{i2}] and ψ_{j} = [ψ_{1j}, ψ_{2j}]. Then

(ΦD^{(6)}(ψ_{k}ũ)/r̃_{kk})_{ij} = (ϕ_{i1}δ_{1j}^{(6)}ψ_{1k} + ϕ_{i2}δ_{2j}^{(6)}ψ_{2k})ũ_{j}/r̃_{kk} ,

where δ_{1j}^{(6)} and δ_{1j}^{(6)} are the scaling factors from the operator ∂_{j}^{(6)}. From the definition of V, using the fact that l̃_{i} ≐ r̃_{ik}/r̃_{kk}, this can be shown to be

(ΦD^{(6)}(ψ_{k}ũ)/r̃_{kk})_{ij} = δ̂_{ij}^{(6)}v_{ik}b_{ik}^{−1}l̃_{i}ũ_{j} ,

where |δ̂_{ij}^{(6)}| ≤ max_{j=1,2} |δ_{kj}^{(6)}|. In matrix form, we obtain

ΦD^{(6)}(ϕ_{k}ũ)/r̃_{kk} = Δ̂ ∘ diag{v_{ik}/b_{ik}}l̃ũ ,

where Δ̂ and subsequent Δ̂^{(·)} are matrices with elements bounded in magnitude by ϵ. Similarly, all the other terms can be expressed as either

(i) an elementwise product of Δ̂^{(·)} and a normal product of l̃ũ and matrices derived from B or V, or

(ii) an elementwise product of the form Δ̂^{(·)} ∘ B^{I} ∘ V′ ∘ L_{:,k+1:n}U_{k+1:n,:}.

When this is done, the result follows.

The next theorem uses Lemma 3.3 to obtain an elementwise bound for |G|.

Theorem 3.5 Let H_{k} be as in Theorem 3.4. Then

|G| ≤ c_{1}b_{min}^{−1}Δ̂^{(1)} ∘ |L̂||U| + c_{2}b_{min}^{−1}|L||Û| ∘ Δ̂^{(2)} + c_{3}b_{min}^{−1}Δ̂^{(3)} ∘ |L|diag{v_{kk}^{(k)}}|U|
+ c_{4}|B^{I}| ∘ Δ̂^{(4)} ∘ ∑_{k=2}^{n} |R̂′_{k}|

where b_{min} is the minimum modulus of the elements of B, L̂ = [v_{:k}^{(k)}]_{k=1}^{n} ∘ L, Û = U ∘ [v_{k:}^{(k)}]_{k=1}^{n}, and R̂′_{k} = V^{(k)} ∘ L_{:,k:n}U_{k:n,:}.

Proof. G is evaluated by carrying out the summation in (19), and using the identities ∑_{i=k}^{n} a_{:k}b_{k:} = AB and ∑_{i=k}^{n} x_{k}a_{:k}b_{k:} = A diag{x_{k}}B.

We now apply the last step in the above methodology to derive an expression for ∥E∥.

Theorem 3.6 Let E be the backward error E = L̃Ũ − R in the factorization of R using the GKO algorithm, let L̂, Û, R̂, B and V be as above. Then ∥E∥ is bounded by

∥E∥ ≤ ϵ (c_{5} \frac{b_{max}}{b_{min}} g_{1} + c_{6}ng_{2}) ∥L∥∥U∥ ,  (35)

where the Frobenius norm is used, b_{max} and b_{min} are the maximum and minimum moduli of the elements of B, c_{5} and c_{6} are small constants, and g_{1} and g_{2} are generator growth factors, defined by

g_{1} = c_{7} \frac{∥L̂∥}{∥L∥} + c_{8} \frac{∥Û∥}{∥U∥} + c_{9}∥diag\{v_{kk}^{(k)}\}∥ ,  (36)
g_{2} = max_{k=2,...,n}{|R̂_{k}∥/∥R_{k}∥} ,  (37)

with c_{7}, c_{8}, c_{9} < 1.

Proof. From step 5 of the above methodology, we essentially invert the Sylvester equation (30) to derive an expression for E. To do this we apply (18) in Corollary 3.2. This can be written in matrix form

E = B ∘ G

so

|E| = |B| ∘ |G| ≤ c_{1} \frac{b_{max}}{b_{min}} Δ̂^{(5)} ∘ |L̂||U| + c_{2} \frac{b_{max}}{b_{min}} ∥L||Û| ∘ Δ̂^{(6)} +
c_{3} \frac{b_{max}}{b_{min}} Δ̂^{(7)} ∘ |L|diag\{v_{kk}\}|U| + c_{4}Δ^{(8)} ∘ ∑_{k=2}^{n} |R̂′_{k}| .  (38)

We now define g_{2} ≡ max_{k=2,...,n} ∥R̂^{(k)}∥/∥R^{(k)}∥, g_{4} ≡ ∥L̂∥/∥L∥, g_{5} ≡ ∥Û∥/∥U∥ and g_{6} ≡ ∥diag{v_{kk}^{(k)}}∥. These can be considered to be generator growth factors — they are functions of the V^{(k)}, which from the definition (29) are the ratio of the products of the magnitudes of the generators to the products of the generators. We will see in §5 that these growth factors can sometimes be large.

Taking the Frobenius norm of (38), we can easily show that

∥E∥ ≤ c_{1}δ_{1} \frac{b_{max}}{b_{min}} g_{3}∥L∥∥U∥ + c_{2}δ_{2}g_{4} \frac{b_{max}}{b_{min}} ∥L∥∥U∥ + c_{3}δ_{3} \frac{b_{max}}{b_{min}} g_{5}∥L∥∥U∥ + c_{6}nδ_{4}g_{2}∥L∥∥U∥ . (39)

where 0 ≤ |δ_{1}|, …, |δ_{4}| < ϵ. The result follows by collecting the first three terms of (39).

The following corollary specializes the above result to the case when R is derived from a Toeplitz matrix.

Corollary 3.7 Let R be derived from a Toeplitz matrix T by the transformation (10) in Theorem 2.2, and let c_{1}, c_{2}, g_{1}, g_{2} and E be as defined in Theorem 3.6. Then ∥E∥ is bounded by

∥E∥ ≤ ϵc_{10}g_{3}n∥L∥∥U∥ ,  (40)

where c_{10} = max(2c_{5}/π, c_{6}) and g_{3} = max(g_{1}, g_{2}).

Proof. Recall that B = [1/(t_{i}−s_{j})] is the ordinary Cauchy matrix with displacement operator ∇_{{D_{s},D_{t}}}; from equations (11) in Theorem 2.2, the t_{i} are n equally-spaced points around the unit circle, including one at (1,0), and the s_{j} are also n equally-spaced points around the unit circle, with each s_{j} between two of the t_{i}. Clearly π/n < t_{i} − s_{j} < 2 ∀i, j, so by the definition of B,

\frac{b_{max}}{b_{min}} < 2n/π .  (41)

Substituting (41) in (39), bounding 2c_{5}/π and c_{6} by c_{10}, and bounding g_{1} and g_{2} by g_{3} yields the result.

The above results show that the expressions for the backward error bounds from the GKO algorithm are similar to the ones for Gaussian elimination with partial pivoting (GE/PP) [9], except for the generator growth factors which might arise in particular cases where the Φ^{(k)} and Ψ^{(k)} are large, but not the Φ^{(k)}Ψ^{(k)} or the R_{k}. So there may be some cases where large error growth may occur in the GKO algorithm but not GE/PP. In §5, we give an example where this occurs.

4 Error Analysis of the GKO-Toeplitz Algorithm

Recall that the steps in the GKO-Toeplitz algorithm are (i) compute the generators from the Toeplitz matrix T using (4) and (5), (ii) convert them to generators of a Cauchy matrix using (12) and (iii) compute factors L and U of this Cauchy matrix using the GKO algorithm. The factors of T are then given by (13). There are errors incurred at each of these steps. In this section, we do not consider permutations, as these do not contribute to the error. We will derive a bound for the perturbation matrix E_{T}, defined by

F^{∗}L̃ŨF D = T + E_{T} .

In our development, we show in Theorem 4.1 that E_{T} consists of two components — the first due to the error ∥E∥ incurred in the Cauchy factorization and the second due to the errors incurred in computing the Cauchy generators Φ̃ and Ψ̃. The latter is a Toeplitz-type perturbation ΔT such that T + ΔT transforms exactly to Φ̃ and Ψ̃. We then derive two lemmas needed to derive ΔT, and then present the main result of this section in Theorem 4.4.

4.1 Main components of E_{T}

E_{T} has two main components, as is shown in the following.

Theorem 4.1 Let F and D be as in Theorem 2.2, let Φ̃ and Ψ̃ be the Cauchy generators computed using (4), (5) and (12), and let L̃ and Ũ be the factors computed from Φ̃ and Ψ̃ using the GKO algorithm. Then the perturbed factorization of T satisfies

F^{∗}L̃ŨF D ≡ T + E_{T} = T − F^{∗}EF D + ΔT ,  (42)

where E is as in Theorem 3.6 and ΔT is a Toeplitz-type perturbation of T such that T + ΔT has generators Ω̃ and Γ̃ that transform exactly to Φ̃ and Ψ̃ using (4), (5) and (12).

Proof. Let R̃ be the Cauchy matrix generated by Φ̃ and Ψ̃. We have

R̃ = L̃Ũ + E ,

and we know from (10) that Φ̃ and Ψ̃ are the generators for

R̃ = F (T + ΔT)D^{−1}F^{∗}

where T + ΔT is some Toeplitz-type matrix. From the above two equations we obtain

T + ΔT = F^{∗}R̃F D = F^{∗}(L̃Ũ + E)F D ,

from which the desired result follows.

Thus, by (42), we see that E_{T} has one component with the same norm bound as E, and another which perturbs T to a matrix such that its generators, say Ω̃ and Γ̃, transform exactly to Φ̃ and Ψ̃. Before we derive an expression for ΔT, we need two preliminary results : expressions for Ω̃ and Γ̃, and a method to recover T + ΔT from its generators Ω̃ and Γ̃.

4.2 Estimation of ΔT — preliminary results

The required results are given in the following two lemmas.

Lemma 4.2 Let Ω and Γ be as in (4) and (5), and let Ω̃ and Γ̃ transform exactly to Φ̃ and Ψ̃ using (12). Let [a, b] = Ω̃ − Ω and let [c, d ] = Γ̃^{∗} − Γ^{∗}. Then

a = 0

and ∥b∥, ∥c∥ and ∥d∥ are bounded by

∥b∥ ≤ ϵk_{1}n^{3/2}∥ω_{:2}∥ ,  (43)
∥c∥ ≤ ϵk_{2}n^{3/2}∥γ_{1:}∥ ,  (44)
∥d∥ ≤ ϵ .  (45)

Proof. We first consider the errors incurred in the computation of Φ̃ and Ψ̃. We have

Φ̃ = fl{F̃[e_{1}, ω̃_{:2}]}, where F̃ = fl(F) , ω̃_{:2} = fl(ω_{:2})
= [1, fl(F̃ω̃_{:2})] , where 1 = [1, 1, …, 1]^{T}
= [1, F̃ω̃_{:2} + k_{3}n∥ω̃_{:2}∥δ^{(1)}]

where |δ_{i}^{(1)}| < ϵ , i = 1, …, n. After a few more steps, this becomes

Φ̃ = F [e_{1}, ω_{:2} + b]

where b = Δ^{(7)}ω_{:2} + k_{4}(n + 1)∥ω̃_{:2}∥F^{∗}δ^{(1)}. In a similar way, it can be shown that

Ψ̃^{∗} = F D[γ_{1:}^{∗} + c, e_{n} + d]

where c = k_{5}D^{∗}Δ^{(8)}Dγ_{1:}^{∗} + k_{6}(n + 1)∥γ_{1:}∥δ^{(2)} and d = D^{∗}F^{∗}d_{n}Δ^{(9)}f_{n:}^{T}. Now the expressions in square brackets transform exactly to Ω̃ and Γ̃ respectively, and by taking norms of b, c and d the bounds (43) to (45) can be demonstrated in a few steps.

Lemma 4.3 For any matrix A, let ∇_{{Z_{1},Z_{−1}}}A = B. Then A can be recovered from B using

a_{ij} = ∑_{k=j}^{n} b_{1+(i+k−j) mod n,k} − ∑_{k=1}^{j−1} b_{1+(i+k−j) mod n,k} .  (46)

Proof. From the displacement operator ∇_{{Z_{1},Z_{−1}}}, the following properties of B are easily seen:

b_{ij} = a_{i−1,j} − a_{i,j−1}, 1 < i ≤ n , 1 ≤ j < n ,  (47)
b_{1j} = a_{nj} − a_{i,j+1}, 1 ≤ j < n ,  (48)
b_{in} = a_{i−1,j} + a_{i1}, 1 < i ≤ n and  (49)
b_{1n} = a_{n,n−1} + a_{11} .  (50)

It can be verified that if the elements of A are given by (46), then (47) to (50) are satisfied.

Equation (46) shows that an element a_{ij} is recovered by computing x − y, where x is the sum of elements of B down the diagonal, commencing from b_{i+1,j} and proceeding to the last column, wrapping from the last row to the first if necessary during the summing; y is a similar “wrapped diagonal sum” from the first column to b_{i,j−1}.

4.3 Main result

We now use Theorem 4.1, Lemma 4.2 and Lemma 4.3 to derive a bound for the backward error ∥E_{T} ∥ in the GKO-Toeplitz algorithm.

Theorem 4.4 Let F and D be as in Theorem 2.2, and let L̃ and Ũ be the factors computed from T using the GKO-Toeplitz algorithm. Then the perturbed factorization of T satisfies

F^{∗}L̃ŨF D ≡ T + E_{T} = T + E^{(1)} + E^{(2)} ,  (51)

where E^{(1)} is a general matrix with norm ∥E^{(1)}∥ = ∥E∥, E is as in Theorem 3.6, and E^{(2)} is a Toeplitz-type matrix with norm bounded by

∥E^{(2)}∥ ≤ ϵc_{11}n^{2}(∥t_{1:}∥ + ∥t_{:1}∥) .  (52)

Proof. By comparing (51) and (42), we see that E^{(1)} = −F^{∗}EF D, and because F and D are orthogonal matrices,

∥E^{(1)}∥ = ∥E∥ .  (53)

From the above comparison we also have E^{(2)} = ΔT, a Toeplitz-type perturbation of T such that T + ΔT has generators Ω̃ and Γ̃ that transform exactly to the Cauchy generators Φ̃ and Ψ̃ computed using (4). In the following, we use E^{(2)} for ΔT. From Lemma 4.2, we have

∇(T + E^{(2)}) = Ω̃Γ̃ = ΩΓ + e_{1}c^{∗} + ω_{:2}d^{∗} + be_{n}^{T} ,

where b, c and d are bounded as in (43) to (45). The second-order error term bd^{∗} has been omitted. We then have

∇E^{(2)} = e_{1}c^{∗} + ω_{:2}d^{∗} + be_{n}^{T} ,

and we use (46) to compute E^{(2)}. This yields, after some algebra

|e_{:j}^{(2)}| = C_{j−1}(|c^{R}| + |b^{R}|) + |p_{:j}|

where C_{k} is a matrix which by premultiplication, circularly upshifts a vector k places, x^{R} indicates the reversal of x, and the moduli of p_{:j} are bounded by

|p_{ij}| ≤ |ω_{:2}|^{T} C_{j−i−1}|d|
≤ ∥ω_{:2}∥∥d∥ .  (54)

Using (43), (44), (54) and (45) it is easily seen that

∥e_{:j}^{(2)}∥ ≤ c_{12}n^{3/2}(∥ω_{:2}∥ + ∥γ_{1:}∥) .

From this, using the definitions (4) and (5), we obtain the bound (52) for E^{(2)}. Together with (53), this yields the result.

5 Discussion of Error Bounds

We first discuss the factors in the above error bounds and relate them to what would be expected for Gaussian elimination with partial pivoting (GE/PP). Then we show, for both the Cauchy and Toeplitz variants, that there are some cases where the backward error growth can be large.

5.1 Relation of bounds to those for GE/PP

Consider the backward error E incurred by the Cauchy variant (equation (35)). The term ∥L∥∥U∥ is similar to that obtained for GE/PP [9]. However, the first factor contains the generator growth factors g_{1} and g_{2}. These are given by ratios of norms of the hatted quantities to the unhatted quantities in (36) and (37). The former are derived from the latter by elementwise multiplication by submatrices of the V^{(k)}, which from their definitions (29) are the ratio of the products of the magnitudes of the generators to the products of the generators. For an ordinary Cauchy matrix, v_{ij}^{(k)} = 1 ∀i, j, k because Φ^{(k)} and Ψ^{(k)} have only one column and row respectively. However, for higher displacement-rank Cauchy matrices, there may be significant cancellation in the computation of the denominator of (29), so they may be significant growth in the size of the L̂, Û and R̂_{k} compared to the L, U and R_{k} respectively.

The backward error E_{T} incurred by the Toeplitz variant has two components — one with the same norm as E above, and a Toeplitz-type component with norm bounded as in (52). The latter bound is proportional to n^{2} and contains no growth factors, so it would be expected that the bound would be dominated by the first component.

We next give examples where the generator growth might be expected to be large in the Cauchy and Toeplitz variants.

5.2 Examples of large generator growth

Cauchy case. Here, we can select an example where all the elements of V = V^{(1)} are large. This will occur when significant cancellation occurs in the computation of the ϕ_{i}ψ_{j}. Such an example is

Φ = [a, a + f] , Ψ = [a, −a]^{T} ,

where ∥a∥ is of order unity, and ∥f∥ is very small. Then ΦΨ = −fa, that is, all the elements of ΦΨ are very small compared to those of |Φ||Ψ|. Moreover, because a and f can be arbitrary except for their norms, the original matrix [(t_{i} − s_{j})^{−1}ϕ_{i}ψ_{j}] is in general well-conditioned.

Toeplitz case. The Toeplitz case has an extra constraint on the selection of Φ and Ψ, since it must be generated from Ω and Γ using the transformations (12). Because of this constraint, there is no case where all the elements of V can be made large. However, all of the first column of V can be made large, and this will cause error growth, in spite of the pivoting. This will happen in the following case.

Recall that a_{i−j} = t_{ij} ∀i, j. Select

a_{0} = 1  (55)
and a_{i} = −a_{i−n} , 1 ≤ i ≤ n − 1 ,  (56)

so that Ω = [e_{1}, e_{1}]. Then all of the first column of V will be large if ψ_{11} + ψ_{12} is very small compared to ψ_{11} and ψ_{12}. It can be verified from (12) that if we select a_{1}, …, a_{n−1} to satisfy

∑_{j=1}^{n−1} a_{n−j} exp(iπ(j − 1)/n) = − exp(iπ(n − 1)/n + δ/2  (57)

then ψ_{11} + ψ_{12} = δ. There is a wide variety of choices for the a_{j}. Let n be even, and set

a_{n−j} = 0  (58)

except for a_{n/2−1} and a_{n−1}. Then (57) is satisfied when

a_{n/2−1} = − sin(π/n) + ℑ(δ/2) , a_{n−1} = cos(π/n) + ℜ(δ/2) .  (59)

So if δ is small, and the a_{j} are selected according to (55), (58), (59) and (56), all of the first column of V will be large, with magnitude O(1/δ).

Numerical examples. Order-8 Toeplitz matrices were generated according to (55), (58), (59) and (56), with δ = 10^{−k}, k = 2, …, 16. For each matrix, the system T x = 1 was solved. It was found that the normalized solution error ∥x̃ − x∥/∥x∥ grew as the square of 1/δ, and the normalized residual ∥T x̃ − 1∥/∥b∥ grew linearly with 1/δ. Thus the algorithm is only weakly stable in this case.

6 Modified GKO Algorithm

The problem with the original pivoting strategy is that when all elements of r_{:1} are small and all elements of v_{:1} are large, normal partial pivoting will not stabilize the algorithm. Complete pivoting will do so, but requires O(n^{2}) operations to find the pivot at each major step and O(n^{3}) operations overall. However, a strategy of using the largest element in the first row and column should stabilize the algorithm in most cases, and we see that it does in the above cases.

To incorporate this row-1/column-1 pivoting, it is easy to see that the following steps should be added to the GKO algorithm (Algorithm 2.2):

• Step 1: add substep P′ ← I, where P′ will be the matrix of column interchanges.

• After loop to recover column 1 of R_{k} : add loop to recover row 1 of R_{k}.

• After loop to find maximum max_{1} in column 1 : add loop to find maximum max_{2} in row 1. If max_{1} ≥ max_{2}, carry out row interchanges as in Algorithm 2.2. Otherwise carry out column interchanges by swapping the appropriate elements in s, Ψ^{(k)} and r_{k:}^{(k)}, and the appropriate columns in U and P′.

• After computation of L, U, P and P′, the factors of R are P^{T} LUP^{′T}.

Results. When the modified algorithm was used on the same set of systems as was considered in the previous section, it was found that the normalized solution error ∥x̃ − x∥/∥x∥ grew linearly with 1/δ and the condition number of T, and the normalized residual ∥T x̃−1∥/∥b∥ was approximately constant at about 4 × 10^{−15}, a small multiple of ϵ. Thus the modified algorithm is stable in this case.

7 Conclusions

It has been shown that bound for the backward error in the GKO algorithm is similar to that for partial pivoting, except that extra factors, the generator growth factors, are included. These factors can be large when there is sufficient cancellation in the computation of the generators. Examples of this have been presented, and it was demonstrated that the original GKO algorithm was only weakly stable in these cases. A modified version which uses row 1/column 1 pivoting was then presented; this version was stable in these cases.

It is not known whether there are any cases upon which the modified algorithm will give large errors. Further work needs to be done to ascertain this, and if such cases can be found, the pivot strategy needs to be improved further. The aim is to find the maximum in R, or an element close to the maximum, still in O(n) operations. An extension of the above strategy may be to have a few iterations in the search, i.e. search for the row-1/column-1 maximum, say at r_{1p}, then search along column p for the maximum there, and so on. This may find a better pivot at the expense of some extra work.

A practical strategy is to use the modified algorithm of §6 followed by a check of the residual; in the unlikely event that the residual is large we can resort to a stable O(n^{3}) algorithm.

References

[1] E. Bareiss, “Numerical solution of linear equations with Toeplitz and vector Toeplitz matrices”, Numer. Math. 13 (1969), 404–424.

[2] A. W. Bojanczyk, R. P. Brent and F. R. de Hoog, “Stability analysis of a general Toeplitz system solver”, Numerical Algorithms, to appear. Preliminary version available as TR-CS-93-15, CSL, ANU, August 1993 (revised June 1994) [available by anonymous ftp from

nimbus.anu.edu.au:/pub/Brent/rpb143tr.*].

[3] A. W. Bojanczyk, R. P. Brent, F. R. de Hoog and D. R. Sweet, “On the stability of the Bareiss and related Toeplitz factorization algorithms”, SIAM J. Matrix Analysis Appl. 16 (1995), 40–57.

[4] R. Brent, “Old and new algorithms for Toeplitz systems”, Proceedings SPIE, Volume 975, Advanced Algorithms and Architectures for Signal Processing III (edited by Franklin T. Luk), SPIE, Bellingham, Washington, 1989, 2–9.

[5] J. R. Bunch, “Stability of methods for solving Toeplitz systems of equations”, SIAM J. Sci. Stat. Comp. 6 (1985), 349–364.

[6] T. F. Chan and P. C. Hansen, “A lookahead Levinson algorithm for general Toeplitz systems”, IEEE Proc. Signal Processing 40 (1992), 1079–1090.

[7] J. Chun and T. Kailath, “Fast triangularization and orthogonalization of Hankel and Vandermonde matrices”, Linear Alg. Apps. 151 (1991), 199–228.

[8] I. Gohberg, T. Kailath and V. Olshevsky, “Gaussian elimination with partial pivoting for structured matrices”, preprint, May 1994.

[9] G. H. Golub and C. Van Loan, Matrix Computations, 2nd ed., Johns Hopkins Press, 1989.

[10] M. H. Gutknecht and M. Hochbruck, “Look-ahead Levinson and Schur algorithms for non-Hermitian Toeplitz systems”, IPS Research Rept. 93-11, ETH-Zürich, August 1993.

[11] D. R. Sweet, “The use of pivoting to improve the numerical performance of Toeplitz matrix algorithms”, SIAM J. Matrix Anal. Appl. 14 (1993), 468–493.