Stability analysis of fast Toeplitz linear system solvers^{∗}

A. W. Bojanczyk School of Electrical Engineering Cornell University Ithaca, NY 14853-5401

R. P. Brent^{†}

Computer Sciences Laboratory Australian National University Canberra, ACT 2601

F. R. de Hoog Division of Mathematics and Statistics CSIRO, Canberra, ACT 2601

Abstract We present a numerical stability analysis of the Bareiss and Levinson algorithms for

solving a symmetric positive definite Toeplitz system of linear equations.

0 Introduction

Our interest is in analysis of numerical stability of algorithms for solving a linear system

T x = b  (0.1)

where T is an n × n positive definite Toeplitz matrix and b is an n × 1 vector. We solve this system in floating point arithmetic with relative precision ϵ.

We adopt from [14] the following defintions of numerical stability and good-behavior.

Definition 0.1: We say that an algorithm for solving the equation (0.1) is numerically stable if the computed solution x̃ satisfies

||x − x̃|| ≤ c_{1}ϵcond(T)||x̃||  (0.2)

where the constant c_{1} may depend on the dimension n of the system, and cond(T) = ||T || ||T^{−1}|| is the condition number of T.

Definition 0.2: We say that an algorithm for solving the equation (0.1) is well-behaved if the computed solution x̃ satisfies

||T x̃ − b|| ≤ c_{2}ϵ||T || ||x̃||  (0.3)

where the constant c_{2} may depend on the dimension n of the system.

The concepts of numerical stability and good behavior are related to those of weak and strong stability as presented in [7]. The reader is referred to [17] for in depth treatment of roundoff analysis of matrix algorithms.

It is known that an algorithm is well-behaved iff there exists a matrix ΔT such that

(T + ΔT)x̃ = b , ||ΔT || ≤ c_{3}||T || .  (0.4)

∗Appeared as Report CMA-MR17-91, CMA, ANU, August 1991. Copyright \begin{matrix}c \\ ◯ \end{matrix} 1991, the authors.

†E-mail address: rpb@cslab.anu.edu.au rpb126tr typeset using LaTEX

Note that in our definitions we do not require that for a Toeplitz matrix T the perturbation ΔT be also Toeplitz.

It is also well known that good-behavior implies numerical stability but not vice versa. This is manifested by the size of the residual vector. Most direct methods for solving a general system of linear equations are well-behaved while most iterative methods are only numerically stable [14].

The roundoff error analysis of Toeplitz systems solvers was presented by Cybenko in [9], and by Sweet in [18]. Cybenko showed that Levinson-Durbin algorithm produces the residual which, under the condition that all so called reflection coefficients are positive, is of comparable size to that produced by the well behaved Cholesky method. Next he hypothesised that the same is true even if the reflection coefficients are not all positive. This would indicate that numerical quality of the Levinson-Durbin algorithm is comparable to that of the Cholesky method.

In his PhD thesis Sweet [18] presented a roundoff error analysis of a variant of Bareiss method [2] and concluded that the method is well-behaved. In this paper we generalize both earlier results on stability of Levinson, and Bareiss algorithms.

We show that that a class of Bareiss type algorithms for solving a positive definite Toeplitz system of linear equations is well-behaved. We also show that Levinson algorithm is numerically stable. Numerical experiments suggest that a standard implementation of Levinson algorithm is not well behaved [19].

Our approach is the following. We treat Bareiss and Levinson algorithms as matrix decomposition algorithms and derive them as sequences of downdating operations [3]. Next we present a first order analysis by bounding the first term in an asymptotic expansion for the error in powers of ϵ, [17]. By analyzing the propagation of the error of the first order in the sequence of downdatings that define the algorithms we obtain bounds on the perturbations of the factors in the decompositions. The bounds conform with the Definition 1.1 and Definition 1.2 for Levinson and Bareiss algorithms respectively.

1 Notation.

Unless it is clear from the context, all vectors will be real and of dimension n. Likewise, all matrices are real and their default dimension is n × n. For a ∈ ℜ^{n}, ∥a∥ denotes the usual Euclidean norm and if T ∈ ℜ^{n×n},

∥T ∥ = max_{a∈ℜ^{n}} \frac{∥T a∥}{∥a∥}.

Our primary interest is in symmetric positive definite Toeplitz matrices T whose i, jth entry is

t_{ij} = t_{|i−j|}.

We denote by e_{k}, k = 1, …, n the unit vector whose elements, with the exception of the kth element which is one, are all zero.

Z = ∑_{k=1}^{n−1} e_{k+1}e_{k}^{T}

J = ∑_{k=0}^{n−1} e_{k+1}e_{n−k}^{T}
Z_{0} = I

Z_{k} = Z^{k}
J_{k} = JZ_{n−k}
E_{k} = Z_{n−k}^{T}Z_{n−k}

The matrix H(θ) ∈ ℜ^{2×2} is defined by

H(θ) = \frac{1}{cos θ} [ \begin{matrix}1 & − sin θ \\ − sin θ & 1 \end{matrix} ]  (1.1)

and we note that it has eigenvalues

λ_{1}(θ) = λ_{2}^{−1}(θ) = sec θ − tan θ.

2 Elementary Downdating.

In this Section we introduce the concept of elementary downdating. In the subsequent sections Bareiss and Levinson algorithms will be derived in terms of a sequence of downdating steps. The numerical properties of the two algorithms will then be related to the properties of the sequence of elementary downdating steps.

Suppose that we have u_{k}, v_{k} ∈ ℜ^{n} that satisfy

