Perturbation Analyses for the Cholesky Downdating Problem.
New perturbation analyses are presented for the block Cholesky downdating problem These show how changes in R and X alter the Cholesky factor U. There are two main cases for the perturbation matrix $\D R$ in R: (1) $\D R$ is a general matrix; (2)$\D R$ is an upper triangular matrix. For both cases, first-order perturbation bounds for the downdated Cholesky factor U are given using two approaches --- a detailed "matrix--vector equation" analysis which provides tight bounds and resulting true condition numbers, which unfortunately are costly to compute, and a simpler "matrix equation" analysis which provides results that are weaker but easier to compute or estimate. The analyses more accurately reflect the sensitivity of the problem than previous results. As $X\rightarrow 0$, the asymptotic values of the new condition numbers for case (1) have bounds that are independent of $\kappa_2(R)$ if $R$ was found using the standard pivoting strategy in the Cholesky factorization, and the asymptotic values of the new condition numbers for case (2) are unity. Simple reasoning shows this last result must be true for the sensitivity of the problem, but previous condition numbers did not exhibit this.}
Introduction
. Let A 2 R n\Thetan be a symmetric positive definite matrix. Then
there exists a unique upper triangular matrix R 2 R n\Thetan with positive diagonal elements
such that A = R T R. This factorization is called the Cholesky factorization,
and R is called the Cholesky factor of A.
In this paper we give perturbation analyses of the following problem: given an
upper triangular matrix R 2 R n\Thetan and a matrix X 2 R k\Thetan such that R T
is positive definite, find an upper triangular matrix U 2 R n\Thetan with positive diagonal
elements such that
This problem is called the block Cholesky downdating problem, and the matrix U
is referred to as the downdated Cholesky factor. The block Cholesky downdating
problem has many important applications, and the case for k=1 has been extensively
studied in the literature (see [1, 2, 3, 8, 11, 12, 16, 17, 18]).
Perturbation results for the single Cholesky downdating problem were presented
by Stewart [18]. Eld'en and Park [10] made an analysis for block downdating. But
these two papers just considered the case that only R or X is perturbed. More
complete analyses, with both R and X being perturbed, were given by Pan [15] and
Sun [20]. Pan [15] gave first order perturbation bounds for single downdating. Sun [20]
gave strict, also first order perturbation bounds for single downdating and first order
perturbation bounds for block downdating.
The main purpose of this paper is to establish new first order perturbation results
and present new condition numbers which more closely reflect the true sensitivity of
This research was partially supported by NSERC of Canada Grant OGP0009236.
y School of Computer Science, McGill University, Montreal, Quebec, Canada, H3A 2A7,
(chang@cs.mcgill.ca), (chris@cs.mcgill.ca).
the problem. In Section 2 we will give the key result of Sun [20], and a new result using
the approach of these earlier papers. In Section 3 we present new perturbation results,
first by the straightforward matrix equation approach, then by the more detailed and
tighter matrix-vector equation approach. The basic ideas behind these two approaches
were discussed in Chang, Paige and Stewart [6, 7]. We give numerical results and
suggest practical condition estimators in Section 4.
Previous papers implied the change \DeltaR in R was upper triangular, and Sun [20]
said this, but neither he nor the others made use of this. In fact a backward stable
algorithm for computing U given R and X would produce the exact result U
for nearby data R+ \DeltaR and X it is not clear that \DeltaR would be upper
triangular - the form of the equivalent backward error \DeltaR would depend on the
algorithm, and if it were upper triangular, it would require a rounding error analysis
to show this. Thus for completeness it seems necessary to consider two separate cases
- upper triangular \DeltaR and general \DeltaR. We do this throughout Sections 3-4, and
get stronger results for upper triangular \DeltaR than in the general case.
In any perturbation analysis it is important to examine how good the results
are. In Section 3.2 we produce provably tight bounds, leading to the true condition
numbers (for the norms chosen). The numerical example in Section 4 indicates how
much better the results of this new analysis can be compared with some earlier ones,
but a theoretical understanding is also desirable. By considering the asymptotic case
as 0, the results simplify, and are easily understandable. We show the new
results have the correct properties as X ! 0, in contrast to earlier results.
Before proceeding, let us introduce some notation. Let n\Thetan , then
up(B), sut(B), slt(B) and diag(B) are defined by
and
2. Previous perturbation results, and an improvement. The condition
number for downdating presented by Pan [15] involves the square of the condition
number of R, - 2
proposed new condition
numbers which are simple and proportional to - 2 (R). The condition number for the
block Cholesky downdating problem proposed in [20] is
(\Gamma) is the smallest singular value
of \Gamma. Notice that for fixed 1. Now we use a
similar approach to derive a new bound, from which Sun's bound follows.
First we derive some relationships among U , R, X and \Gamma.
From (1.1) obviously we have
From (1.1) it follows that
so that taking the 2-norm gives
From (1.1) we have
which, combined with (2.2), gives
From (1.1) we have
which, combined with (2.2), gives
s
By (2.2) we have
Finally from (2.4) we see
To derive first order perturbation results we consider the perturbed version of
where U , U+ \DeltaU and R are upper triangular matrices with positive diagonal elements.
when \DeltaR and \DeltaX are sufficiently small, (2.7) has a unique solution \DeltaU .
Multiplying out the two sides of (2.7) and ignoring second order terms, we obtain a
linear matrix equation for the first order approximation d
\DeltaU to \DeltaU :
U T d
\DeltaU T
In fact it is straightforward to show d
U(0), the rate of change of U(-) at
so d
\DeltaU also has a precise meaning. From (2.8), we have
d
(R
Notice since d
\DeltaU U \Gamma1 is upper triangular, it follows, with (1.2), that
d
(R
But for any symmetric matrix B,
F \Gamma2 (b 2
Thus from (2.9) we have
\DeltaU
pkU \GammaT (R T \DeltaR
which, combined with (2.2) and (2.3), gives
\DeltaU k F -
resulting in the new perturbation bound for relative changes
\DeltaU k F
pp
which leads to the condition numbers for the Cholesky downdating problem:
pp
for U with respect to relative changes in R and X, respectively. Notice from (2.1)
. So we can define a new overall condition number
Rewriting (2.10) as
\DeltaU
and combining it with (2.5) and (2.6), gives Sun's bound
\DeltaU k F
We have seen the right hand side of (2.10) is never worse than that of (2.12), so
Although fi 2 is a minor improvement on fi 1 , it is still not what we want. We can
see this from the asymptotic behavior of these "condition numbers". The Cholesky
factorization is unique, so as X ! 0, U ! R, and X T \DeltaX ! 0 in (2.8). Now for any
upper triangular perturbation \DeltaR in R, \DeltaU ! \DeltaR, so the true condition number
should approach unity. Here
(R). The next section shows how we can
overcome this inadequacy.
3. New perturbation results. In Section 2 we saw the key to deriving first
order perturbation bounds for U in the block Cholesky downdating problem is the
equation (2.8). We will now analyze it in two new approaches. The two approaches
have been used in the perturbation analyses of the Cholesky factorization, the QR factorization
(see Chang, Paige and Stewart [6, 7]), and LU factorization (see Chang [4]
and Stewart [19]). The first approach, the refined matrix equation approach, gives a
clear improvement on the previous results, while the second, the matrix-vector equation
approach, gives a further improvement still, which leads to the true condition
numbers for the block Cholesky downdating problem.
3.1. Refined matrix equation analysis. In the last section we used (2.8) to
produce the matrix equation (2.9), and derived the bounds directly from this. We
now look at this approach more closely.
Let D n be the set of all n \Theta n real positive definite diagonal matrices. For any
U . Note that for any matrix B we have
First with no restriction on \DeltaR we have from (2.9)
d
U
so taking the F-norm gives
\DeltaU
It is easy to show for any B 2 R n\Thetan (see Lemma 5.1 in [7])
g. Thus from (3.1) we have
\DeltaU
(kU \GammaT R T \DeltaR -
U
(using (2:2); (2:3))
which is an elegant result in the changes alone. It leads to the following perturbation
bound in terms of relative changes
\DeltaU
Although here it would be simpler to just define an overall condition number, for later
comparisons it is necessary for us to define the following two quantities as condition
numbers for U with respect to relative changes in R and X, respectively (here subscript
G refers to general \DeltaR, and later the subscript T will refer to upper triangular \DeltaR):
c RG (R; X)
where
c RG (R; X;D) j
Then an overall condition number can be defined as
c G (R; X;D);
where
Obviously we have
c G (R;
Thus with these, we have from (3.3) that
\DeltaU k F
if we take become (2.10), and
It is not difficult to give an example to show fi 2 can be arbitrarily larger than c G (R; X),
as can be seen from the following asymptotic behaviour.
It is shown in [7, x5.1, (5.14)] that with an appropriate choice of D,
has a bound which is a function of n only, if R was found using the standard pivoting
strategy in the Cholesky factorization, and in this case, we see the condition number
c G (R; X) of the problem here is bounded independently of - 2 (R) as X ! 0, for general
\DeltaR. At the end of this section we give an even stronger result when X ! 0 for the
case of upper triangular \DeltaR. Note in the case here that fi 2 in (2.11) can be made as
large as we like, and thus arbitrarily larger than c G (R; X).
In the case where \DeltaR is upper triangular , we can refine the analysis further. From
we have
d
Notice with (1.3) and (1.4)
U \GammaT R T \DeltaRU
But for any upper triangular matrix T we have
so that if we define T
up[diag(U \GammaT R T
Thus from (3.12), (3.13) and (3.14) we obtain
d
\DeltaU
As before, let
U , where . From (3.15) it follows that
\DeltaU
U \GammaT \DeltaR T
Then, applying (3.2) to this, we get the following perturbation bound
\DeltaU k F
(3.
Comparing (3.16) with (3.3) and noticing (2.3), we see the sensitivity of U with respect
to changes in X does not change, so c X (R; X) defined in (3.4) can still be regarded as
a condition number for U with respect to changes in X. But we now need to define
a new condition number for U with respect to upper triangular changes in R, that is
(subscript T indicates upper triangular \DeltaR)
c RT (R; X)
c RT (R; X;D);
where
Thus an overall condition number can be defined as
where
c T (R;
Obviously we have
With these, we have from (3.16) that
\DeltaU k F
What is the relationship between c T (R; X) and c G (R;
n \Theta n upper triangular matrix observe the following two facts:
are the eigenvalues of T , so that
which gives
(Note: In fact we can prove a slightly sharper inequality ksut(T
Therefore
c RT (R;
n)
n)
using (2:2))
n)c RG (R; X;D);
so that
c RT (R; X)
n)c RG (R; X):
Thus we have from (3.8) and (3.18)
n)c G (R; X):
On the other hand, c T (R; X) can be arbitrarily smaller than c G (R; X). This can be
seen from the asymptotic behaviour, which is important in its own right. As
since
so for upper triangular changes in R, whether pivoting was used in finding R or not,
Thus when X ! 0, the bound in (3.19) reflects the true sensitivity of the problem.
For the case of general \DeltaR, if we do not use pivoting it is straightforward to make
c G (R; X) in (3.7) arbitrarily large even with
3.2. Matrix-vector equation analysis. In the last subsection, based on the
structure of \DeltaR, we gave two perturbation bounds using the so called refined matrix
equation approach. Also based on the structure of \DeltaR, we can now obtain provably
sharp, but less intuitive results by viewing the matrix equation (2.8) as a large matrix-vector
equation. For any matrix C n\Thetan , denote by c (i)
j the
vector of the first i elements of c j . With this, we define ("u" denotes "upper")
c (1)c (2):
c (n)
It is the vector formed by stacking the columns of the upper triangular part of C into
one long vector.
First assume \DeltaR is a general real n \Theta n matrix. It is easy to show (2.8) can be
rewritten into the following matrix-vector form (cf [7])
WU uvec( d
2 \Theta n(n+1)
r 11
and YX 2 R n(n+1)
Since U is nonsingular, WU is also, and from (3.22)
uvec( d
U YX vec(\DeltaX);
so taking the 2-norm gives
\DeltaU
resulting in the following perturbation bound
\DeltaU k F
-RG (R; X)
where
Now we would like to show
Before showing this, we will prove a more general result. Suppose from (2.8) we are
able to obtain a perturbation bound of the form
\DeltaU
- ff R
(3.
where ff R and ff X , two functions of R and X, are other measures of the sensitivity of
the Cholesky downdating problem with respect to changes in R and X. Let
Then from (3.23) and (3.28) we have
U ZR vec(\DeltaR)k 2
- ff R
Notice \DeltaR can be any (sufficiently small) n \Theta n real matrix, so we must have
which gives
Similarly, we can show
Notice since (3.9) is a particular case of (3.28), (3.27) follows. Thus we have from
(3.8) and (3.26)
The above analysis shows for general \DeltaR, -RG (R; X) and -X (R; X) are optimal
measures of the sensitivity of U with respect to changes in R and X, respectively, and
thus the bound (3.24) is optimal. So we propose -RG (R; X) and -X (R; X) as the true
condition numbers for U with respect to general changes in R and X, respectively,
and - CDG (R; X) as the true overall condition number of the problem in this case.
It is easy to observe that if X ! 0, - CDG (R; X)
just WU with each entry u ij replaced by r ij . If R was found using the standard
pivoting strategy in the Cholesky factorization, then kW \Gamma1
R ZR k 2 has a bound which
is a function of n alone (see [5] for a proof). So in this case our condition number
- CDG (R; X) also has a bound which is a function of n alone as
Remark 1: Our numerical experiments suggest c G (R; X) is usually a good approximation
to - CDG (R; X). But the following example shows c G (R; X) can sometimes be
arbitrarily larger than - CDG (R; X).
where ffl is a small positive number. It is not difficult to show
But c G (R; X) has an advantage over - CDG (R; X) - it can be quite easy to estimate
- all we need do is choose a suitable D in c G (R; X;D). We consider how to do this
in the next section. In contrast - CDG (R; X) is, as far as we can see, unreasonably
expensive to compute or estimate.
Now we consider the case where \DeltaR is upper triangular. (2.8) can now be rewritten
as the following matrix-vector form
WU uvec( d
2 \Theta n(n+1)
2 and YX 2 R n(n+1)
2 \Thetakn are defined as before, and WR 2
R n(n+1)
2 \Theta n(n+1)
2 is just WU with each entry u ij replaced by r ij . Since U is nonsingular,
WU is also, and from (3.30)
uvec( d
U YX vec(\DeltaX);
so taking the 2-norm gives
\DeltaU
which leads to the following perturbation bound
\DeltaU k F
where
Note -X (R; X) is the same as that defined in (3.25).
As before, we can show that for the case where \DeltaR is upper triangular, - RT (R; X)
and -X (R; X) are optimal measures of the sensitivity of U with respect to changes in
R and X, respectively, and thus the bound (3.32) is optimal. In particular, we have
In fact - RT (R; X) -RG (R; X) can also be proved directly by the fact that the columns
of WR form a proper subset of the columns of ZR , and the second inequality has been
proved before. Thus we have from (3.8), (3.18), (3.26) and (3.33)
By the above analysis, we propose - RT (R; X) and -X (R; X) as the true condition
numbers for U with respect to changes in R and X, respectively, and - CDT (R; X) as
the true overall condition number, in the case that \DeltaR is upper triangular.
If as well X ! 0, then since U ! R, W \Gamma1
So in this case the Cholesky downdating problem becomes very well conditioned no
matter how ill conditioned R or U is.
Remark 2: Numerical experiments also suggest c T (R; X) is usually a good approximation
to - CDT (R; X). But sometimes c T (R; X) can be arbitrarily larger than
- CDT (R; X). This can also be seen from the example in Remark 1. In fact, it is not
difficult to obtain
Like - CDG (R; X), - CDT (R; X) is difficult to compute or estimate. But c T (R; X) is easy
to estimate, which is discussed in the next section.
4. Numerical tests and condition estimators. In Section 3 we presented new
first order perturbation bounds for the the downdated Cholesky factor U using first
the refined matrix equation approach, and then the matrix vector equation approach.
We defined - CDG (R; X) for general \DeltaR, and - CDT (R; X) for upper triangular \DeltaR, as
the true overall condition numbers of the problem. Also we gave two corresponding
practical but weaker condition numbers c G (R; X) and c T (R; X) for the two \DeltaR cases.
We would like to choose D such that c G (R; X;D) and c T (R; X;D) are good approximations
to c G (R; X) and c T (R; X), respectively. We see from (3.5), (3.6) and
(3.17) that we want to find D such that
its infimum.
By a well known result of van der Sluis [21], - 2 (D \Gamma1 U) will be nearly minimal when
the rows of D \Gamma1 U are equilibrated. But this could lead to a large i D . So a reasonable
compromise is to choose D to equilibrate U as far as possible while keeping i D - 1.
Specifically, take
we use a standard condition estimator to estimate
Notice from (2.4) we have oe n
2 . Usually k, the number of
rows of X, is much smaller than n, so oe n (\Gamma) can be computed in O(n 2 ). If k is not
much smaller than n, then we use a standard norm estimator to estimate kXR
in O(n 2 ). Similarly kUk 2 and kRk 2 can be estimated in O(n 2 ). So finally c G (R; X;D)
can be estimated in O(n 2 ). Estimating c T (R; X;D) is not as easy as estimating
c G (R; X;D). The part kdiag(RU (R; X;D) can easily be computed in
O(n), since diag(RU
c RT (R; X;D) can roughly be estimated in O(n 2 ), based onp
F
and the fact that kRU \Gamma1 k F can be estimated by a standard norm estimator in O(n 2 ).
The value of kXU (R; X;D) can be calculated (if k !! n) or estimated
by a standard estimator in O(n 2 ). All the remaining values kRk 2 , kXk 2 and kUk 2
can also be estimated by a standard norm estimator in O(n 2 ). Hence c RT (R; X;D),
c X (R; X;D), and thus c T (R; X;D) can be estimated in O(n 2 ). For standard condition
estimators and norm estimators, see Chapter 14 of [14].
The relationships among the various overall condition numbers for the Cholesky
downdating problem presented in Section 2 and Section 3 are as follows.
n)c G (R; X)
Now we give one numerical example to illustrate these. The example, quoted from
Sun [20], is as follows.
0:240
2:390C C C C C A
. The results obtained using MATLAB are shown in

Table

4.1 for various values of - :

Table
c G (R; X;D) 3.60e+03 3.61e+02 3.79e+01 1.79e+01 1.78e+01 1.78e+01
c T (R; X;D) 2.12e+03 2.12e+02 1.79e+01 1.07e+00 1.00e+00 1.00e+00
Note in Table 4.1 how fi 1 and fi 2 can be far worse than the true condition numbers
- CDG (R; X) and - CDT (R; X), although fi 2 is not as bad as fi 1 . Also we observe
that c G (R; X;D) and c T (R; X;D) are very good approximations to - CDG (R; X)
and - CDT (R; X), respectively. When X become small, all of the condition numbers
decrease. The asymptotic behavior of c G (R; X;D), c T (R; X;D), - CDG (R; X) and
our theoretical results - when
- CDG (R; X) will be bounded in terms of n since here R corresponds to the Cholesky
factor of a correctly pivoted A, and c T (R; X); - CDT (R; X) ! 1.

Acknowledgement

. We would like to thank Ji-guang Sun for suggesting to us
that the approach describe in (a draft version of) [6] might apply to the Cholesky
downdating problem.



--R

Analysis of a recursive least squares hyperbolic rotation algorithm for signal processing
Accurate downdating of least squares solutions
A note on downdating the Cholesky factorization
PhD Thesis
A perturbation analysis for R in the QR factorization
New perturbation analyses for the Cholesky factorization
Perturbation analyses for the QR factorization

Block downdating of least squares solutions
Perturbation analysis for block downdating of a Cholesky decomposition
Numerical computations for univariate linear models
Methods for modifying matrix factor- izations
Matrix computations
Accuracy and Stability of Numerical Algorithms
A perturbation analysis of the problem of downdating a Cholesky factorization
Least squares modification with inverse factorizations: parallel implications

The effects of rounding error on an algorithm for downdating a Cholesky factor- ization
On the Perturbation of LU and Cholesky Factors
Perturbation analysis of the Cholesky downdating and QR updating problems
Condition numbers and equilibration of matrices
--TR
asymptotic condition;downdating;cholesky factorization;perturbation analysis;sensitivity;condition