Stability Issues in the Factorization of Structured Matrices.
This paper provides an error analysis of the generalized Schur algorithm of Kailath and Chun [SIAM J. Matrix Anal. Appl., 15 (1994), pp. 114--128]---a class of algorithms which can be used to factorize Toeplitz-like matrices, including block-Toeplitz matrices, and matrices of the form $T^{T}T$, where $T$ is Toeplitz. The conclusion drawn is that if this algorithm is implemented with hyperbolic transformations in the factored form which is well known to provide numerical stability in the context of Cholesky downdating, then the generalized Schur algorithm will be stable.  If a more direct implementation of the hyperbolic transformations is used, then it will be unstable.  In this respect, the algorithm is analogous to Cholesky downdating; the details of implementation of the hyperbolic transformations are essential for stability.  An example which illustrates this instability is given.  This result is in contrast to the ordinary Schur algorithm for which an analysis by Bojanczyk, Brent, De Hoog, and Sweet [SIAM J. Matrix Anal. Appl., 16 (1995), pp. 40--57] shows that the sta- bility of the algorithm is not dependent on the implementation of the hyperbolic transformations.
Introduction
. The Schur algorithm is a popular and fast method for the
Cholesky factorization of a square, positive definite Toeplitz matrix, T . It performs
reliably and in [3] it was shown to be stable in the sense that if the algorithm runs
to completion and -
C is the computed Cholesky factor,
Ck is guaranteed
to be small. This paper will perform a similar stability analysis which applies to
several special cases of the generalized Schur algorithm, [10]. In its full generality,
the generalized Schur algorithm can be adapted to the factorization of a wide variety
of structured matrices. The analysis given here is primarily of interest for block-
Toeplitz and Toeplitz-block matrices, as well as for matrices of the form T T T , where
T is rectangular and Toeplitz. The key notion behind the general algorithm is the
concept of displacement rank, [10].
One of the most significant examples is the Cholesky factorization of T T T . This
factor is also the factor R in the QR factorization of the rectangular Toeplitz matrix
T and the obvious application of this fact to the solution of Toeplitz least squares
problems is explored in [1]. However, the analysis given there assumes the use of the
algorithm presented in [2] rather than the generalized Schur algorithm. The basic
idea is to obtain R without bothering about finding Q, thus avoiding any problems
associated with the loss of orthogonality which are common to all fast Toeplitz QR
algorithms. The method of semi-normal equations, possibly with iterative refinement,
can then be used to find the least squares solution. The resulting equations are
and a weak stability result can be given concerning their solution
provided that for the computed R, kR As will be shown
here, this can be done using the generalized Schur algorithm as well as the algorithm
in [2]. Consequently, the observations concerning Toeplitz least squares problems in
[1] can be adapted to an alternate approach-that of finding R from the generalized
Schur algorithm and then solving the semi-normal equations to obtain the solution.
A comparison of the analysis in this paper can also be made with the analysis of
the Schur algorithm given in [3]. Despite some major similarities, the conclusions of
This work was supported by ARPA (grant 60NANB2D1272) and NSF (grant CCR-9209349)
y (stewart@sccm.stanford.edu).
z Dept. Mathematical Engineering, Universit'e Catholique de Louvain, Louvain-la-Neuve, Bel-
gium, (Vandooren@anma.ucl.ac.be).
this paper will be, in one very important respect, quite different: the implementation
details of the Schur algorithm are less critical to stability than those of the generalized
Schur algorithm. Both can be implemented using hyperbolic transformations of the
with where the transformations are applied with straightforward matrix
multiplication; however there are reasons to be concerned about the stability of this
approach. In many algorithms, such as Cholesky downdating, [4], this is not a good
idea. A downdating algorithm based on such transformations will not be stable unless
they are implemented in factored form,
ae
It is worth noting that the manner in which the elements of these factors are computed
is also important for stability. Details can be found in [4] and a correct implementation
can be found in Section 3. Such a seemingly minor difference in implementation makes
the difference between stability and instability for Cholesky downdating, but the
Schur algorithm is more robust-it is stable using either approach. Unfortunately,
the generalized Schur algorithm will be shown to be analogous to the downdating
problem in this regard; the factored transformations are essential for stability. A
proof of the stability of the factored hyperbolic approach will be given in Section 4
while an example which illustrates the instability of the other approach will be given
in Section 5. The stability proof in [1] also requires use of the factored hyperbolic
transformations. This is not surprising since the presentation in [2] recasts the problem
in the form of a downdating problem.
Finally, in Section (5), some comments will be made about the nature of the
instability. The propagation of errors is essentially stable regardless of which implementation
is used; the difference is in the size of the local errors. From the fact that
the error propagation is stable, it is possible to show that there will be many instances
in which the instability does not completely destroy the accuracy of the computed
factors. This is true even for quite ill-conditioned problems. Numerical
examples will be given to support this claim.
2. Displacement Rank and Stable Computation of the Initial Genera-
tors. Given an n \Theta n symmetric, positive-definite matrix A, which is not necessarily
Toeplitz, we define the displacement of A to be
(1)
where ZA is strictly lower triangular(a matrix with zeros on the diagonal). The
restriction of symmetry on A can easily be relaxed to allow A to be Hermitian, but
for simplicity we deal only with the real symmetric case here.
In [10], the only additional restriction on Z is that it be a fast, O(n) rather than
process to multiply a vector by Z. Otherwise, the factorization algorithm will
have a complexity of O(n 3 ) instead of O(n 2 ). In addition to this, we will also make
two assumptions which will guarantee the stability of the algorithm and simplify the
analysis. For stability it will be necessary to assume that kZk 2 - 1. Otherwise,
repeated multiplication by Z will magnify errors. Further, to simplify the analysis,
it will be assumed that a vector can be multiplied by Z without incurring any error.
These two assumptions essentially limit Z to be a shift or a block shift. Since these
are the forms that Z takes for Toeplitz least squares problems and block-Toeplitz
factorization respectively, these assumptions are reasonable, even if they
do remove a good deal of generality. The second assumption is not essential if the
multiplication is done in a stable manner, but there don't seem to be any common
examples to which the analysis could be made to apply by dropping it.
To be more specific about typical examples of Z, if A is square and Toeplitz or
has the form T T T for some rectangular Toeplitz matrix T , we choose (Z
would require ZB
to be a block shift, with the ones on the subdiagonal replaced with identity matrices.
The significance of the analysis here will be for algorithms based on displacements
involving Z T and ZB .
It is well known that if T is a symmetric Toeplitz matrix, normalized to have ones
on the diagonal, will be an indefinite
rank two matrix and
where
and \Sigma
The vectors u and v are referred to as generators of T . Clearly, the matrix T
is uniquely determined by these generators. Performing operations on the genera-
tors, rather than on T directly is what allows efficient, O(n 2 ) algorithms for Toeplitz
systems.
Assume that for some ZA satisfying the previously specified restrictions, A belongs
to a class of positive-definite matrices for which DA as defined by Equation (1) has
rank ff. Clearly DA will have a decomposition of the form
where
and \Sigma A = I p \Phi I q . Here we take u 11 and v 11 to be scalars and u 21 and v 21 to be
column vectors of size respectively, with
In Section 3, we will give a summary of the general algorithm. More details
and many special cases can be found in [10]. But first, since the numerical stability
of the overall factorization will depend on the stability of the process of finding the
generators for A, we will briefly discuss how this can be done in a stable manner.
As already seen, the generators of a Toeplitz matrix can be found trivially with
no error. For a symmetric block-Toeplitz matrix, the process is only slightly more
complicated and the errors in the initial generators will not cause a problem with
stability. The generators for a Toeplitz least squares problem can also be found in
a reliable manner with a minimal amount of computation, and again the process is
stable.
In general, the numbers p and q correspond respectively to the number of positive
and negative eigenvalues of D. When this decomposition cannot be obtained trivially,
as in the Toeplitz or block-Toeplitz case, it can be obtained through an eigenvalue
decomposition or through Gaussian elimination with a symmetric pivoting strategy
to obtain an LDL T factorization of the displacement, [7].
Of course, computing a full eigenvalue decomposition will slow down the overall
algorithm. In the absence of any specific knowledge of the form which the generators
will take, and in the case in which D can be computed with O(n 2 ) complexity,
computing ff steps of the LDL T decomposition will give a set of generators without
destroying the O(n 2 ) complexity of the algorithm. However it is important to note
that when truncating such a decomposition, the pivoting scheme which is chosen can
be more critical to stability than when the decomposition is carried out completely.
The fact that the Bunch-Parlett strategy is a backward stable method for computing
LDL T is well known, [7]. The stability of the Bunch-Kaufman procedure has been
established more recently in [9]. However, it is simple to verify numerically that the
17)=8 is a constant chosen to minimize element growth for symmetric
pivoting and ffl ? 0 is very small, leads to a large error if two steps of the Bunch-
Kaufman approach are used to obtain a rank 2 LDL T factorization. The Bunch-
Parlett procedure will not have this problem. The reason for the difference can be
found by looking at the sensitivity of the scalar Schur complement which is truncated
by these two algorithms. An analysis of pivoted low rank Cholesky factorization which
highlights the relevance of Schur complement sensitivity in the case of semi-definite
matrices was given in [8]; the issues are the same for indefinite matrices, and the
extension to the more complicated pivoting strategies is direct. The result can be
summarized as follows: backward stability guarantees that if the computed Schur
complement after two stages of the factorization process is small, then no large errors
will be incurred in dropping it and terminating the factorization. Rank deficiency
ensures that the exact value of the Schur complement will be zero, but if it is highly
sensitive to perturbations in the original matrix, then the backward error can cause
the computed value to become large. It is shown in [8] that the sensitivity of the
Schur complement will depend on the size of the elements in L. A careful scrutiny of
the two pivoting strategies shows that, although both are backward stable and give
bounds on the norms of the Schur complements, only Bunch-Parlett gives bounds
on the elements of L. However, this digression is included only for completeness,
to show that the process of obtaining the generators can always be done in a stable
fashion. In practice, there is generally a simpler means of obtaining the generators-
for Toeplitz matrices, they can be obtained with no computation at all. A general
factorization method would only be needed in the rare case in which D does not have
a zero structure which makes the choice of generators obvious, or at least easy to
compute. An example is the case of general non-Toeplitz matrices for which ff equals
two. The class of matrices for which broader than the class of Toeplitz
matrices, and DA will not always have a zero structure which makes the process of
finding the generators trivial.
3. The Generalized Schur Algorithm. Assuming that we have managed to
find the generators for a structured matrix A, we can apply the generalized Schur
algorithm to find a factorization Much of the power and generality of
this approach is due to the fact that it is a straightforward and intuitively clear
generalization of the Schur algorithm; there is little essential difference between the
cases 2. Since the former case is well known, the presentation here of
the latter case can be kept brief. A more leisurely description can be found in [10].
A . Since Z is strictly lower triangular, 0, and from
this it follows that
or
(2)
For any transformation J such that J T we find that JG is also a set of
generators for A. This follows immediately, since G T J T \SigmaJ
From Equation (2) and the positivity of A, we see that
means that if a transformation, J 1 , of the form
are orthogonal, is applied to G to give
U 22
then j-u Therefore it is possible to take ae = \Gamma-v 11 =-u 11 so that
~
U 22
Note that J T
so that J 2 J 1 G is also a set of generators for
A. We have proven that we can always find a set of generators for which the first
column has only one nonzero element. Such generators are said to be proper.
This form of the generators is very useful. Equation (2) implies that the first row
of A is not different from the first row of G T \SigmaG. When the generators are in proper
form, this means that
\Theta
~ u 11 ~
Clearly, the first row of the Cholesky factor of A will be
\Theta
~
If we let
~
- \Theta
~
be the Schur complement of A, we see that
~
~
- \Theta
~
~
~
- \Theta
~
This means that merely post-multiplying the first row of the generators,
\Theta ~
by Z T will give the generators of the Schur complement.
Additional unitary and hyperbolic transformations can now be used to introduce
zeros into the second elements of these vectors to obtain the first row of the Cholesky
factor of the Schur complement-the second row of the Cholesky factor of A. The
process can be continued recursively to obtain the complete Cholesky factor of A. At
each stage, the positivity of the Schur complements guarantee that an appropriate ae
can be found.
The resulting algorithm, with the hyperbolic transformations implemented in the
stable form in MATLAB c
code is as follows,
function [C]=gschur(G,n,p,q)
for i=1:n,
for
for
G(p+1,i+1:n)=[-s c]*G(1:p:p+1,i+1:n);
It was mentioned in Section 1 that the computation of the elements in the factors
is significant for stability. It is also important that the leading element of the first
generator be computed separately as shown here. The analysis in [4] does not apply
to other implementations.
Although the algorithm given above will suffice to find the factors of a block-
Toeplitz matrix, it is significantly different from some of the approaches given in the
literature. The numerical properties are not different, but a more block oriented
perspective on block-Toeplitz factorization can be found in [6].
In the 2-generator case, the matrix J 2 which zeros the first elements of all the
generators except the first is almost uniquely defined by the constraints
x-
\Sigma. The only possible variation is that each row of H can be multiplied
by \Gamma1. This is not the case when there are more generators. The obvious example
is the possibility of applying an arbitrary hyperbolic transformation acting on the
last generators after the appropriate zeros have been introduced into their
leading elements. This paper will show that when J is formed from a block diagonal
orthogonal matrix and an appropriately implemented hyperbolic transformation, the
algorithm will be stable.
4. Stability Analysis. A stability analysis of the Schur algorithm for Toeplitz
factorization is presented in [3]. The analysis can be broken into two parts: one which
shows how local errors propagate through the algorithm and one which bounds the
local errors. The propagation of errors is essentially the same for the generalized
Schur algorithm, but the problem of bounding local errors is slightly more difficult.
We will first modify results from [3] to apply to the general algorithm and then later
develop new inequalities which will complete the analysis.
At the end of the kth stage of the algorithm, let the generators stored in the
computer be
~
~
U 22;k
~
U 2;k
~
Using MATLAB c
G k;Z be defined by
~
G k;Z (1;
and
~
Also, let G k and G k;Z be the generators which would result from exact computations.
The computed generators will satisfy
~
where ffl is the machine precision and fflF k are errors incurred locally in computing the
new form of the generators at step k through the use of orthogonal and hyperbolic
transformations.
Since it is possible that G 1 will already be in proper form, we will assume that it
is and treat any errors which appear in it separately from the other errors. Summing
this equation from grouping the terms related to the top rows
of ~
G k and ~
G k;Z in the left hand side gives
~
G
G
After simplifying the terms that cancel in the right hand side, we obtain
~
R T ~
R T ~
where
~
~
~
~
~
R T ~
Using this to find the displacement
R T ~
R T ~
and applying Equation (2) gives
R T ~
where the sum over k has been reduced by noting that Z T
This shows that if we can bound errors in computation of the initial
generators and the local errors, then the algorithm will be stable. The errors in the
initial generators are not a problem, since from D the generators can be obtained in
a backward stable manner using Bunch-Parlett pivoting to compute LDL T or, more
typically, in a more direct fashion. Since the methods for obtaining the generators
may vary, in the analysis and error bounds which follow we will ignore this source of
error and only concern ourselves with the effects of local errors due to the unitary
and hyperbolic transformations.
The local errors are given by the expression
Because any bounds on the errors produced by the transformations will depend on
the norm of the generators, it is essential to bound the generators.
Theorem 4.1. When the generators are computed by applying a sequence of
plane rotations and a hyperbolic transformation, they satisfy
F
1 be the two generators u 1 and v 1 after orthogonal transformations
have been performed on the generators from step k \Gamma 1. Then
ae
ae
- u
Taking Frobenius norms gives
To see why the final inequality is true, let look for the minimum value of
on 0 - x - 1. If
then
and
It is easy to see that f assumes lower values at point is
the only possible maximum. Since the plane rotations don't affect the norm of the
generators and since kZk 2 - 1,
This inequality can be expanded recursively to give
Finally, for an arbitrary positive semi-definite, rank with a factorization
This follows from the fact that the Frobenius norm squared equals the sum of squares
of the singular values and from a standard inequality relating the vector 2-norm and
the vector 1-norm. This completes the proof of the theorem.
To complete a stability analysis all that is needed is to show that the orthogonal
and hyperbolic transformations produce a local error, fflF k , which is proportional to
the norm of the generators. Note that this does not necessarily refer to norm of
the generators prior to the transformation. In fact, in the case of the hyperbolic
transformation, it is necessary to look at the norm of one of the generators which is
produced by the hyperbolic transformation.
An error analysis of hyperbolic transformations is given in [4]. The result is that
if the transformations are applied in factored formp
ae
then
~
\Deltau T
\Deltav T
The mixed error vectors c
\Deltau and c
\Deltav satisfy
\Deltau T
c
\Deltav T
F
F
where ffl is the unit roundoff.
In addition to this, there is a result in [11] concerning the application of plane
rotations. In particular, there will exist orthogonal -
such that
~
~
U 2;k
where
~
U 2;k
F
and
~
F
If we let
\Theta
then clearly
Also, if we let
d
\Deltau c
\Deltav
then the error bounds, (5) and (7), can be used to show that
and
c
\Deltau T
c
\Deltau T
c
\Deltav T
c
\Deltav T
are standard basis vectors. This gives
\Theta
~
d
\DeltaG k
corresponding to a bound
Since Theorem 4.1 shows that
we get a bound on the local error of
ip
F
From (3), we see that
\Gammap
F
5. An Unstable Implementation. In the last section, the stability of the generalized
Schur algorithm was proven for the case in which the hyperbolic transformations
are applied in factored form. It has already been noted that the implementation
is not unique. There are many transformations which introduce the needed zeros in
the generators and they can be applied in multiple ways: the obvious example is the
application of the hyperbolic transformation in factored form, (4), versus the direct
multiplication approach.
It turns out that the factored form of the hyperbolic transformation is crucial to
the stability of the algorithm. To see this, take and the generators
For small j these generators correspond to an ill-conditioned A for which 4. If the
Schur algorithm with the hyperbolic transformations applied in factored form is used,
then a small error A will be achieved, independent of the size of j. However,
if the hyperbolic transformations are multiplied directly, then the backward error will
depend to a significant extent on j.
To see why, note that this example gives, after two steps of the generalized Schur
algorithm, the following generators
In fact, the example was obtained by performing a reversal of the algorithm on Equation
(12). The construction of the matrix in Equation (11) was not really necessary
since the essential mechanism of instability is represented in Equation (12); this
construction merely removes any objection as to whether the generators (12), corresponding
to a displacement rank 2 matrix, could actually occur in the factorization
of a displacement rank 4 matrix. If not, then one might argue that the instability
shown when applying the Schur algorithm to Equation (12) is not really relevant to
the stability of the factorization of displacement rank 4 matrices. Since the matrix
in Equation (11) corresponds to four linearly independent generators which reduce to
Equation (11) in two steps of the Schur procedure, this objection can be dismissed.

Table

5 shows the dependence of the errors on j. It was necessary to go to very
ill-conditioned matrices to demonstrate a significant loss of accuracy. The increase
of the backward error by a factor 10 5 is more than seems reasonable in a backward
stable algorithm, yet the increase in error is quite modest for an unstable algorithm
in which the condition number is increased by a factor of It turns out that this
is primarily a result of the fact that instability is due to large local errors rather than
unstable propagation of errors.
In contrast, over this very wide of condition numbers, the error in the Cholesky
factors computed by the stable form of the algorithm are all around 10 \Gamma15 . This is
consistent with the conclusion of the last section that the algorithm is stable.
It is interesting to note that the algorithm for stable even when the
hyperbolic transformations are applied directly. This is proven in [3]. The reason is
that whenever a hyperbolic transformation of large norm which might magnify errors
is performed, it can be proven that the norm of the generators drops drastically.
Since the error bounds, [4], for the unstabilized form of the hyperbolic transformation
are proportional to the norm of the new generators multiplied by the norm of the
transformation, the large norm of the hyperbolic transformation is canceled by a
proportional decrease in the norm of the generators.
The reason for the decrease in norm of the generators is that whenever ae is very
close to negative one, the action of a hyperbolic transformation is close to just taking
the difference between the two generators and scaling the result by 1=
. For
the case whenever the leading nonzero elements are O(j) apart, leading to a ae
very close to one, it can be proven that all the other components of the generators are
within O(j) of each other. As can be seen in the example here, this is clearly not the
case for more than two generators: the leading elements from which the hyperbolic
transformation is computed are within O(j), but the other components of the two
generators are only within O( p j). Their difference isn't small enough to completely
cancel out the large norm of the hyperbolic transformation.

Table
Errors for Different Values of j
In

Table

5, the backward error seems to display a proportionality to the square
root of the condition number. This is a property which is suggested by Theorem
4.1 and the fact that the overall backward error is just a sum of local errors. These
imply that neither the generators nor the effect of previous errors will be unduly
magnified by the hyperbolic transformations. This makes it possible to concentrate
our attention on just one stage of the algorithm. We do not expect the results in the
overall factorization to be significantly worse than the worst possible loss of accuracy
in a single stage.
Note that we have already proven that when applying a hyperbolic transforma-
tion, the result is always a set of generators which satisfy the bound
or, just looking at the two generators on which the hyperbolic transformation has
acted,
As a direct result of the analysis in [4], the local error introduced by the hyperbolic
transformation which produces u 1;k and v 1;k will be proportional to the norm of the
generators and the norm of the transformation, which can be bounded as follows
ip
F
Although there does not appear to be a general theory relating the condition number
of an arbitrary low ff-rank matrix to the values ae k , it is proven in [5] that for a positive
definite Toeplitz matrix,
Y
If similar inequalities hold for matrices of displacement rank ff, then this is enough to
suggest that the errors will tend to be proportional to the square root of the condition
We can get a somewhat more conclusive result. If we assume that the factorization
goes to completion, we can write of the leading elements of - u
which ae is computed,
We have assumed without loss of generality that - v
this is not the case, then either generator can be multiplied by \Gamma1
without causing any problems.
This means
So the error incurred during the hyperbolic transformation will be at worst proportional
top
ip
F
This explains the results in Table 5: because of the stability of the error propagation
in the algorithm, as well as the bound on the generators, the errors will be at
worst a modest multiple of p
ffl. Although the backward errors are dependent on the
reflection coefficients, there is a limit to this dependence.
6.

Summary

. This paper has proven the stability of a method for Cholesky factorization
of low ff-rank matrices. The stability result is dependent on the particular
implementation. Since the result depends on having a bound on the norm of the gener-
ators, it is important to choose transformations which permit such a bound. Working
with transformations which can be factored as a product of a sequence of plane rotations
and one hyperbolic transformation gives such a bound. In the case in which the
algorithm is implemented with unfactored hyperbolic transformations, this bound on
the generators is not sufficient for stability, however bounds were given which support
the assertion that the stability of the error propagation keeps the instability from
being as bad as might otherwise be expected.
Although not discussed here, it is worth noting that the approach in this paper
can be used in manner similar to that of [1] to provide a fast, weakly stable solution
of Toeplitz least squares problems. The conclusions would be essentially the same as
those in [1].



--R

A Weakly Stable Algorithm for General Toeplitz Systems

On the Stability of the Bareiss and Related Toeplitz Factorization Algorithms
A Note on Downdating the Cholesky Factorization
The Numerical Stability of the Levinson-Durbin Algorithm for Toeplitz Systems of Equations
High Performance Algorithms for Toeplitz and Block Toeplitz Matrices
Matrix Computations
Analysis of the Cholesky decomposition of a semi-definite matrix

Generalized Displacement Structure for Block-Toeplitz
The Algebraic Eigenvalue Problem
--TR

--CTR
Paolo Foschi , Erricos J. Kontoghiorghes, Estimation of VAR Models: Computational Aspects, Computational Economics, v.21 n.1-2, p.3-22, February -April
P. Alonso , J. M. Bada , A. M. Vidal, An Efficient Parallel Algorithm to Solve Block-Toeplitz Systems, The Journal of Supercomputing, v.32 n.3, p.251-278, June      2005
toeplitz matrices;structured matrices;stability;schur algorithm