e_{j}^{T} u_{k} = 0 , j < k  (2.1a)
e_{j}^{T} v_{k} = 0 , j ≤ k  (2.1b)

and we wish to find u_{k+1}, v_{k+1} ∈ ℜ^{n} that satisfy

u_{k+1}u_{k+1}^{T} − v_{k+1}v_{k+1}^{T} = Zu_{k}u_{k}^{T} Z^{T} − v_{k}v_{k}^{T}  (2.2a)

and

e_{j}^{T} u_{k+1} = 0 , j < k + 1  (2.2b)
e_{j}^{T} v_{k+1} = 0 , j ≤ k + 1  (2.2c)

We refer to this problem as the elementary downdating problem. This problem is a special subproblem of a more general problem that arises in Cholesky factorization of a positive definite difference of two outer product matrices, see [6], [5], [1].

Clearly, the elementary downdating problem will have a unique solution (up to obvious sign changes) if

|e_{k}^{T} u_{k}| > |e_{k+1}^{T}v_{k}|

in which case

( \begin{matrix}u_{k+1}^{T} \\ v_{k+1}^{T} \end{matrix} ) = H (θ_{k}) ( \begin{matrix}u_{k}^{T} Z^{T} \\ v_{k}^{T} \end{matrix} )  (2.3)

where

sin θ_{k} = e_{k+1}^{T}v_{k}/e_{k}^{T} u_{k}  (2.4a)
cos θ_{k} = \sqrt{1 − sin^{2} θ_{k}}  (2.4b)

It is easy to verify that

e_{j}^{T} u_{k+1} = 0 , j < k + 1  (2.5a)
e_{j}^{T} v_{k+1} = 0 , j ≤ k + 1  (2.5b)

The calculation of u_{k+1}, v_{k+1} via (2.3), (2.4) can be performed in the obvious manner. Following common usage such algorithms will be referred to as hyperbolic downdating. However, some computational advantages may be obtained by rearranging (2.3).

For example, we may rewrite (2.3) as

v_{k+1} = (v_{k} − sin θ_{k}Zu_{k})/ cos θ_{k}  (2.6a)
u_{k+1} = − sin θ_{k}v_{k+1} + cos θ_{k}Zu_{k}  (2.6b)

which is obtained by using the second component of (2.3) to eliminate v_{k} in the first component of (2.6). It is of course possible to construct an alternative algorithm by using the first component of (2.3) to eliminate Zu_{k} in the second component of (2.3). It is easily seen that a single step of (2.3) or (2.6) requires roughly 4(n − k) + O(1) multiplications.

The relations (2.4a,b), (2.6a,b) form the basis for an alternative mixed algorithm for solving the elementary downdating problem. Although the computational cost of the two approaches are the same, it turns out that algorithms based on (2.4) (2.6) have superior stability properties to algorithms based on the obvious implementation of (2.3) (2.4) (see [[5] ]). However, while this superior stability is crucial in some applications such the downdating of a Cholesky decomposition, it is of lesser importance in the construction of efficient algorithms for the solution of symmetric positive definite Toeplitz system of equations.

It turns out that a substantial increase in efficiency can be achieved by considering the following modified downdating problem. Given α_{k}, β_{k} and w_{k}, x_{k} ∈ ℜ^{n} that satisfy

e_{j}^{T} w_{k} = 0 , j < k  (2.7a)
e_{j}^{T} x_{k} = 0 , j ≤ k  (2.7b)

find α_{k+1}, β_{k+1} and w_{k+1}, x_{k+1} ∈ ℜ^{n} that satisfy

α_{k+1}^{2}w_{k+1}w_{k+1}^{T} − β_{k+1}^{2}x_{k+1}x_{k+1}^{T} = α_{k}^{2}Zw_{k}w_{k}^{T} Z^{T} − β_{k}^{2}x_{k}x_{k}^{T}  (2.8a)

and

e_{j}^{T} w_{k+1} = 0 , j < k + 1  (2.8b)
e_{j}^{T} x_{k+1} = 0 , j ≤ k + 1  (2.8c)

If we make the identification

u_{k} = α_{k}w_{k}
v_{k} = β_{k}x_{k}

we find that the modified elementary downdating problem is equivalent to the elementary downdating problem. However, the extra parameters in the problem can be chosen judiciously to eliminate some multiplications. For example, if we take α_{k} = β_{k}, α_{k+1} = β_{k+1},

sin θ_{k} = e_{k+1}^{T}x_{k}/e_{k}^{T} w_{k}  (2.9a)
α_{k+1} = α_{k}/ cos θ_{k}  (2.9b)

and

w_{k+1} = Zw_{k} − sin θ_{k}x_{k}  (2.10a)
x_{k+1} = − sin θ_{k}Zw_{k} + x_{k}  (2.10b)

we obtain a basis for an algorithm which clearly requires about half the number of multiplications as an algorithm based on (2.3) (2.4). Similarly, we can obtain a mixed modified algorithms via

sin θ_{k} = β_{k}e_{k+1}^{T}x_{k}/α_{k}e_{k}^{T} w_{k}  (2.11a)
α_{k+1} = α_{k} cos θ_{k}  (2.11b)
β_{k+1} = β_{k}/ cos θ_{k}  (2.11c)

x_{k+1} = x_{k} − \frac{sin θ_{k}α_{k}}{β_{k}} Zw_{k}  (2.12a)
w_{k+1} = −\frac{sin θ_{k}β_{k+1}}{α_{k+1}} x_{k+1} + Zw_{k}  (2.12b)

The stability of modified mixed algorithms based on (2.11), (2.12) has been investigated in [13].

3 Cholesky Factorization.

