On Positive Semidefinite Matrices with Known Null Space.
We show how the zero structure of a basis of the null space of a positive semidefinite matrix can be exploited to determine a positive definite submatrix of maximal rank. We discuss consequences of this result for the solution of (constrained) linear systems and eigenvalue problems. The results are of particular interest if A and the null space basis are sparse. We furthermore execute a backward error analysis of the Cholesky factorization of positive semidefinite matrices and provide new elementwise bounds.
Introduction
. The Cholesky factorization A
exists for any symmetric positive semidenite matrix A. In fact, R is the upper triangular
factor of the QR factorization of A 1=2 [11, x10.3]. R can be computed with the
well-known algorithm for positive denite matrices. However, zero pivots may appear.
As zero pivots come with a zero row/column in the reduced A, a zero pivot implies a
zero row in R. To actually compute a numerically stable Cholesky factorization of a
positive semidenite matrix one is advised to apply diagonal pivoting [11].
A semidenite matrix A may be given implicitly, in factored form
pn is of full row rank that does not need to be a
exposes the singularity of A explicitly as In this
case both the linear system and the eigenvalue problem can be solved eciently and
elegantly by working directly on the matrix F , never forming the matrix A explicitly.
In fact, in some applications, not assembling the matrix A but its factor F is the most
important step in the overall process of the numerical computation. One obvious
reason is that the (spectral) condition number of F is the square root of the condition
number of A. In nite element computation, F is the so called natural factor of the
stiness matrix A [2]. In the framework of linear algebra, every symmetric positive
semidenite matrix is the Gram matrix of some set of vectors, the columns of F .
Another possibility to have the singularity of A explicit is to have available a
basis of its null space N (A). This is the situation that we want to investigate in this
note. We will see that knowing a basis of N (A) allows to determine a priori when the
zero pivots will occur in the Cholesky factorization. It also permits to give a positive
denite submatrix of A right away. These results are of particular interest if A and the
null space basis are sparse. This is the case in the application from electromagnetics
that prompted this study [1]. There, a vector that is orthogonal to the null space
corresponds to a discrete electric eld that is divergence-free.
Our ndings permit to work with the positive denite part of A and to compute
a rank revealing Cholesky factorization A = R T R where the upper trapezoidal R
has full row rank. What is straightforward in exact arithmetic amounts to simply
replacing by zero potentially inaccurate small numbers. We analyze the error that is
introduced by this procedure.
Swiss Federal Institute of Technology (ETH), Institute of Scientic Computing, CH-8092 Zurich,
Switzerland (arbenz@inf.ethz.ch)
y University of Zagreb, Departement of Mathematics, Bijenicka 30, HR-10000 Zagreb, Croatia
(drmac@math.hr). The work of this author was supported by the Croatian Ministry of Science and
Technology grant 037012.
We complement this note with some implications of the above for solving eigenvalue
problems and constrained systems of equations.
2. Cholesky factorization of a positive semidenite matrix with known
null space. In this section we consider joint structures of a semidenite matrix A
and its null space.
Theorem 2.1. Let A = R T R be the Cholesky factorization of the positive
semidenite matrix A 2 R nn . Let Y 2 R nm with R(Y
m. These are the only zero entries on the diagonal of R.
Proof. Notice that the assumptions imply that Y := [y has full rank.
By Sylvester's law of inertia R has precisely m zeros on its diagonal. Further,
(Ry
whence r n i n i
If only n 1  n 2      nm , Y ,
ipped upside-down, can be transformed into
column-echelon form in order to obtain strong inequalities.
The Cholesky factor R appearing in Theorem 2.1 is an n  n upper triangular
matrix with m zero rows. These rows do not aect the product R T R. Therefore, they
can be removed from R to yield an (n m)  n matrix b
R with b
If the numbers n i are known, it is convenient to permute the rows of Y and
accordingly the rows and columns n i of A to the end. Then Theorem 2.1 can be
applied with m+i. The last m rows of R in Theorem 2.1 vanish. So, b
R is
upper trapezoidal.
After the just mentioned permutation the lowest mm block of Y is non-singular,
in fact, upper triangular. This consideration leads to an alternative formulation of
Theorem 2.1.
Theorem 2.2. Let A = R T R be the Cholesky factorization of the positive
semidenite matrix A 2 R nn . Let Y 2 R nm with R(Y (A). If the last
m rows of Y are linearly independent, then the leading principal (n m)  (n m)
submatrix of A is positive denite and R can be taken (n m)  n upper triangular.
Proof. Let
O Y 2
consists of the last m rows of Y . W is therefore invertible. Applying a congruence
transformation with W on A gives
A 11 A 12
A 21 A 22
O Y 2
A 11 O
O O
By Sylvester's law of inertia A 11 must be positive denite.
Let A
R 11 be the Cholesky factorization of A 11 . Then, the Cholesky
factor of the matrix in (2.2) is [R 11 ; O] 2 R (n m)n . Therefore, the Cholesky factor
of A is [R
Theorem 2.2 is applicable as long as the last m rows of Y form an invertible
matrix. of Y are linearly independent, we can permute Y such that
these rows become the last ones. In particular, if we want A 11 to be as sparse as
possible, we may choose to be the m most densely populated rows/columns
of A with the following greedy algorithm: If we have determined choose
i k+1 to be the index of the densest column of A such that rows
are linearly independent. In this way we can hope for an A 11 with sparse Cholesky
factors.
Remark 2.1. The equation
in
@
in a simply connected
domain
is satised by all constant functions u. The discretization
of (2.3) with nite elements of Lagrange type [4] leads to a positive semidenite
matrix A with a one dimensional null space spanned by the vector e with all entries
equal to 1. Theorem 2.1 now implies that, no matter how we permute A, in the
factorization the single zero on the diagonal of R will not appear before the
very last elimination step.
Example 2.1. Let A and Y be given by
O. As the last two rows of Y are linearly independent, Theorem 2.2 states
that the principal 33 submatrix of A is positive denite and that its Cholesky factor
is 3  5 upper triangular. In fact,
R =4
Let P be the permutation matrix, that exchanges 2nd with 4th and 3rd with 5th
entry of a 5-vector. Then,
Now we have according to Theorem 2.1 the Cholesky factor R 1
of A 1 has zero diagonal elements at positions 3 and 5. Indeed,
pp
3. Consistent semidenite systems. In this section we discuss how to solve
where A, R, and Y are as in Theorem 2.1. Without loss of generality, we can assume
that m. We split matrices and vectors in (3.1),
A 11 A 12
A T
R TR T
with is obtained from A by deleting rows and
columns m. The factorization (3.2) yields
Although A 11 is invertible, its condition number can be arbitrarily high. To reduce
ll-in during factorization [8] any symmetric permutations can be applied to A 11
without aecting the sequel. As R T has full rank, O or
diagonal
elements. Because the right side b of (3.1) has to satisfy
It is now easy to show that a particular solution of (3.1) is given by x with components
In fact, employing (3.3){(3.5) the second block row in (3.2) is
A T
The manifold S of the solutions of (3.1){(3.2) is
The vector a can be determined such that the solution x satises some constraints
In particular, if
then x is perpendicular to the null space of A.
Let now A be given implicitly as a Gram matrix
and Y 2 R nm be as above. (This may require renumbering the columns of F .) As
and as Y 2 is nonsingular, the block F 2 depends linearly on F 1 . Therefore, the QR
factorization of F has the form
R 11 R 12
O O
RR T , the factor equals the upper trapezoidal
Cholesky factor in (3.2).
4. Error Analysis. In this section we give backward error analyses for the
semidenite Cholesky factorization and for the null space basis.
4.1. Semidenite Cholesky factorization. The
oating-point computation
of the Cholesky factorization of a semidenite matrix is classied as unstable by
Higham [11, x10.3.2]. The principal problem is the determination of the rank of the
matrix.
If we assume, as we do in this note, that a basis of the null space of the matrix
under consideration is known a priori then, of course, its rank is known. Let A be
partitioned as in (3.2). We assume that A 11 2 R rr is positive denite numerically,
i.e. that the Cholesky factorization does not break down in
oating point arithmetic
with round-o unit u. Due to a result by Demmel [5] (see also [11, Thm.10.14]) this
is the case if,
min
where  min () denotes the minimal eigenvalue, kk is the spectral norm, and
If (4.1) does not hold, A 11 is not numerically denite. Note that
positive denite with unit diagonal. The assumption on  min can be relaxed
if, for instance, we use double precision accumulation during the factorization. Then
f(r) can be replaced by a small integer for all r not larger than 1=u. We assume,
however, that 2rf(r)u < 1.
The Cholesky decomposition of A is computed as indicated in (3.3). The Cholesky
factor of A 11 is computed rst. Then the matrix R 12 is obtained as the solution of
the matrix equation R T
Let e
R 11 denote the computed
oating-point Cholesky factor of A 11 . Then the
following two important facts are well-known.
(1) There exists a symmetric A 11 such that A 11
R
R 11 and
1i;jr
f(r)u: (4.2)
This is the backward error bound by Demmel [5], [11, Theorem 10.5].
(2) Let imply that the
Frobenius norm of
assumption (4.1) implies show [7] that there
exists an upper triangular matrix such that
e
s
s
<p:
Let e
R 12 be the
oating-point solution of the matrix equation e
R T
is the computed approximation of the exact Cholesky factor
Let e
R T e
R be partitioned conforming with (3.2). Since A+A is positive
semidenite and of rank r by construction, the equation e
A
A 1e
A 12 holds.
If we compute e
R 12 column by column, then, using Wilkinson's analysis of triangular
linear systems [11, Theorem 8.5],
R
6 PETER ARBENZ AND ZLATKO DRMA
where the matrix absolute values and the inequality are to be understood entry-wise.
Thus, we can write e
R 12 as
e
R T
R 12 j: (4.3)
Also, if we
R T
R T
11 j, we have
e
R T
R T
Further, from the inequality j e
using the M-matrix
property of
I
we obtain
Hence, relations (4.2), (4.3), (4.5) imply that the backward error for all (i; j) in the
in (3.2) is bounded by
A 22
We rst observe that k j j k
and that
s k.
Note that our assumptions imply that
r
2:
It remains to estimate the backward error in the (2; 2) block of the partition (3.2).
Using relation (4.4), we compute A 22 = e
R
R
R T
R 1e
R T
Using the inequalities from relations (4.4), (4.5) we obtain, for all (i; j),
We summarize the above analysis in the following
Theorem 4.1. Let A be a nn positive semidenite matrix of rank r with block
partition (3.2), where the r  r matrix A 11 is positive denite with the property (4.1).
Then the
oating-point Cholesky factorization with roundo u will compute an upper
trapezoidal matrix e
R of rank r such that e
R T e
A is a symmetric
backward perturbation with the following bounds:
A ii A jj
A ii A jj
~
A ii A jj ; r <
In the last estimate,
s k. Further, if e
is the exact Cholesky factor of A, then
e
R 11 R
R
~
Here, the matrix is upper triangular and  is to the rst order j j
.
Further, let the Cholesky factorization of A 11 be computed with pivoting so that
(R 11 ) ii
k=i (R 11
Then, the error R
R 11 R 11 is also
row-wise small, that is
Remark 4.1. Note that Theorem 4.1 also states that in the positive denite case
the Cholesky factorization with pivoting computes the triangular factor with small
column- and row-wise relative errors. This aects the accuracy of the linear equation
solver (forward and backward substitutions following the Cholesky factorization) not
only by ensuring favorable condition numbers but also by ensuring that the errors in
the coecients of the triangular systems are small.
4.2. Null space error. We now derive a backward error for the null space
Y of A. We seek an n  (n r) full rank matrix e
Y such that Y is
small and e
A e
As the null space and the range of A change simultaneously
(being orthogonal complements of each other), the size of Y necessarily depends on
a certain condition number of A; and the relevant condition number will depend on
the form of the perturbation A.
The equation that we investigate is e
equivalently, e
RY . If e
R is suciently close to R (to guarantee invertibility of e
RR
RR
RR
Though simple this equation is instructive. First of all, only the components of the
columns of R that lie in the null space N (A) aect the value of Y . Also, Y
keeps the full column rank of Y . Finally, Y T
kY k= min (Y ). It is easy to modify Y such that  min (Y )  1, e.g., if Y
Thus, kY k measures the angle between the true null space and the null space of the
perturbed matrix e
A. In the sequel we try to bound kY k.
If we rewrite (4.7) as
we get, after some manipulations,
Proposition 4.2. Let D be nonsingular matrix and let
If kR 0 (R
Here, y denotes the orthogonal projection onto the
null space of A.
We will discuss choices for D later. The Proposition indicates that the crucial
quantity for bounding kY k is kR 0 Y k. The following two examples detail this fact.
Example 4.1. Let  be big, of the order of 1=u, and let
pp
The null space of A is spanned by
which means that deleting
any row and column of A leaves a nonsingular 22 matrix. Let's choose it be the last
one, and let us follow the algorithm. For the sake of simplicity, let the only error be
committed in the computation of the (1; 1) entry of e
R 11 which is
instead of
3. Then we solve the lower triangular system for e
R 12 and obtain
e
Thus,
If we take perform the computation in Matlab where u  2:22  10
O() such that
the angle between Y and Y is small.
Example 4.2. We alter the (1; 1) entry
3 of R of the previous example to get ,
Again, we delete the last row and column of A and
proceed as in Example 4.1. Let us again assume that the only error occurs in the
entry of R 11 which becomes =(1
e
and
Again, O(1). But now also kY In fact, in computations with
Matlab, we observe an angle as large as O(10 2 ) between Y and Y .
Remark 4.2. Interestingly, if we set Example 4.1, the Matlab function
chol() computes the Cholesky factor
e
R =4
It is clear that the computed and stored A is a perturbation of the true A. Therefore,
numerically, it can be positive denite. It is therefore quite possible to know the rank
r < n of A exactly, to have a basis of the null space of A and a numerically stored
positive denite
oating-point A. Strictly speaking, this is a contradiction. Certainly,
from an application or numerical point of view, it is advisable to be very careful when
dealing with semideniteness.
In Examples 4.1 and 4.2 we excluded the largest diagonal entry of A. In fact, we
can give an estimate that relates the error in R 12 to the size of the deleted entries.
Suppose we managed the deleted diagonal entries of A to be the smallest
ones. Can we then guarantee that the relevant error in R will be small, and can we
check the stability by a simple, inexpensive test?
According to Theorem 4.1, the matrix R 11 is computed with row-wise small relative
error, provided that the Cholesky factorization of A 11 is computed with pivoting.
If that is the case, then it remains to estimate the row-wise perturbations of R 12 . If
is as in Theorem 4.1, then the inequality
holds for all
with some  i 2 (0; =2]. The angle  i has a nice interpretation. Let A = F T F be any
factorization of A, with full column rank and F T
Then  i is the angle between F 1 e i and the span of fF 1 e g. (This is easily
seen from the QR factorization of F 1 .)
The following Proposition states that well-conditioned dominance
of A 11 over A 22 ensure accurate rows of the computed matrix e
R.
Proposition 4.3. With the notation of Theorem 4.1, let A (and accordingly Y )
be arranged such that
If the Cholesky factorization of A 11 is computed with (standard) pivoting, then
sin  i
where sin  i is dened in (4.10).
Proof. This follows from relations (4.6), (4.9), (4.10) and the assumption (4.11).
We only note that in (4.9) and (4.10) we can replace
Remark 4.3. If spans N
partition of A s satises condition (4.11). If we apply the preceeding analysis to A s
and SY , we get an estimate for Y in the elliptic norm generated by S.
Note that Proposition 4.2 is true for any diagonal D as long as k(R 0 moderately
big and kR 0 k is small. We have just seen that R 0 is nicely bounded if we
choose
has an inverse nicely bounded
independent of A 11 because [11, x10]
Here the function h(r) is in the worst case dominated by 2 r and in practice one usually
observes an O(r) behaviour. In any case, k(D 1 R 11 is at most r times larger than
sophisticated pivoting can make sure that the behaviour of h(r)
is not worse than Wilkinson's pivot growth factor. We skip the details for the sake of
brevity.
To conclude, if the Cholesky factorization of A 11 is computed with pivoting and
relation (4.11) holds, then the backward error in Y can be estimated using (4.8)
and (4.12), where
4.3. Computation with implicit A. We consider now the backward stability
of the computation with A given implicitly as
rank r. Thus, the Cholesky factorization of A is accomplished by computing the QR
factorization of F .
In the numerical analysis of the QR factorization we use the standard, well-known
backward error analysis which can be found e.g. in [11, x18]. The simplest form of
this analysis states that the backward error in the QR factorization is column-wise
small. For instance, if we compute the Householder (or Givens) QR factorization of
F in
oating point arithmetic with roundo u, then the backward error F satises
n) is a polynomial of moderated degree in the matrix dimensions.
Our algorithm follows the same ideas as in the direct computation of R from
A. The knowledge of a null space basis admits that we can assume that F is in the
the p  r matrix F 1 is of rank r, see section 3. We then apply r
Householder re
ections to F which yields, in exact arithmetic, the matrix
R 11 R 12
O R 22
where R 11 2 R rr is upper triangular and nonsingular. If
conforming with F , then F is the QR factorization of F 1 .
In
oating point computation, R 22 is unlikely to be zero. Our algorithm simply
sets to zero whatever is computed as approximation of R 22 . As we shall see, the
backward error (in F ) of this procedure depends on a certain condition number of the
Theorem 4.4. Let F 2 R pn have rank r and be partitioned in the form
pr has the numerically well determined full rank r. More
specically, if is obtained from F 1 by scaling columns to have unit Euclidean
norm, then we assume that
Let the QR factorization of F be computed as described above, and let e
be the computed upper trapezoidal factor.
Then there exist a backward perturbation F and an orthogonal matrix b
Q such
that F
R is the QR factorization of F +F . The matrix F +F has rank r.
are partitioned as F , and Q 1 := b
then
c k;
e
R 11 R
bounds the roundo.
Proof. Let e
F (r) be the matrix obtained after r steps of the Householder QR
factorization. Then there exist an orthogonal matrix b
Q and a backward perturbation
F such that
e
R 11
e
R 12
O e
R 22
Our assumption on the numerical rank of F 1 implies that F 1
e
R 11 is the
QR factorization with nonsingular e
R 11 . Now, setting e
R 22 to zero is, in the backward
error sense, equivalent to the QR factorization of a rank r matrix,
R 11
e
R 12
O O
O O
O e
R 22
It remains to estimate b
First note that F
the i-th column of R 12 has the same norm as the corresponding column of F 2 . Then,
and we can write
To estimate Q 1 , we rst note that F
e
R 11 imply that
(R 11
e
and that
R
R
R
Thus, e
R
11 is the Cholesky factor of I +E, where
Now, by [7], kEkF < 1=2 implies that e
R
I +, where is upper triangular
and
<p:
Hence, R 11
e
Finally note that e
R 11 R
We remark that
e
which means that we can nicely bound R 12
R 12 R 12 . We have, for instance,
If we use entry-wise backward analysis of the QR factorization (jF
then we can also write
where the matrix absoulute values and inequalities are understood entry-wise, and " 2
is dened similarly as " 1 .
From the above analysis we see that the error in the computed matrix e
R is
bounded in the same way as in Theorem 4.1. Also, the QR factorization can be
computed with the standard column pivoting and R 11 can have additional structure
just as in the Cholesky factorization of A 11 . Therefore, the analysis of the backward
null space perturbation based on e
R T holds in this case as well. However, the bounds
of Theorem 4.4 are sharper than those of Theorem 4.1.
5. Constrained systems of equations. Let again be
Y 2 R nm having full rank. Let C 2 R nm be a matrix with full rank. Systems of
equations of the form
A C
x
y
c
appear at many occasions, e.g. in mixed nite element methods [3], or constrained
optimization [12]. They have a solution for every right side if R
which is the case if H := Y T C is nonsingular. In computations of Stokes [3] or
Maxwell equations [1] the second equation in (5.1) with imposes a divergence-free
condition on the
ow or electric eld, respectively.
To obtain a solution of (5.1) we rst construct a particular solution of the rst
block row. Pre-multiplying it by Y T yields As b Cy 2 R(A) we
can proceed as in section 3 to obtain a vector ~
x with A~ Cy. The solution x
of (5.1) is obtained by setting determining a such that C T
Thus,
x).
This procedure can be described in an elegant way if a congruence transformation
as in (6.2) is applied. Multiplying (5.1) by W T  I m , cf. (2.2), yields4
A 11 O C 1
O O H
a
n;r b:
Notice that ~
From (5.2) we read that
This geometric approach diers from the algebraic one based on the factorization4
A 11 A 12 C 1
A T
R T
O O
R T
O I m54
R
O O C
O C T
where the LU factorization of C
employed to solve (5.1). In the
geometric approach the LU factorization of H is used instead. Of course, there is a
close connection between the two approaches: Using (3.4) we get C T
. Notice that the columns of C or Y can be scaled such that the condition
numbers of H or C 2 R T
are not too big. Notice also that Y can be chosen
such that Y in which case C T
perturbation
analysis of (5.1){(5.3) remains to be done in our future work.
Golub and Greif [9] use the algebraic approach to solve systems of the form (5.1)
if the positive semidenite A has a low-dimensional null space. As they do not have
available a basis for the null space they apply a trial-and-error strategy for nding a
permutation of A such that the leading rr principal submatrix becomes nonsingular.
They report that usually the rst trial is successful. This is intelligible because n
the basis of the null space is dense which is often the case.
If the null space of A is high-dimensional then Golub and Greif use an augmented
Lagrangian approach. They modify (5.1) such that the (1; 1) block becomes positive
denite,
x
y
c
Here,  is some symmetric positive denite matrix, e.g. a multiple of the identity.
A+CC T is positive denite if Y T C is nonsingular. The determiniation of a good
is dicult. Golub and Grei thoroughly discuss how to choose  and how the 'penalty
aects the condition of the problem. In contrast to this approach where
a term is added to A that is positive denite on the null space of A, N (A) can be
avoided right away if a basis of it is known.
6. Eigenvalue problems. Let us consider the eigenvalue problem
where A is symmetric positive semidenite with
positive denite. We assume that the last m rows of Y are linearly independent such
that W in (2.1) is nonsingular. Then,
A 11 O
O O
where
Using the decomposition
O
O H
I O
I
with the Schur complement S := M 11 C
1 , and noting that P T W T
is easy to see that the positive eigenvalues of (6.1) are the eigenvalues of
A 11
14 PETER ARBENZ AND ZLATKO DRMA
Notice that S is dense, in general, whence, in sparse matrix computations, it should
not be formed explicitly.
If y is an eigenvector of (6.4) then
y
y
is an eigenvector of (6.1). By construction, C T
to the null space of A.
We now consider the situation when A and M are given in factored form,
such that the rank of F 1
equals the rank of A. Let us nd an implicit formulation of the reduced problem (6.4).
With W from (2.1) we have As before, A
R 11 is computed by the QR factorization of F 1 . It remains to compute a Cholesky
factor of the Schur complement S, but directly from the matrix B. To that end we
employ the QL factorization ('backward' QR factorization) of BW ,
whence, with (6.3),
Straightforward calculation now reveals that
Thus, the eigenvalues of the matrix pencil are the squares of the generalized
singular values [10] of the matrix pair (R equivalently, the squares of the
singular values of R 11 L 1
11 . An eigenvector y corresponds to a right singular vector
y. The blocks L 21 and L 22 come into play when the eigenvectors of (6.1) are to
be computed: using (6.7) equation (6.5) becomes
y
22 L 21 y
It is known that the GSVD of (R 11 ; L 11 ) can be computed with high relative
accuracy if the matrices (R 11 ) c and (L 11 ) c are well conditioned [6]. Here, (R 11 ) c and
are obtained by R 11 and L 11 , respectively, by scaling their columns to make
them of unit length. Obviously,  2 ((R 11 is the spectral
condition number. It remains to determine  2 ((L 11 ) c ). From (6.6) we get
whence
. Let the diagonal matrix D 1 be such that (B 1 1has columns of unit length. Further, let (B 1 be the QR factorization of
As Q 1 is orthogonal we have k(L 11
where  is the largest principal angle [10] between

r minD=diagonal  2 (L 11 D) [13] [11, Thm.7.5], we have
So, we have identied condition numbers that do not depend on column scalings and
that have a nice geometric interpretation. If the perturbations are column-wise small,
then these condition number are the relevant ones.
7. Concluding remarks. In this paper we have investigated ways to exploit the
knowledge of an explicit basis of the null space of a symmetric positive semidenite
matrix.
We have considered consistent systems of equations, constrained systems of equations
and generalized eigenvalue problems. First of all, the knowledge of a basis of the
null space of a matrix A permits to extract a priori a maximal positive semidenite
submatrix. The rest of the matrix is redundant information and is needed neither for
the solution of systems of equations nor for the eigenvalue computation. The order
of the problem is reduced by the dimension of the null space. In iterative solvers it is
not necessary to complement preconditioners with projections onto the complement
of the null space.
Our error analysis shows that a backward stable positive semidenite Cholesky
factorization exists if the principal r  r submatrix,
This does however not mean that the computed Cholesky factor ~
R has a null space
that is close to the known null space of R, A = R T R. We observed that the backward
error in the null space is small if the error in the Cholesky factor is (almost) orthogonal
to the null space of A. We show that this is the case if the positive denite principal
r r submatrix after scaling is well conditioned and if its diagonal elements dominate
those of the remaining diagonal block.
For systems of equations and eigenvalue problems, we considered the case when
is rectangular. This leads to interesting variants of the original
algorithms and most of all leads to more accurate results.
What remains to be investigated is the relation between extraction of a positive
denite matrix and ll-in during the Cholesky factorization. In future work we will
use the new techniques in applications and, if possible, extend the theory to matrix
classes more general than positive semidenite ones.



--R

A comparison of solvers for large eigenvalue problems originating from Maxwell's equations

Mixed and Hybrid Finite Element Methods
The Finite Element Method for Elliptic Problems
On oating point errors in Cholesky


Computer Solution of Large Sparse Positive De
Techniques for solving general KKT systems
Matrix Computations
Accuracy and Stability of Numerical Algorithms

Condition numbers and equilibration of matrices
--TR
positive semidefinite matrices;null space basis;cholesky factorization