Let us consider the following sequence of calculations. Given u, v ∈ ℜ^{n} with

v^{T} e_{1} = 0,

let

u_{1} = u, v_{1} = v

and for k = 1, …, n − 1 solve the elementary downdating problem defined by (2.1) (2.2), which we assume for the moment has a solution for each k. Now let

U = ∑_{k=1}^{n} e_{k}u_{k}^{T}  (3.1)

and

T = U^{T} U.  (3.2)

Equation (2.2a) is

u_{k+1}u_{k+1}^{T} − v_{k+1}v_{k+1}^{T} = Zu_{k}u_{k}^{T} Z^{T} − v_{k}v_{k}^{T}

and, on summing k from 1 to n − 1 we obtain

T − ZT Z^{T} = u_{1}u_{1}^{T} − v_{1}v_{1}^{T} − (Zu_{n}u_{n}^{T} Z^{T} − v_{n}v_{n}^{T})

If we now observe from (2.4.) that

Zu_{n} = v_{n} = 0

we obtain

T − ZT Z^{T} = uu^{T} − vv^{T}  (3.3)

In the terminology of Kailath, Kung and Morf [15], the matrix T has (+)− displacement rank 2. Furthermore, from (3.2), T is symmetric positive definite while from (2.1), (3.1) U is upper triangular. Thus (3.2) is the Cholesky decomposition of T.

We now return to the question of existence of a solution to the elementary downdating problem for each k = 1, …, n − 1. It is easy to verify that if T is positive definite and satisfies (3.3) with

e_{1}^{T} v = 0

then

|e_{k}^{T} u_{k}| > |e_{k+1}^{T}v_{k}|, k = 1, …, n − 1.

Consequently, the elementary downdating problem then has a solution for each k = 1, …, n − 1.

We have derived, albeit in a rather indirect manner, the basis of an algorithm for calculating the Cholesky decomposition of a positive definite matrix T that satisfies (3.3). In summary, for a positive definite matrix T satisfying

T − ZT Z^{T} = uu^{T} − vv^{T}
e_{1}^{T} v = 0

let

u_{1} = u, v_{1} = v

and for k = 1, …, n − 1 we calculate u_{k+1}, v_{k+1} such that for k = 1, …, n − 1

u_{k+1}u_{k+1}^{T} − v_{k+1}v_{k+1}^{T} = Zu_{k}u_{k}^{T} Z^{T} − v_{k}v_{k}^{T}  (3.4a)
e_{k+1}^{T}v_{k+1} = 0  (3.4b)

Then,

T = U^{T} U

U = ∑_{k=1}^{n} u_{k}u_{k}^{T} .

It is clear from Section 2 that the above algorithm requires only O(n^{2}) operations. Note that if T is a Toeplitz matrix, then (3.3) holds with

u^{T} = (t_{0}, t_{1}, …, t_{n−1})^{T} /\sqrt{t_{0}}
v^{T} = (0, t_{1}, …, t_{n−1})^{T} /\sqrt{t_{0}}

and thus the above yields an efficient means of calculation the Cholesky decomposition of a positive definite Toeplitz matrix. On closer examination it is not difficult to see that the above formulas are essentially equivalent to Bareiss algorithm [2] for calculating the Cholesky decomposition of a positive definite Toeplitz matrix. In fact we have a class of the Bareiss type Toeplitz solvers where each solver corresponds to a particular way of realizing the elementary downdating steps. For example, the connection with the modified elementary downdating problem is straightforward. On making the identification

u_{k} = α_{k}w_{k}
v_{k} = β_{k}x_{k}

we obtain

α_{k+1}^{2}w_{k+1}w_{k+1}^{T} − β_{k+1}^{2}x_{k+1}x_{k+1}^{T} = α_{k}^{2}Zw_{k}w_{k}^{T} Z^{T} − β_{k}^{2}x_{k}x_{k}^{T}
e_{k+1}^{T}x_{k+1} = 0

which is the modified elementary downdating problem. Then,

T = W^{T} D^{2}W

where

W = ∑_{k=1}^{n} e_{k}w_{k}^{T}

D = ∑_{k=1}^{n} α_{k}e_{k}e_{k}^{T} .

As noted previously, the modified elementary downdating problem has computational advantages and therefore leads to more efficient algorithms for the factorization of a positive definite Toeplitz matrix. However, in the sequel we will not dwell in the precise details of algorithms. Rather we shall relate algorithms based on the modified elementary downdating problem to those based on the elementary downdating problem. Thus for example if α̃_{k}, β̃_{k} w̃_{k} and x̃_{k} are the result of (2.9), (2.10) or of (2.11), (2.12) we define

ũ_{k} = α̃_{k}w̃_{k}
ṽ_{k} = β̃_{k}x̃_{k}

and then analyse the properties of

Ũ = ∑_{k=1}^{n} e_{k}ũ_{k}^{T} .

4 Analysis of Factorization Algorithms.

Here we present an analysis of the Cholesky factorization algorithm (3.4) for a positive definite matrix satisfying (3.3). Next the result of the analysis is applied to the case when the matrix is Toeplitz.

Let ũ_{k}, ṽ_{k} be the values of u_{k}, v_{k} that are essentially computed. We may assume that ũ_{k}, ṽ_{k} will satify a perturbed version of (3.4). That is,

ũ_{k+1}ũ_{k+1}^{T} − ṽ_{k+1}ṽ_{k+1}^{T} = Zũ_{k}ũ_{k}^{T} Z^{T} − ṽ_{k}ṽ_{k}^{T} + ϵG_{k} + O(ϵ^{2})  (4.1)

where ϵ is the machine relative precision and G_{k} depends on the precise specification of the algorithm used. It is not difficult to show that

ũ_{k} = u_{k} + O(ϵ), ṽ_{k} = v_{k} + O(ϵ)

and the aim of this section is to provide a first order analysis of the error. By a first order analysis we mean that we are bounding the first term in an asymptotic expansion for the error in powers of ϵ. Thus, one has to think of ϵ → 0 while the problem remains fixed, see [17].

On summing (4.1) from k between 1 and n − 1 we obtain

T̃ − ZT̃Z^{T} = ũ_{1}ũ_{1}^{T} − ṽ_{1}ṽ_{1}^{T} − (Zũ_{n}ũ_{n}^{T} Z^{T} − ṽ_{n}ṽ_{n}^{T}) + ϵ ∑_{k=1}^{n−1} G_{k} + O(ϵ^{2}).

where

T̃ = Ũ^{T} Ũ
Ũ = ∑_{k=1}^{n} e_{k}ũ_{k}^{T}

Since

Zũ_{n} = O(ϵ), ṽ_{n} = O(ϵ)

we find that

T̃ − ZT̃Z^{T} = ũ_{1}ũ_{1}^{T} − ṽ_{1}ṽ_{1}^{T} + ϵ ∑_{k=1}^{n−1} G_{k} + O(ϵ^{2}).  (4.2)

Now define

Ẽ = T̃ − T.  (4.3)

Then using (3.3) and (4.2)

Ẽ − ZẼZ^{T} = ũ_{1}ũ_{1}^{T} − uu^{T} + ṽ_{1}ṽ_{1}^{T} − vv^{T} + ϵ ∑_{k=1}^{n−1} G_{k} + O(ϵ^{2}).

from which it follows that

Ẽ = ∑_{j=0}^{n−1} Z_{j}[(ũ_{1}ũ_{1}^{T} − u_{1}u_{1}^{T}) + (ṽ_{1}ṽ_{1}^{T} − v_{1}v_{1}^{T})]Z_{j}^{T} + ϵ ∑_{j=0}^{n−1} ∑_{k=1}^{n−1} Z_{j}G_{k}Z_{j}^{T} + O(ϵ^{2}).  (4.4)

We see therefore that the error consists of two parts – the first associated with initial errors and the second associated with the fact that (4.1) contains an inhomogeneous term.

Now

∥ũ_{1}ũ_{1}^{T} − uu^{T} ∥ ≤ 2∥u∥ ∥ũ_{1} − u∥ + O(ϵ^{2})
∥ṽ_{1}ṽ_{1}^{T} − vv^{T} ∥ ≤ 2∥v∥ ∥ṽ_{1} − v∥ + O(ϵ^{2})

Furthermore, from (3.3)

T r(T) − T r(ZT Z^{T}) = ∥u∥^{2} − ∥v∥^{2} > 0

and hence

∥ ∑_{j=0}^{n−1} Z_{j}[ũ_{1}ũ_{1}^{T} − uu^{T} + ṽ_{1}ṽ_{1}^{T} − vv^{T}]Z_{j}^{T} ∥ ≤ 2n∥u∥[∥ũ_{1} − u∥ + ∥ṽ_{1} − v∥] + O(ϵ^{2})  (4.5)

This demonstrates that initial errors do not propagate unduly. To investigate the second term we require a preliminary result.

Lemma 4.1

∥Z_{j}v_{k}∥ ≤ ∥Z_{j+1}u_{k}∥

Proof Let

T_{k} = T − ∑_{l=1}^{k} u_{l}u_{l}^{T}

= ∑_{l=k+1}^{n} u_{l}u_{l}^{T}

Then, it is easy to verify that

T_{k} − ZT_{k}Z^{T} = Zu_{k}u_{k}^{T} Z^{T} − v_{k}v_{k}^{T}

and since T_{k} is positive semi-definite

T r[Z_{j}T_{k}Z_{j}^{T} − Z_{j+1}T_{k}Z_{j+1}^{T}] = ∥Z_{j+1}u_{k}∥^{2} − ∥Z_{j}v_{k}∥^{2} ≥ 0

□

We now demonstrate stability when the inhomogeneous term G_{k} satisfies

∥Z_{j}G_{k}Z_{j}^{T} ∥ ≤ c{∥Z_{j+1}u_{k}∥^{2} + ∥Z_{j}v_{k}∥^{2} + ∥Z_{j}u_{k+1}∥^{2} + ∥Z_{j}v_{k+1}∥^{2}}  (4.6)

and c is a positive constant. This bound is satisfied by a mixed downdating strategy such as (2.6) or (2.11) (2.12) but does not hold in general. For further details of the precise error analysis of mixed downdating algorithms, see [13], [5].

Theorem 4.1 Let ũ_{k}, ṽ_{k} satisfy (4.1.) where G_{k} satisfies (4.5). Then

∥T − Ũ^{T} Ũ∥ ≤ 2n∥u∥{∥ũ_{1} − u∥ + ∥ṽ_{1} − v∥} + 4ϵc ∑_{j=0}^{n−1} T r(Z_{j}T Z_{j}^{T}) + O(ϵ^{2})

Proof Using lemma 4.1,

∥Z_{j}G_{k}Z_{j}^{T} ∥ ≤ 2c{∥Z_{j+1}u_{k}∥^{2} + ∥Z_{j}u_{k+1}∥^{2}}.

Furthermore, since

T r(Z_{j}T Z_{j}) = ∑_{k=1}^{n} ∥Z_{j}u_{k}∥^{2},

it follows that

∥ ∑_{j=0}^{n−1} ∑_{k=1}^{n} Z_{j}G_{k}Z_{j}^{T} ∥ ≤ 4c ∑_{j=0}^{n−1} T r(Z_{j}T Z_{j}^{T})  (4.7)

The result now follows on combining (4.2) (4.3) and (4.5).

□

Remark A somewhat less general version of Theorem 4.1 has been obtained by Sweet in [18].

□

In general, an estimate of the form (4.5) does not hold. A weaker bound on G_{k} which does hold for all schemes outlined in section 2 is

∥Z_{j}G_{k}Z_{j}^{T} ∥ ≤ \frac{c}{4}∥H(θ_{k})∥(∥Z_{j+1}u_{k}∥ + ∥Z_{j}v_{k}∥)(∥Z_{j}u_{k+1}∥ + ∥Z_{j}v_{k+1}∥)  (4.8)

which, by Lemma 4.1. simplifies to

∥Z_{j}G_{k}Z_{j}^{T} ∥ ≤ c∥H(θ_{k})∥ ∥Z_{j+1}u_{k}∥ ∥Z_{j}u_{k+1}∥

The essential difference between (4.4) and (4.6) is the occurence of the multiplier ∥H(θ_{k})∥ which can be quite large. This term does cause difficulties in other applications such as the downdating of a Cholesky decomposition (see [5]) but, because of the special structure of the matrix T is of lesser importance here. The key result which allows us to obtain stability when (4.7) holds is:

Lemma 4.2

∥H(θ_{k})∥ ∥Z_{j}u_{k+1}∥ ≤ 2(n − k − j)∥Z_{j+1}u_{k}∥.

Proof It is easy to verify from (2.3) that

\frac{1 ∓ sin θ_{k}}{cos θ_{k}} \{u_{k+1} ∓ v_{k+1}\} = Zu_{k} ∓ v_{k}

and from (1.1) that

∥H(θ_{k})∥ = \frac{1 + | sin θ|}{cos θ}.

Thus,

∥H(θ_{k})∥ ∥Z_{j}u_{k+1}∥ ≤ ∥H(θ_{k})∥ ∥Z_{j}v_{k+1}∥ + ∥Z_{j+1}u_{k}∥ + ∥Z_{j}v_{k}∥
≤ ∥H(θ_{k})∥ ∥Z_{j+1}u_{k+1}∥ + 2∥Z_{j+1}u_{k}∥

where the last inequality was derived via lemma 4.1. Thus

∥H(θ_{k})∥ ∥Z_{j}u_{k+1}∥ ≤ 2 ∑_{l=j+1}^{n−k} ∥Z_{l}u_{k}∥

and the result follows.

□

Remark: Lemma 4.2 does not hold for the computed quantities unless we introduce an O(ϵ) term. However in a first order analysis we only need it to hold for the exact quantities.

Theorem 4.2 Let ũ_{k}, ṽ_{k} satisfy (4.1) where G_{k} satisfies (4.6). Then

∥T − Ũ^{T} Ũ∥ ≤ 2nϵ∥u∥{∥ũ_{1} − u∥ + ∥ṽ_{1} − v∥} + 2ϵc ∑_{j=1}^{n−1}(n − j)T r(Z_{j}T Z_{j}^{T})

Proof Applying Lemma 4.2 to (4.6) yields

∥Z_{j}G_{k}Z_{j}^{T} ∥ ≤ 2c(n − j − 1)∥Z_{j+1}u_{k}∥^{2}

and hence

∥ ∑_{j=0}^{n−1} ∑_{k=1}^{n−1} Z_{j}G_{k}Z_{j}^{T} ∥ ≤ 2c ∑_{j=1}^{n−1} ∑_{k=1}^{n−1}(n − j)∥Z_{j}u_{k}∥^{2}

≤ 2c ∑_{j=1}^{n−1}(n − j)T r(Z_{j}T Z_{j}^{T})  (4.9)

The result now follows from (4.2) (4.3) and (4.7).

□

Note that when T is Toeplitz,

T r(Z_{j}T Z_{j}^{T}) = (n − j)t_{0} .

Hence from Theorem 4.1 and 4.2 we obtain the main result on the stability of the Bareiss type algorithms for solving a symmetric positive definite Toeplitz system:

Corolarry 4.1: The Bareiss type Toeplitz solvers produce an upper triangular matrix Ũ such that

T = Ũ^{T} Ũ + ΔT

where ∥ΔT ∥ ≤ ϵt_{0}n^{2} when mixed downdating scheme is used, or ∥ΔT ∥ ≤ ϵt_{0}n^{3} when hyperbolic downdating scheme is used instead.

□

However, we have been unable to obtain numerical results that demonstrate that a ‘mixed downdating strategy’ yields superior results to ‘hyperbolic downdating’.

5 Preliminary Results on Basic Recurrences.

In this section we analyse properties of a nonlinear recurrence that is related to the Levinson recurrence for solving Yule-Walker problem. We then derive bounds on the growth of the magnitude of the terms in the recurrence. These bounds are applied in the next section in the roundoff error analysis of Levinson type algorithms.

Let us first consider the recurrence

( \begin{matrix}p_{k+1}^{T} \\ q_{k+1}^{T} \end{matrix} ) = H(θ_{k}) ( \begin{matrix}p_{k}^{T} Z^{T} \\ q_{k}^{T} \end{matrix} ) + ( \begin{matrix}f_{k}^{T} \\ g_{k}^{T} \end{matrix} )  (5.1)

where p_{1}, q_{1}, f_{1}, g_{1} are n-dimensional vectors. Then, it follows from (2.3) and (5.1) that

p_{k+1}u_{k+1}^{T} − q_{k+1}v_{k+1}^{T} = Zp_{k}u_{k}^{T} Z^{T} − q_{k}v_{k}^{T} + f_{k}u_{k+1}^{T} − g_{k}v_{k+1}^{T}

and hence that

∑_{k=1}^{n} p_{k}u_{k}^{T} − Z ∑_{k=1}^{n} p_{k}u_{k}^{T} Z^{T} = p_{1}u^{T} − q_{1}v^{T} + ∑_{k=1}^{n−1}(f_{k}u_{k+1}^{T} − g_{k}v_{k+1}^{T})

Thus,

P U = ∑_{j=0}^{n−1} Z_{j}(p_{1}u^{T} − q_{1}v^{T})Z_{j}^{T} + ∑_{j,k} Z_{j}(f_{k}u_{k+1}^{T} − g_{k}v_{k+1}^{T})Z_{j}^{T}

where

P = ∑_{k=1}^{n} p_{k}e_{k}^{T}

As a result,

p_{n} = {∑_{j=0}^{n−1} Z_{j}(p_{1}u^{T} − q_{1}v^{T})Z_{j}^{T} + ∑_{j,k} Z_{j}(f_{k}u_{k+1}^{T} − g_{k}v_{k+1}^{T})Z_{j}^{T} }r_{n}  (5.2)

where

r_{n} = U^{−1}e_{n}

is the last column of U^{−1}.

Our motivation for examining (5.1) is its close connection with the recurrence

y_{k+1} = Zy_{k} − sin θ_{k}J_{k}y_{k} + \sqrt{c_{k+1}}∥H(θ_{k})∥h_{k}  (5.3)

where

J_{k} = JZ_{n−k}
c_{k} = t_{0} ∏_{j=1}^{k−1} cos^{2} θ_{j}
e_{j}^{T} y_{1} = 0, j > 1
e_{j}^{T} h_{k} = 0, j > k + 1.

It is easy to verify that the solution of (5.3) satisfies

e_{j}^{T} y_{k} = 0, j > k

and hence

J_{k+1}J_{k}y_{k} = Zy_{k}
J_{k+1}Zy_{k} = J_{k}y_{k}

If we now define

p_{k} = c_{k}^{−1/2}y_{k}
q_{k} = c_{k}^{−1/2}J_{k}y_{k}

we find that p_{k}, q_{k} satisfy (5.1) with

f_{k} = ∥H(θ_{k})∥h_{k}
g_{k} = ∥H(θ_{k})∥J_{k+1}h_{k}

Hence, it follows from (5.2) that

y_{n} = {t_{0}^{−1/2} ∑_{j=0}^{n−1} Z_{j}y_{1}(u − v)^{T} Z_{j}^{T}
+ ∑_{j,k} ∥H(θ_{k})∥Z_{j}(h_{k}u_{k+1}^{T} − J_{k+1}h_{k}v_{k+1}^{T})Z_{j}^{T} \}\sqrt{c_{n}}r_{n}  (5.4)

Now let a_{k} satisfy the homogeneous equation

a_{k+1} = Za_{k} − sin θ_{k}J_{k}a_{k}  (5.5a)
a_{1} = e_{1}  (5.5b)

Then, from (5.4)

a_{n} = t_{0}^{−1/2} ∑_{j=0}^{n−1} Z_{j}e_{1}(u − v)^{T} Z_{j}^{T} \sqrt{c_{n}}r_{n}

We now consider only the case when T is Toeplitz. Then

u − v = t_{0}^{1/2}e_{1}  (5.6)

and hence

a_{n} = \sqrt{c_{n}}r_{n}  (5.7)

Hence from (5.4) and (5.7) the solution of (5.3) is

y_{n} = {∑_{j=0}^{n−1} Z_{j}y_{1}e_{j+1}^{T} + ∑_{j,k} ∥H(θ_{k})∥Z_{j}(h_{k}u_{k+1}^{T} − J_{k+1}h_{k}v_{k+1}^{T})Z_{j}^{T} }a_{n}  (5.8)

Remark It is easy to obtain solutions of (5.3) when the restrictions on y_{1} and h_{k} are relaxed by using (5.8) and superposition.

Lemma 5.1 If (5.6) holds then the solution of (5.3) satisfies

∥y_{n}∥ ≤ {|e_{1}^{T} y_{1}| + 4n^{2}t_{0}^{1/2}(∑_{k=1}^{n−1} ∥h_{k}∥_{1}^{2})^{1/2}}∥a_{n}∥

Proof It is straightforward to see that

∥ ∑_{j=0}^{n−1} Z_{j}y_{1}e_{j+1}^{T}∥ = |e_{1}^{T} y_{1}|  (5.9)

In addition, it is easy to verify that

∥ ∑_{j=0}^{n−k−1} Z_{j}h_{k}e_{j+1}^{T}∥_{1} = ∥ ∑_{j=0}^{n−k−1} Z_{j}h_{k}e_{j+1}^{T}∥_{∞} = ∥h_{k}∥_{1}.  (5.10)

and we now use (5.10) to estimate the second term in (5.8). We have

∥ ∑_{j=0}^{n−1} Z_{j}h_{k}u_{k+1}^{T}Z_{j}^{T} ∥ = ∥ ∑_{j=0}^{n−k−1} Z_{j}h_{k}u_{k+1}^{T}Z_{j}^{T} ∥

≤ ∥ ∑_{j=0}^{n−k−1} Z_{j}h_{k}e_{j+1}^{T}∥ ∥ ∑_{j=0}^{n−k−1} Z_{j}u_{k+1}e_{j+1}^{T}∥

≤ n^{1/2}∥h_{k}∥_{1}∥u_{k+1}∥

Similarly, using Lemma 4.1

∥ ∑_{j=0}^{n−1} Z_{j}J_{k+1}h_{k}v_{k+1}^{T}Z_{j}^{T} ∥ ≤ n^{1/2}∥h_{k}∥_{1}∥u_{k+1}∥

Thus,

∥ ∑_{j,k} ∥H(θ_{k})∥Z_{j}(h_{k}u_{k+1}^{T} − J_{k+1}h_{k}v_{k+1}^{T})Z_{j}^{T} ∥

≤ 2n^{1/2} ∑_{k=1}^{n−1} ∥H(θ_{k})∥ ∥h_{k}∥_{1}∥u_{k+1}∥

≤ 2n^{1/2}(∑_{k=1}^{n−1} ∥H(θ_{k})∥^{2}∥u_{k+1}∥^{2})^{1/2}(∑_{k=1}^{n−1} ∥h_{k}∥_{1}^{2})^{1/2}

≤ 4n^{3/2}(∑_{k=1}^{n−1} ∥Zu_{k}∥^{2})^{1/2}(∑_{k=1}^{n−1} ∥h_{k}∥_{1}^{2})^{1/2}  (5.11)

where we have used Lemma 4.2 to establish the last inequality. In addition,

∑_{k=1}^{n−1} ∥Zu_{k}∥^{2} = T r(ZT Z^{T}) = (n − 1)t_{0}  (5.12)

and the result now follows from (5.8), (5.9), (5.11) and (5.12).

□

6 The Levinson–Durbin Algorithm.

In this section we consider the problem of calculating the last column of the inverse of a Toeplitz matrix. This is a slight variation of the usual Levinson–Durbin algorithm which is often associated with the solution of the Yule–Walker equations. However a simple modification of the analysis presented here enables us to analyse this and other modifications of the basic algorithm.

We begin by recalling (5.5) (5.6). That is, the solution of

a_{k+1} = (Z − sin θ_{k}J_{k})a_{k}, a_{1} = e_{1}  (6.1)

is

a_{n} = \sqrt{c_{n}}r_{n} = \sqrt{c_{n}}U^{−1}e_{n}

where

c_{k} = t_{0} ∏_{j=1}^{k−1} cos^{2} θ_{j}

and

T = U^{T} U

Thus,

T a_{n} = \sqrt{c_{n}}e_{n}^{T} u_{n}e_{n}

Furthermore, it is easy to verify from (3.4) that

e_{n}^{T} u_{n} = \sqrt{c_{n}}.

Hence,

T a_{n} = c_{n}e_{n}  (6.2)

and consequently c_{n}^{−1}a_{n} is the last column of T^{−1}. Clearly, the recurrence (6.1) provides an effective means of calculating a_{n} when the coefficients sin θ_{k} are known. Incidently, − sin θ_{k} are often referred to as the reflection or partial correlation coefficients. In Sections 2 and 3, we considered the calculation of sin θ_{k} via the auxiliary vectors u_{k}, v_{k} and this is essentially how they are determined in Bareiss’s algorithm. It is however possible to calculate them independently as follows

sin θ_{k} = v̂^{T} Za_{k}/c_{k}
c_{k+1} = (1 − sin^{2} θ_{k})c_{k}

where

v̂^{T} = (0, t_{1}, …, t_{n−1})^{T}
= t_{0}^{1/2}v^{T}

Thus, we obtain

a_{1} = e_{1} , c_{1} = t_{0}  (6.3a)
sin θ_{k} = v̂^{T} Za_{k}/c_{k}  (6.3b)
c_{k+1} = (1 − sin^{2} θ_{k})c_{k}  (6.3c)
a_{k+1} = (Z − sin θ_{k}J_{k})a_{k}  (6.3d)

and the last column of T^{−1} is given by c_{n}^{−1}a_{n}. Suppose that the computed values of a_{k} satisfy

ã_{k+1} = (Z − s̃_{k}J_{k})ã_{k} + ϵ\sqrt{c_{k+1}}∥H(θ_{k})∥ξ_{k}  (6.4)

where s̃_{k} = sin θ̃_{k} is the computed value of sin θ_{k}.

Then on multiplying (6.4) by E_{k+1}T we obtain

E_{k+1}T ã_{k+1} = (Z − s̃_{k}J_{k})E_{k}T ã_{k}
+ v̂^{T} Zã_{k}e_{1} − s̃_{k}v̂^{T} Zã_{k}e_{k+1}
+ ϵ\sqrt{c_{k+1}}∥H(θ_{k})∥E_{k+1}T ξ_{k}

If we now write

E_{k}T ã_{k} = c̃_{k}e_{k} + δ_{k}  (6.5)

where c̃_{k} is the calculated value of c_{k} and note that

δ_{k} = O(ϵ),

we obtain, to first order

δ_{k+1} = (Z − s_{k}J_{k})δ_{k} + ϵ\sqrt{c_{k+1}}∥H(θ_{k})∥\{α_{k}e_{1} + β_{k}e_{k+1} + E_{k+1}T ξ_{k}\}

ϵ\sqrt{c_{k+1}}∥H(θ_{k})∥α_{k} = v̂^{T} Zã_{k} − s̃_{k}c̃_{k}
ϵ\sqrt{c_{k+1}}∥H(θ_{k})∥β_{k} = c̃_{k} − c̃_{k+1} − s̃_{k}v̂^{T} Zã_{k}

Rounding errors that could be obtained from (6.3) and (6.4) are of the form

∥ξ_{k}∥ ≤ K∥r_{k}∥, r_{k} = U^{−1}e_{k}
|α_{k}| ≤ Kk∥r_{k}∥
|β_{k}| ≤ Kk∥r_{k}∥.

and error estimates for δ_{n} can now be estimated via Lemma 5.1.

Let A = (ã_{1}, ã_{2}, ⋯, ã_{n}). It follows from (6.1)-(6.4) that A is upper triangular such that

T A = L + ΔL  (6.6)

where L is lower triangular, ΔL = (δ_{1}, δ_{2} ⋯, δ_{n}) is upper triangular and the columns of ΔL are bounded as indicated by Lemma 5.1, that is,

∥δ_{k}∥ ≤ ϵ(∑_{i=1}^{k−1} ∥r_{i}∥^{2})^{1/2}∥a_{k}∥  (6.7)

Note now that

A^{T} T A = A^{T} L + A^{T} ΔL  (6.8)

The matrix A^{T} T A is symmetric so A^{T} L + A^{T} ΔL also has to be symmetric. But A^{T} L + A^{T} ΔL is ”almost” lower triangular thus

A^{T} L + A^{T} ΔL = D + ΔD , ∥ΔD∥ ≤ ∥A∥∥ΔL∥ ≤ ϵ∥U^{−1}∥∥A∥^{2}  (6.9)

and also

T = A^{−T} DA^{−1} + A^{−T} ΔDA^{−1}  (6.10)

Now the computed solution x̃ is obtained as x̃ = AD^{−1}A^{T} b. Hence

x − x̃ = (I − AD^{−1}A^{T} T)x = −AD^{−1}ΔDA^{−1}x  (6.11)

or

\frac{∥x − x̃∥}{∥x∥} ≤ ∥AD^{−1}ΔDA^{−1}∥  (6.12)

References

[1] S.T. Alexander, C-T. Pan and R.J. Plemmons, “Analysis of a Recursive least Squares Hyperbolic Rotation Algorithm for Signal Processing”, Linear Algebra and Its Applications, vol 98, pp 3-40, 1988.

[2] J.R. Bareiss, “Numerical Solution of Linear Equations with Toeplitz and Vector Toeplitz Matrices”, Numerische Mathematik, vol 13, pp 404-424, 1969.

[3] A.W. Bojanczyk, R.P. Brent and F.R. De Hoog, “QR Factorization of Toeplitz Matrices”, Numerische Mathematik, vol 49, pp 81-94, 1986.

[4] A.W. Bojanczyk, R.P. Brent and F.R. De Hoog, “Stability Analysis of Fast Toeplitz Linear System Solvers” Technical Report CMA-MR17-ANU-91, Center for Mathematical Analysis, The Australian National University, August 1991.

[5] A.W. Bojanczyk, R.P. Brent, P. Van Dooren and F.R. De Hoog, “A Note on Downdating the Cholesky Factorization”, SIAM JSSC , vol 8, pp 210-220, 1987.

[6] A.W. Bojanczyk and A. Steinhardt, “Matrix Downdating Techniques for Signal Processing”, Proceedings of the SPIE Conference on Advanced Algorithms and Architectures for Signal Processing III, vol 975, pp 68-75, 1988.

[7] J.R. Bunch, “The Weak and Strong Stability of Algorithms in Numerical Linear Algebra”, Linear Algebra and Its Applications, vol 88/89, pp 49-66, 1987.

[8] J. Chun, T. Kailath and H. Lev-Ari, “Fast Parallel Algorithms for QR and Triangular Factorization”, SIAM JSSC , vol 8, pp 899-913, November 1987.

[9] G. Cybenko, “The Numerical Stability of the Levinson-Durbin Algorithm for Toeplitz Systems of Equations”, SIAM JSSC , vol 1, no 3, pp 303-319, September 1980.

[10] G. Cybenko, “The Numerical Stability of the Lattice Algorithm for Least Squares Linear Prediction Problems”, BIT ,

[11] I. Gohberg, I. Koltracht and D. Xiao, “On the solution of the Yule-Walker equation”, Proceedings of the SPIE Conference on Advanced Algorithms and Architectures for Signal Processing IV, vol 1566, July 1991.

[12] G.H. Golub and C. Van Loan, Matrix Computations, John Hopkins Press, Baltimore, MD, 1983.

[13] R. Fletcher and M.J.D. Powell, “On the Modification of LDL^{T} Factorizations”, Mathematics of Computation, vol 28, no 128, pp 1067-87, October 1974.

[14] M. Jankowski and H. Wozniakowski, “Iterative Refinement Implies Numerical Stability”, BIT , vol 17, pp 303-311, 1977.

[15] T. Kailath, S.Y. Kung and M. Morf, “Displacement Ranks of Matrices and Linear Equations”, J. Math. Anal. Appl., vol 68, pp 395-407, 1979.

[16] I. Koltracht and P. Lancaster, “Condition Numbers of Toeplitz and Block Toeplitz Matrices”, in Operator Theory: Advances and Applications vol 18, Birkhauser Verlag Basel, 1986.

[17] W. Miller and C. Wrathall, Software for Roundoff Analysis of Matrix Algorithms, Academic Press, 1980.

[18] D. Sweet, “Numerical Methods for Toeplitz Matrices”, PhD thesis, University of Adelaide, 1982.

[19] J.M.Varah, “The Prolate Matrix: A Good Toeplitz Test Example“, SIAM Conference on Control, Signals and Systems, November 5-8, San Francisco, 1990.