Parallel Complexity of Numerically Accurate Linear System Solvers.
We prove a number of negative results about practical (i.e., work efficient and numerically accurate) algorithms for computing the main matrix factorizations. In particular, we prove that the popular Householder and Givens methods for computing the QR decomposition are P-complete, and hence presumably inherently sequential, under both real and floating point number models. We also prove that Gaussian elimination (GE) with a weak form of pivoting, which aims only at making the resulting algorithm nondegenerate, is likely to be inherently sequential as well. Finally, we prove that GE with partial pivoting is P-complete over GF(2) or when restricted to symmetric positive definite matrices, for which it is known that even standard GE (no pivoting) does not fail. Altogether, the results of this paper give further formal support to the widespread belief that there is a tradeoff between parallelism and accuracy in numerical algorithms.
Introduction
strongly nonsingular
Mauro Leoncini Giovanni Manzini Luciano Margara
August 8, 1997
Parallel Complexity of Numerically Accurate Linear
System Solvers.
This work merges preliminary results presented at ESA '96 and SPAA '97.
Dipartimento di Informatica, Universit'a di Pisa, Corso Italia 40, 56125 Pisa, Italy, and IMC-CNR.
via S. Maria 46, 56126 Pisa, Italy. Email: leoncini@di.unipi.it. Supported by Murst 40% funds.
Dipartimento di Scienze e Tecnologie Avanzate, Universit'a di Torino, Via Cavour 84, 15100 Alessan-
dria, Italy. Email: manzini@unial.it.
Dipartimento Scienze dell'Informazione, Universit'a di Bologna, Piazza Porta S. Donato 5, 40127
Italy. Email: margara@cs.unibo.it.
Matrix factorization algorithms form the backbone of state-of-the-art numerical libraries
and packages, such as LAPACK and MATLAB [2, 14]. Indeed, factoring a matrix is
almost always the first step of many scientific computations, and usually the one which
places the heaviest demand in terms of computing resources. In view of their importance,
some authors have investigated the parallel complexity of the most popular matrix factor-
izations, namely the ( ) and (5) decompositions (see Appendix A for definitions
and simple properties). A list of positive known results follows.
decomposition is in arithmetic , whenever it exists, i.e., provided that the
leading principal minors of the input matrix are nonsingular (in this case we will
say that the matrix is ) [16, 18].We prove a number of negative results about practical (i.e., work efficient and
numerically accurate) algorithms for computing the main matrix factorizations. In
particular, we prove that the popular Householder and Givens' methods for computing
the QR decomposition are P-complete, and hence presumably inherently
sequential, under both real and floating point number models. We also prove that
Gaussian Elimination (GE) with a weak form of pivoting, which only aims at making
the resulting algorithm nondegenerate (but possibly unstable), is likely to be
inherently sequential as well. Finally, we prove that GE with partial pivoting is
P-complete when restricted to Symmetric Positive Definite matrices, for which it
is known that even plain GE does not fail. Altogether, the results of this paper
give further formal support to the widespread belief that there is a tradeoff between
parallelism and accuracy in numerical algorithms.
ffl2 1log
stable
Minimal
LU
PLU NC
QR A NC
LU QR
LU QR
PLU
O n
O n O n
O n n O
Not to be confused with the analogous LFMIS problem of graph theory, which is known to be P-complete
[10].
decomposition is in arithmetic for matrices with full column rank, since it
easily reduces to decomposition of strongly nonsingular matrices [16].
decomposition is in arithmetic for nonsingular matrices [7]. The algorithm
for finding a permutation matrix such that is strongly nonsingular builds
upon the computation of the Lexicographically First Maximal Independent Subset
of the rows of a matrix, which is in [3] .
5 factorization of an arbitrary matrix is in arithmetic [7]. A permutation
5 such that the leftmost submatrix of has full column rank,
can be found by computing LFMIS of sets of (column) vectors.
Unfortunately, none of the above algorithms has proved to be numerically accurate
with respect to a realistic model of arithmetic, such as double precision floating point.
Actually, finding a numerically stable NC algorithm to compute the (or ) decomposition
of a matrix can be regarded as one of the major open problems in parallel
computation theory [10].
That a positive solution to the above problem is not just around the corner is confirmed
by the negative results that can be proved of the algorithms in practical use for
computing the and decompositions. Already in 1989 Vavasis proved that Gaussian
Elimination with Partial pivoting (GEP), which is the standard method for computing
the decomposition, is P-complete over the real or rational numbers [20]. Note that,
strictly speaking, membership in P could not be defined for the real number model. When
dealing with real matrices and real number computations we then assume implicitly that
the class P be defined to include the problems solvable in polynomial time on such models
as the real RAM (see [17]). The result in [20] was proved by showing that a decision
problem defined in terms of GEP's behavior was P-complete. For parallel complexity
theory the P-completeness result implies that GEP is likely to be inherently sequential,
i.e., admitting no implementations . One of the authors proved then
that GEP is probably even harder to parallelize, in the sense that no (
implementation can exist unless all the problems in admit polynomial speedup [12].
In this paper we prove new negative results about the classical factorization algo-
rithms. We consider the methods of Householder's reflections and of Givens' rotations
to compute the decomposition of a matrix. The main use of the decomposition
is within iterative methods for the computation of the eigenvalues of a matrix and to
compute least squares solutions to overdetermined linear systems [9]. Moreover, both
Householder's and Given's methods (hereafter denoted HQR and GQR, respectively) are
potential competitors of Gaussian Elimination to solve systems of linear equations stably
in parallel. To date, the fastest stable parallel solver is based on GQR and is characterized
by ( ) parallel time on an ( ) processor PRAM [19], while both GEP and HQR run
in ( log processors. Also, GQR is especially suitable for solving
large sparse systems, given its ability to annihilate selected entries of the input matrix at
very low cost.
We also consider the application of Gaussian Elimination with Partial pivoting to
special classes of matrices and a weaker form of pivoting which we will call
R
pivoting (GEM). Under minimal pivoting, the pivot chosen to annihilate a given column
is the first nonzero on or below the main diagonal. Minimal pivoting is especially suitable
for systolic-like implementations of linear system solvers (see, e.g., [11], although it is not
called this way). Minimal pivoting can be regarded as the minimum modification required
for Gaussian Elimination to be nondegenerate on arbitrary input matrices.
We prove the following results.
1. HQR and GQR are P-complete over the real or floating point numbers. We exhibit
reductions from the NAND Circuit Value Problem (NANDCVP) with fanout 2.
In particular, what we prove to be P-complete is to decide the sign of a given diagonal
element of the upper triangular matrices computed by either HQR or GQR. Our
reductions seem to be more intricate than the simple one in [20]. This is probably
a consequence of the apparently more complex effect of reflections and rotations
with respect to the linear combinations of Gaussian Elimination. We would like to
stress that the P-completeness proofs for the case of floating point arithmetic apply
directly, and have been checked with, the algorithms available in the state-of-the-art
package Matlab using the IEEE 754 standard for floating point arithmetic. In other
words, the negative results apply to widely "in-use" algorithms.
2. We extend Vavasis' result proving that GEP is P-complete on input strongly non-singular
matrices. This class includes matrices which are important in practical
applications, namely the diagonally dominant and symmetric positive definite ones.
Note that plain Gaussian Elimination (no pivoting) is guaranteed not to fail on
input a strongly nonsingular matrix. However, since it is usually unstable, one still
uses GEP.
3. We prove that GEM is P-complete on general matrices.
4. We show that the known algorithm for computing a PLU decomposition of a
nonsingular matrix corresponds to GE with a nonstandard pivoting strategy which
only slightly differs from Minimal Pivoting. Also, we prove that GE with such a
nonstandard strategy is P-complete on input arbitrary matrices, which somehow
accounts for the difficulties of finding an algorithm to compute the PLU decomposition
of possibly singular matrices.
The results of this paper give further evidence of the pervasiveness of a phenomenon
that has been observed also by numerical analysts (from a more practical perspective).
Namely that there is a "tradeoff" between the degree of parallelism, on the one hand, and
nondegeneracy and accuracy properties of numerical algorithms, on the other [5].
The rest of this paper is organized as follows. In Section 2 we introduce a little notation
and give some preliminary definitions. In Section 3 we describe the key ideas that are
common to the P-completeness proofs for all the factorization methods considered in the
paper. In Sections 4 and 5 we address QR decomposition via Householder's reflections and
Givens' rotations, respectively. In Section 6 we prove our negative results about Gaussian
Elimination with Partial and Minimal pivoting. In section 7 we show a correspondence
between a known decompostion algorithm and Gaussian Elimination. We
conclude with some further considerations and open problems. In Appendix A we discuss
Preliminaries
F
column
minor
principal
row
row echelon
a A
x y x
A a a
A A A
I O I i j
O
x
A n
LU A L; U L
A LU A LU
A
PLU A
PLU
QR A Q; R Q
R A QR QR
the algorithms considered in this paper. In Appendix B we give some basic definitions
about the floating point number representation. Clearly, more material about these well
known algorithms and the computer arithmetic can be found in many excellent textbooks
(in particular, see [9]). Finally, we include one technical proof in Appendix C.
The notations adopted here for matrices and matrix-related concepts are the standard
ones (see [9]). Matrices are denoted by capital letters. The ( ) entry of a matrix
is referred to by either or [ ] . Vectors are designated by lower case letters, usually
taken from the end of the alphabet, e.g., , , etc. Note that the notation refers to a
vector, i.e., an 1 matrix, for some 1.
The -th row (resp., column) of a matrix is denoted by (resp. A of a
matrix is any submatrix of . A minor is any square submatrix of formed
by the same set of row and column indices.
The symbols and denote the identity matrix (with [
otherwise) and the zero matrix (such that [ respectively. The zero vector is
denoted using the symbol 0. The transpose of is the matrix such that = . is
usually denoted by . A matrix is orthogonal when = . A permutation matrix
is a matrix which is zero everywhere except for just one 1 in each row and column.
Any permutation matrix is orthogonal. The transpose of a (column) vector is the
vector , i.e., a matrix of size 1 , for some .
Let be a square matrix of order .
The decomposition of is a pair of matrices , such that is lower
diagonal elements, is upper triangular,
. For an arbitrary (even nonsingular) matrix the decomposition
might not be defined. A sufficient condition for its existence (and unicity) is that
be strongly nonsingular.
The decomposition of is a triple of matrices such that and are
as above, is a permutation matrix,
The decomposition is always defined but not unique.
The decomposition of is a pair of matrices , such that is orthogonal,
is upper triangular, and = . The decomposition always exists.
In all the above cases, if is , with , and when the factorization exists, we get
a matrix that is properly said to be in form (rather than upper triangular).
Its leftmost minor is upper triangular while its rightmost ( ) minor is in
general a dense submatrix. However, when no confusion is possible, we will always speak
of the triangular factor of a given factorization.
A detailed description of the algorithms considered in this paper can be found in


Appendix

A. The details, however, are not necessary to understand the common structure
of the reductions. Some details will be required in the proofs of Theorems 4.3 and 5.3,
which deal with the floating point version of the algorithms. Except for these, the
following general description is sufficient. It defines a class of matrix factorization
F
proper embedding
A
p3
Input
Output
3 A framework for reductions to
kk
A A
A A
k a a
a a i a a k
A
A A A A k n A k k
A
A A
O B
A A
O B
R A A
A A a
algorithms that includes, among the others, the classical QR algorithms and Gaussian
Elimination.
Let be the input matrix. The algorithms in bring to upper triangular form by
applying a series of transformations that introduce zeros in the strictly lower triangular
portion of , from left to right. The notation is usually adopted to indicate the
matrix obtained after 1 transformations and its elements are referred to by .
is zero for min . A transformation is applied during one of the algorithm.
Every algorithm satisfies the following properties.
In other words, if the first
entries in row are zero, then the first stages of do not modify row .
If column has complementary nonzero structure with respect to columns 1 1
(by which we mean
not affected by the first 1 transformations.
This is a property that we will call of of a matrix into a larger
matrix . be a matrix, with - of dimensions ( 1), and
let be the triangular factor computed by on input . Let be a matrix having
as a minor and suppose that, as a consequence of the repeated applicability of
and , the first 1 stages of algorithm on input only affects the rows that
identifies . Then contains as a minor. In other words, the factorization of
, viewed as a part of , is the same as the factorization of alone. Perhaps the
simplest example of proper embedding is =
Stage modifies the entry ( ) only if . In particular, it introduces zeros in
the -th column of without destroying the previously introduced zeros. In
view of this, and to avoid some redundant descriptions, in the rest of this paper we
will use the notation to indicate the submatrix of with elements from
rightward and downward.
Our P-completeness results are all based on reductions from the NANDCVP, a restricted
version of CVP (Circuit Value Problem) which we now briefly recall:
: the description of a -input boolean circuit composed entirely of
gates, and boolean values .
: the value ( ) computed by on input .
NANDCVP is P-complete, as reported in [10]. In order to simplify the proofs we will
further assume, without loss of generality, that each gate of has fanout at most two.
What we shall prove in this section is the following general result, which applies to any
factorization algorithm .
a b
a b
a
a
a
a
a
A
A
A
True False
3.1 Informal description
nand
duplicator
copier wire
N A XR
A
C A N
There is an encoding scheme of logical values and a log-space bounded transducer
with the following properties: given the description of a fanout 2
nand circuit and boolean inputs for , builds a matrix of
order such that, if = is the factorization computed by algorithm
, with upper triangular, then [ ] is the encoding of ( ).
We shall prove (Theorem 3.1) that the transducer exists provided that there are certain
elementary matrices with well defined properties. We will later show that such matrices
do exist for the algorithms considered in this paper.
Unfortunately, a formal description and the proof of correctness of the transducer will
require quite a large amount of details. In spite of this, the idea behind the construction is
easy, namely to repeatedly apply the proper embedding property. Hence we first describe
the reduction in an informal way and only afterward proceed to a formal derivation.
Moreover, we have actually implemented the transducer as a collection of Matlab m-
files. These and the elementary matrices for the floating point versions of Householder's
and Givens' QR algorithms (the most interesting and technical ones) are electronically
available through the authors.
Let and let and denote appropriate numerical encodings of the truth values
. We need three kinds of (square) elementary matrices for .
The first such matrix is the . It has [ apply
to compute the factorization of , we get the encoding of ( ) in the right bottom
entry of the upper triangular factor.
The second elementary matrix is the . It has [ . If we compute
an incomplete factorization of , i.e., apply all but the last transformation of to , we
get 0in the right bottom corner of the incomplete triangular factor.
The third elementary matrix is the or . It has [
compute the factorization of we get in the right bottom entry of the triangular factor.
Using these matrices as the building blocks we can construct a matrix that simulates
the circuit . The structure of is close to block diagonal, with one block for
each nand gate in the circuit . Duplicator blocks are used to simulate fanout 2 nand
gates, and wire blocks to route the computed values according to the circuit's structure.
As the factorization of a block diagonal matrix could be performed by independently factoring
the single blocks, a certain degree of overlapping between the blocks is necessary
to pass the computed values around.
To illustrate how the preceding scheme can work in practice, and to see where the
difficulties may appear, consider first the construction of a submatrix which simulates a
gate. The basic idea is to append a duplicator to a nand block as pictorially
shown in Figure 1 (left). The block is the dark gray area, the block is light gray,
and the white zones contain zeros. The right bottom entry of coincides with the top
left entry of . This is an example of proper embedding. Suppose has order . Then
after 1 stages of the encoding of ( ) is exactly where required, i.e., in the top
left entry of from where it can be duplicated. The light gray area in Figure 1 (right)
3.2 Elementary matrices
A N D W
by the time the algorithm starts
working on column
Ws
Ws
has not been modified by the transformations. The only changes in the submatrix still to
be triangularized occurred in the first row of . We already know that entry ( ) has
been modified properly, but this needs not be the case for the other entries, i.e., the black
colored ones in Figure 1 (right). For the simulation to proceed correctly, it is required
that the black entries store the rest of the first row of
. We cannot rely on their initial contents.

Figure

1: N-D matrix composition: effect of the first 1 stages
As a second example, Figure 2 pictorially describes how a block can be used to
pass a value to a possibly far away place. The block (the dark gray area in Figure 2)
is split across non consecutive rows and columns. More precisely, if is of order , the
top left dark gray area is intended to represent the principal minor of order 1 of .
As before, the white zones contain zeros, while the light gray area is of arbitrary size and
stores arbitrary values. This situation again represents a proper embedding, so that the
first 1 stages of on input (with the order of ) will result in the factorization
of the block. This implies that the (encoding of) the logical value initially in the top
left entry has been copied to the far away right bottom entry.

Figure

2: Splitting of a block
To get an idea of what a complete matrix might looks like, see Figure 3, where the
circuit computing the exclusive or of two bits is considered. The corresponding matrix
has four blocks, one block and four blocks denoted by different gray levels.
Note, however, that Figure 3 does not incorporate yet the solution to the "black entries"
problem mentioned above.
To prove the correctness of the transducer (in Theorem 3.1 below) it is convenient to
introduce a block partitioning of the elementary matrices defined in the previous section.
O
A
A
Nand matrix (N)
a b
c
c
c
I I I
O O O
I M O I O
I O
I
input output
input places
auxiliary
a b
AC
AC
R
x a a x

Figure

3: Circuit computing ( ) (left) and the structure of the corresponding
matrix (right).
Let denote one such matrix. We partition as follows
(1)
where the diagonal blocks , , and are square matrices and and also
diagonal (i.e., with only zeros outside the main diagonal). We will refer to and
as to the and submatrices and let and denote their order, respectively.
Note that an elementary matrix actually defines a set of matrices. In fact, we regard the
diagonal entries of as the . A particular elementary matrix is obtained
by filling the input place(s) with the encoding of some logical value(s).
We can now formally define the "behavior" of the elementary matrices with respect
to the factorization algorithm .
We let denote the order of , which in the
block partitioning (1), and set [ is the factorization
computed by , we have
where is the encoding of ( ). If, as often required in practice, overwrites
the input matrix, the value will replace , and this is the reason why we defined
in (1) as the output submatrix. The same remark applies to the other elementary
matrices. We also require that, for any real value , an
(0 0 ) can be defined such that Using the auxiliary
D D D
d
actual input
A
A
F
3.3 Construction and
Duplicator matrix (D)
a
a
a
Wire matrix (W)
a
a
A
D XD R
R
x z
a a x; z ; a
W XW R
R
x a a x
vectors we can solve the "black entries" problem outlined in Section 3.1. Intuitively,
by appending auxiliary vectors to the right of in (instead of simply zeros, as
in

Figure

1 (left)), we can obtain the desired values in the black colored entries of

Figure

(right). As we will see in the proof of Theorem 3.1, the initial zeros in the
auxiliary vector prevents the problem from pumping up in the construction of .
We let denote the order of , which in the
block partitioning (1), and set [ is the incomplete factorization
computed by , i.e., represents the first 2 transformations applied to by
, we obtain
We also require that, for any pair of real numbers and , an auxiliary vector
can be defined such that
We let denote the order of , which in the block
partitioning (1), and set [ is the factorization computed by
, we get
We require again that, for any real number , an auxiliary vector
(0 ) can be defined such that
Elementary matrices (including auxiliary vectors) exist for both Householder's and
Givens' methods and for Gaussian Elimination as well.
In this section we present our main result (Theorem 3.1) on the existence of a single reduction
scheme that works for any algorithm in . We still require a couple of definitions.
- Suppose that has input variables and let be the number of places (inputs
to the gates) where the variables are used. is the number of inputs counting
multiplicities. Let the gates of be sorted in topological order. Given a specific
assignment of logical values to the input variables we may then refer to the -th
as to the -th value from the input set that is required at some gate,
I
A
z -C AB @
z -C A
a
a
a
Theorem 3.1
Proof.
a b
a b
N D W
A
A n O n n C
A XR R
A
O n
C z z n
A N A A n O
A
A A A
A A g
A
has an input row at position
Let elementary matrices , , and be given for . For any fanout
nand circuit with input variables and any truth assignment to ,
we can build a square matrix such that the following holds
(a) has order , where is the number of gates in .
(b) If is the factorization computed by , then is the encoding of
(c) has a number of input rows which equals the number of actual inputs to
; the -th such row has either the structure (2) or (3) depending on whether the
-th actual input to enters the first or the second input of some gate.
(d) Any actual input affects only through one input place.
The construction can be done by using work space.
- When we say that , we intend that initially row
is either
or
where is the encoding of one of the actual inputs to .
(log
We prove the result by induction on . Let , for , be the actual
inputs to and let and be the encodings of and , respectively. The case
is easy. We only have to set = , with [
Property (b) follows from the definition of . Properties (c) and (d) are easily verified as
well. In particular, has exactly 2 input rows at positions 1 and 2, with structure (2)
and (3), respectively, and this clearly matches the number of actual inputs to . Finally,
the actual inputs affect only through the input places [
Now suppose that the number of gates in is 1, and let be a topological
ordering of the DAG representing . Clearly, all the inputs to are actual inputs to .
Let be the circuit with removed and any of its outputs replaced with , the first
input variable. Since has 1 gates, we may assume that can be constructed
which satisfies the induction hypothesis. To build we simply extend to take
into account. There are two cases, depending on the fanout of . We fully work out only
the case of fanout 1, the other being similar but just more tedious to develop (and for the
reader to follow) in details.
1. Let be the gate connected to the output of . Suppose, w.l.o.g., that provides
the first input to . By the induction hypothesis (in particular, by (c)) has the
following structure.
a
a
a
a b
a
I I
I I I
I
O O O
O O O
I I O
I I I
O O O
I I
A
y
Z z z Z Z Z
x
A
N N N O O O O O
N N N A O a O O O
N N N a ff
O O W W O W A A O
O O O X x x
W W W a a
y
O O O Z z z Z Z Z
A A w
O
A a
a ff W W
N N W a a
A
O
W W W a a
y
O Z z z Z Z Z
0has an input row at some position corresponding to the first input to gate
and is the encoding of . Note that, by property (d), the actual logical value
encoded by only affects the definition of through the entry ( ). Using
and and elementary matrices we define as follows
. The minor enclosed in boxes is a set of 1
( is the order of ) auxiliary column vectors for that we choose such
where is the matrix that factorizes . Observe that only the -th row of has
been modified by replacing ,
In what follow we regard as a block 8 10 matrix, and when we refer to the -th
row (or column) we really intend the -th block of rows (columns). Nonetheless,
is square, if is, with order plus the size of . Using part (a) of the
induction hypothesis we then see that It is easy to prove that enjoys
properties (b) through (d) as well. Assume has actual inputs. Since has
fanout 1 has 1 actual inputs and, by induction, has 1 input rows.
Now, by the above construction has exactly (
which proves (c). Property (d) also easily holds. To prove (b) we use the properties
of . By , the application of 1 stages of to only affects the first 3
(block of) rows. Thus (including its auxiliary vectors) is properly embedded in
and hence after the first 1 stages of we get
c
c
I
I I I
I I
I I
a a N N
A w
A
y
Z z z Z Z Z
A A N
A - d w
A
A n
A
A
O n N D W
O n
is the encoding of ( ). The submatrix enclosed in boxes is a
set of 2 auxiliary vectors for that we choose such that
where is the transformation matrix that triangularizes . Note that the entries
corresponding to the first elements of auxiliary vectors contain zero, as required. If
the first element (in the definition) of auxiliary vectors were not zero, we would be
faced with the additional problem of guaranteeing that the first 1 stages would
set these entries to the required values.
It is again easy to see that (including its auxiliary vectors) is properly embedded
in so that additional 1 stages of leads to
. The correctness now follows from the induction hypothesis and
property .
2. The full description of the fanout 2 case is definitely more tedious but introduces
no new difficulties. extends by means of an initial block, followed by a
block, followed by two blocks. Taking the partial overlappings into account,
it immediately follows that the order of is plus the order
of .
The construction of matrix can be done in space proportional to log by simply
reversing the steps of the above inductive process. That is, instead of constructing
and using it to build , which would require more than logarithmic work space, we
compute and immediately output the first in case of a
first gate with fanout 2) rows and columns. We also compute and output row (or the
two rows where the output of a fanout 2 gates has to be sent). All of this can be done in
space (log ) essentially by copying the elementary matrices , , and to the output
medium. The only possible problem might be the computation of , but this is not the
case. In fact, for any 1 , let ( ) be the number of fanout 2 nand gates preceding
gate in the linear ordering of . This information can be obtained, when required, by
repeatedly reading the input and only using (log ) work space for counting. It easily
follows from the above results that the index of the top left entry ( ) of the -th
block is
Hence will be either ( ) or depending on whether the value under consideration
is the first or second input to .
E2
a b ab
nn
qr
True False
c
9 9 9 102
4 Householder's QR decomposition algorithm
Theorem 4.1
Proof.
Theorem 4.2
Proof.
N a b
c
c
a
A n
H A A n n A QR
O n
O n
x
D QR R R R R R z R x
z x
is in P under both exact and floating point arithmetic.
is logspace hard for P under the real number model of arithmetic.
The Matlab program that implements the transducer is indeed log-space bounded. It
only uses the definition of the blocks and simple variables (whose contents never exceed
the size of ) in magnitude. No data structure depending on is required. Clearly, as it
is implemented using double precision IEEE 754 arithmetic, it can properly handle only
the circuits with up to approximately 2 gates.
In this section we prove that HQR is presumably inherently sequential under both exact
and floating point arithmetic. This is done by proving that a certain set , defined in
terms of HQR's behavior is logspace complete for P.
is the factorization computed by HQR, and
Note that by HQR we intend the classical Householder's algorithm presented in many
numerical analysis textbooks. In particular we refer to the one in [9]. This is also the
algorithm available as a primitive routine in scientific libraries (such as LINPACK's
[6]) and environments (like Matlab's [14])
We begin by the ready to hand result about the membership in P.
Follows from standard implementations, which perform ( ) arithmetic operations
computations (see, e.g., [9]).
According to the result of Section 3, to prove that is also logspace hard for P it
is sufficient to exhibit an encoding scheme and the elementary matrices required in the
proof of Theorem 3.1. As we will see, however, the floating point case asks for additional
care to rule out the possibility of fatal roundoff error propagations.
We simply list the three elementary matrices required by Theorem 3.1. For each
elementary matrix , the corresponding auxiliary vector is shown as an additional
column of . That the matrices enjoy the properties defined in Section 3.2 can be automatically
checked using any symbolic package, such as Mathematica .
It is the 9 10 matrix of Figure 4, where 1 1 are the encoding of logical
values (1 for and 1 for ) and is an arbitrary real number. Performing
8 steps of HQR on input
= is the arithmetization of ( ) under the selected encoding.
It is the 6 7 matrix shown in Figure 5 (left). Performing 4 steps of HQR on input
where and are arbitrary real numbers.
a
x
x z
x z
x
a
Theorem 4.3
Proof.
A
computed
is logspace hard for P under finite precision floating point arithmetic.
More precisely, to the best possible approximations of the blocks under the particular machine
arithmetic.

Figure

4: The block for HQR.

Figure

5: The and blocks for HQR.
It is the 2 3 matrix of Figure Figure 5 (right). Performing 1 step of HQR on input
for an arbitrary real number.
Applying a floating point implementation of HQR to any single block defined in Theorem
4.2 results in approximate results. For instance, we performed the decomposition
of the four matrices using the built-in function available in Matlab. We found that
the relative error affecting the computed encoding of ( ) ranged from a minimum
of 0 5 to a maximum of 3 . Here is the roundoff unit and equals
IEEE 754 standard arithmetic. These might appear insignificant errors. However, for a
matrix containing an arbitrary number of blocks, the roundoff error may accumulate to
a point where it is impossible to recover the exact (i.e., 1 or 1) result. Clearly, direct
error analysis is not feasible here, since it should apply to an infinite number of reduction
matrices. Our solution is to control the error growth by "correcting" the intermediate
results as soon as they are "computed" by nand blocks. Note that, by referring to the
values by a certain elementary matrix , we properly intend the non zero values
one finds in the last row of the triangular factor computed by HQR on input (including
the auxiliary vectors). Analogously, the input values to are the ones computed by the
elementary matrix preceding in .
We take duplicator and wire blocks as in Theorem 4.2, and provide a new definition
for nand blocks so that they always compute exact results. To do this, we have to
consider again the structure of , as resulting from Theorem 3.1.
and
a
a
a a
corr
corr
corr
corr corr
z -
z - z -B B @C C A
l C l
(1) (2) ( 1)
c
a x N N
x
N a
a N
x
Let be the -th gate in the topological ordering of , and let and be the gates
providing the inputs to . Let denote the block of corresponding to ,
according to the construction of Theorem 3.1. To prove the result we maintain the
invariant that the values computed by are exact. This is clearly
true 1. Using the invariant we first verify that the errors affecting the values
computed by can be bounded by a small multiple of the roundoff unit. We then use
the bound to show how to redefine so that it computes exact results, thus extending
the invariant to .
From the proof of Theorem 3.1 we know that the output of (and similarly of ) is
placed in one of the input rows of as a consequence of the factorization of possibly a
followed by a block. It follows that the error affecting the output of is only due
to the above factorizations to the factorization of itself. Since there is a limited
number of structural cases (depending on the fanout of gates and ) and considering
all the possible combinations of logical values involved, the largest error ever affecting the
output of can be determined by direct (but tedious) error analysis or more simply
by test runs. For the purpose of the following discussion we may safely assume that the
relative errors affecting the computed quantities are bounded by , for some constant of
order unit ( is actually smaller than 10). In other words, we may assume that the actual
outputs of are Recall that is the last entry
of the generic auxiliary vector ( ) of after the factorization (see the definition of
is Section 3.2). Here, however, we require that be a machine number (i.e., a rational
number representable without error under the arithmetic under considerations).
Having a bound on the error, we are now ready to show how to "correct" the erroneous
outputs. The new nand block, denoted by , extends with two additional rows and
columns, as shown below
positive integer (to be specified below). Note
that ( 1) is precisely that auxiliary vector for the old that produces 1 as output,
. The auxiliary vector for is (0
first requirement on is thus that the quantity be a computer number. As
the length of the significand of we see that a sufficient condition is that
the length of the significand of does not exceed 1. Now, let us apply HQR
to extended by its auxiliary vector. As is properly embedded in , after 8
stages of HQR we get (using the above result on the error)
A second condition on is that we want that 2 to get rid of the
error . An easy argument shows that this implies log . Thus, recalling the bound
planerot
nn
corr
corr
corr
F
a a
a
a
Theorem 5.1
Proof.
Theorem 5.2
I ;
x
x
x
G A A n n A QR
G
G
5 QR decomposition through Givens' rotations
is in P under both exact and floating point arithmetic.
is logspace hard for P under the real number model of arithmetic.
on and , we see that 5 is sufficient. As a consequence, the length of cannot
exceed 6. The actual reflection matrix applied to is then2
which, by easy floating point computation, gives
Applying one more stage now leads to the correct results and . The above requirement
on is by no means a problem. In fact, the auxiliary values ever required are the non zero
elements in the input rows of the blocks that possibly follow nand elementary matrices,
i.e., and blocks. These are simply 1, 2, and 5 4, all of which can be represented
exactly with a 3 bit significand.
The elementary matrices of Theorem 4.3 are available for the general transducer implemented
in Matlab. In particular, is defined
In this section we prove that the following set
is the factorization computed by GQR, and
is logspace complete for P. The way we present the results of this section closely follows
the metodology of Section 4. Here, however, we have to spend some more words about the
particular algorithm considered. In fact, the computation of the QR decomposition can
be done in various ways using plane (or Givens') rotations. Differently from Householder's
reflections, a single plane rotation annihilates only one element of the matrix to which it
is applied, and different sequences of annihilations result in different algorithms. By the
way, this degree of freedom has been exploited to obtain the currently faster (among the
known accurate ones) parallel linear system solvers [19, 15]. We also outline that there
is no QR algorithm available in Matlab (nor in libraries as LINPACK or LAPACK), but
it just provides the primitive that computes a plane rotations. The hardness
results of this section apply to the particular algorithm that annihilates the subdiagonal
elements of the input matrix by proceding downward and rightward. This choice places
GQR in the class defined in Section 2, with the position that one stage of the algorithm
is the sequence of plane rotations that introduce zeros in one column.
See, e.g., [9]. We only point out that the membership in P holds independently
of the annihilation order.
a
a
a
z
G
True False
Proof.
a b
Theorem 5.377777777
ff
x z
x z
x z
x z
ppp
is logspace hard for P under finite precision floating point
As in Theorem 4.2, we simply list the three elementary matrices extended with
the generic auxiliary vector. The matrices are shown in Figures 6 through 8, where
are encodings of logical values (1 for and 1 for ) and and
are arbitrary real numbers. Again, that the matrices enjoy the properties defined in
Section 3.2 can be verified with the help of a symbolic package.

Figure

The block for GQR.

Figure

7: The block for GQR.

Figure

8: The block for GQR.
We now switch to the more delicate case of finite precision arithmetic.
m22
True
e
Proof.
a
a
a
a
0d e
0d e
R
d; d R
e
a
x y
We apply the same ideas of Theorem 4.3. That is, we extend the definition of
so that it always computes the exact results. Here, however, we cannot reuse the
block adopted for the exact arithmetic case. There is a subtle problem that urges for a
different definition of . Let us see in details. If we apply a floating point implementation
of GQR to we clearly get approximate results (note that the matrices for GQR contain
irrational numbers). In particular, instead of 0
0 , in the bottom right corner of
we get (1
. Even if , , , and are of the order of the roundoff
unit , the fact that is not zero causes the whole construction to fail. Note that the
same kind of approximate results are obtained under HQR, but with no damage there. To
get to the point, suppose that is in column of and let us proceed by considering
stage of the algorithm. In HQR one single transformation annihilates the whole column
so that the contribution of a tiny to the -th transformation matrix is negligible. On
the other hand, in GQR the elements are annihilated selectively and, since is not zero,
one additional plane rotation is required between rows and + 1 to place zero in the
entry ( +1 ). Unfortunately this has the effect of making the element in position ( )
positive, which is a serious trouble since this entry contained the encoding, with a small
perturbation, of a truth value. The result is that, when the subsequent plane rotations
(the ones simulating the routing of the logical value) are applied, the value passed around
is always the encoding of , and the simulation fails in general.
We thus need to replace the duplicator with one that returns a true zero in the entry
( 1) of the incomplete factor , which will clearly exploit the properties of floating
point arithmetic.
noindent Let and denote the length of the significand and the largest exponent
such that 2 can be represented in the arithmetic under consideration, respectively.
For the standard IEEE 1023. The nonzero elements of the new
duplicator are only powers of 2. In this way any operation is either exact or is simply a
no operation.
As the new auxiliary vector we define
Nothe that the only possible assignments to and are 0 and 1 or 1 and 0.
The rest of the proof is now similar to that of Theorem 4.3. We show how to correct the
slightly erroneous values computed by an block assuming that the previous blocks
GEP
GEM
a
a
a
a
a
a
a a
a
a
corr
corr
corr
corr
corr
z -
z -B @C A
@A
6 Gaussian Elimination with Pivoting8
0d e
0d e 0d e
0d e 0d e
d e d e
a
a x a x N
G
x
x
return exact results. Let stand for the nand block adopted for the exact arithmetic
version of GQR (Figure 6). The new nand block is then
As the new auxiliary vector we take
i.e., the first 10 entries of ( ) coincide with ( 2 ). Now, let us apply GQR to
extended by its auxiliary vector. As is properly embedded in , after 9 stages of
GQR we get
where , for some small constant of order unit. The plane rotation to annihilate
the entry (2 1) of is represented by
and its application in floating point gives
The crucial point is that, if can be represented
with no more than 2 significant bits, the alignment of the fraction part performed
during the execution of 2 simply cause the contribute 2 to
be lost. Hence, the computed element in the entry (1 3) will be 2 2 . But then
one more rotation produces the exact values and in the last row. Note that the only
value required in place of is 1.
In this section we consider the algorithm of Gaussian Elimination with partial pivoting,
or simply , a technique that avoids nondegeneracies and ensures (almost always)
numerical accuracy. We also consider the less-known Minimal pivoting technique, ,
one that only guarantees that a PLU factorization is found. Minimal pivoting has been
adopted for systolic-like implementations of Gaussian Elimination [11]. A brief description
of these algorithms is reported in Appendix A.
We prove that GEM is inherently sequential, unless applied to strongly nonsingular
matrices while GEP is inherently sequential even when restricted to strongly nonsingular
matrices.
Theorem 6.1
The set is log-space complete for
6.1 Partial Pivoting
6.2 Minimal Pivoting
The proof we give here builds on the original proof in [20], and hence does not share
the common structure of the other reductions in this paper. Essentially we show that,
with little additional effort with respect to Vavasis' proof, we can exhibit a reduction
in which the matrix obtained is strongly nonsingular. As already pointed out, strongly
nonsingular matrices are of remarkable importance in practical applications. This class
contains symmetric positive definite (SPD) and diagonally dominant matrices, which often
arise from the discretization of differential problems. Observe that, on input such matrices,
plain GE (no pivoting) is nondegenerate, but it is not stable in general and hence is not
the algorithm of choice.
As in [20], that GEP is inherently sequential follws from the proof that the folllowing
set is P-complete.
strongly nonsingular and, on input , GEP uses row
to eliminate column .
We postpone the technical proof of Theorem 6.1 to the Appendix C, but give an
example that shows the way the matrix given in [20] is modified. Figure 9 depicts the
reduction matrix which would be obtained according to the rules in [20] on input the
description of the circuit of Figure 3. The matrix is nonsingular; however, it can be seen
that the leading principal minor of order 2 is singular. The matrix we obtain, according
to Theorem 6.1, is shown in Figure 10. It can be easily seen that our matrix is strongly
diagonally dominant by rows, and hence strongly nonsingular.

Figure

9: Matrix corresponding to the exclusive or circuit. The symbol stands for
a zero entry.
The technique of minimal pivoting, i.e., selecting as the pivot row at stage the first
one with a nonzero entry (below or on the main diagonal) in column , is probably
the simplest modification that allows GE to cope with degenerate cases. However, such
False True
A
Theorem 6.2
Proof.
a
b a
x
x
The set is logspace complete for under both real and finite precision
floating point arithmetic.

Figure

10: The matrix for the computation of ( ).
a simple technique is sufficient to make the Gaussian Elimination algorithm inherently
sequential. Note that, even if no formal error analysis is available for GEM, it is not
difficult to exhibit matrices (that can plausibly appear in real applications) such that the
error incurred by GEM is very large. Actually, GEM is likely to be as unstable as plain
GE.
Consider the following set.
is the factorization computed by GEM, and
The set is clearly in , as GEM runs in time ( ) under both models of arith-
metic. We first show that is also logspace hard for when the input matrices are
singular, and then show how to restrict the input set. As GEM belongs to the class , to
prove the hardness of we simply list the three elementary matrices required by Theorem
3.1. Note that the matrices are the same for both models of arithmetic, as the operations
performed by GEM in floating point are exact. The encoding of logical values here is 0
for and 1 for . The matrices are depicted in Figures 11 and 12.

Figure

11: The nand (left) and wire (right) blocks for GEM.
x
z
z
U
A
A
k A A

Figure

12: The block for GEM.
is clearly singular. Now consider the following matrix of order 2 , where
is the order of .
where is the matrix with 1 on the antidiagonal and 0 elsewhere. The determinant of
can be easily proved to be 1. Moreover, if = is the factorization computed
by GEM, then . Note then that what we prove to be P-complete is not exactly
, but a set with a little more complicate definition (which is still in terms of GEM's
behavior).
As usual, in the following the notation will be used to denote the matrix obtained
after 1 stages of GEM on input , and considering only the entries ( ) with .
However, by writing we intend the submatrix obtained after 1 stages of GEM
(on input ), and considering only the entries ( ) such that . With this
position, to prove that the output of can be read off entry ( ) of the factor of
, we show that the executions of GEM on input and on input result in identical
submatrices and , for 0 . The proof is by induction. Initially,
the equality follows from the definition of . Consider stage 1. If column of
contains a nonzero element below or on the main diagonal, say at row index , then the
selected pivot row is the -th under both executions. The results follows then from the
induction hypothesis and the fact that exactly the same operations are performed on the
elements of the submatrices. If no nonzero element is found in column of , then
stage of the first execution has no effect, and hence = . Under the second
execution (the one on input ) by construction the pivot is taken from row 2 + 1.
However, the pivot is the only nonzero element in row thus the effect
of this step is simply the exchange of rows and 2 + 1. But then once more
.
We conclude by observing that the set of Theorem 6.2 is clearly NC computable
when the input set is the class of strongly nonsingular matrices. In fact, in this case,
GEM and plain Gaussian Elimination behave exactly the same.
Theorem 7.1
7 On algorithms for the
l l
Computing the PLU factorization returned by GEMS on input a nonsingular
matrix is in arithmetic .
A n A
O
R
A
In this section we show that a known algorithm for computing a PLU decomposition
of a nonsingular matrix (see [7]) corresponds to GE with a nonstandard pivoting strategy
which is only a minor variation of Minimal pivoting. This result seems to be just a
"curiosity"; however, we can prove that the same strategy is inherently sequential on input
arbitrary matrices, which can be seen as a further evidence of the difficulties of finding
an algorithm to compute the PLU decomposition of possibly singular matrices.
The new strategy will be referred to as Minimal pivoting with circular Shift, and
the corresponding elimination algorithm simply as GEMS. The reason for its name is
that GEMS, like GEM, searches the current column (say, column ) for the first nonzero
element. Once one is found, say in row , a circular shift of row through is performed
to bring row in place of row (and the latter in place of row 1).
We consider the algorithm of Eberly [7]. Given , nonsingular of order , let
denote the matrix formed from the first columns of
the set of indices of the lexicographically first maximal independent subset of the rows of
since has full column rank. Moreover,
Note that the computation of all the is in (see [3]). Now,
. Then a permutation such that has
factorization is simply
where is the th unit (column) vector. Clearly, once has been determined, computing
the factorization of can be done in polylogarithmic parallel time using known
algorithms. We now show by induction on the column index that is the same permutation
determined by GEMS. The basis is trivial, since is the index of the first nonzero
element in column 1 of . Now, for 1, let
be the (partial) factorization computed by GEMS, where is upper triangular with
nonzero diagonal elements (since is nonsingular) and the unit vectors
extend to form a permutation matrix. Clearly, Minimal Pivoting ensures that
. Now, the next pivot row selected by GEMS is the one corresponding to
the first nonzero element in the first column of . Let the index of
the pivot row. Since Gaussian Elimination does nothing but linear combinations between
rows, it follows that the initial matrix satisfies
This in turn implies that = , i.e. that = .
nn
F
Theorem 7.2
Proof.
8 Conclusions and Open Problems
a x; z
a
is logspace hard for P.
except for the auxiliary vector od D
We now show that GEMS is inherently sequential by proving that the set
is the factorization computed by GEMS, and
is P-complete. Clearly, that is in P is obvious, so what remains to prove is the
following.
Once more GEMS is in the class . So we simply give the elementary matrices.
This is very easy. Everything is the same as in the first part of the proof of Theorem 6.2,
. The new definition for ( ) is
It is an easy but interesting exercise to understand why the second part of Theorem
6.2, which extend the P-completeness result to nonsingular matrices, here does not work
(we know that it cannot work, in view of Theorem 7.1).
The matrices corresponding, for both Householder's and Givens' algorithms, to a circuit
are singular, in general. More precisely, the duplicator elementary matrix is singular,
so that all the matrices that do not correspond to simple formulas (fanout 1 circuits)
are bound to be singular. All the attempts we made to extend the proofs to nonsingular
matrices failed. The deep reasons of this state of affairs could be an interesting
subject per se. To see that the reasons for these failures might be deeper than simply
our technical inability, we mention a result of Allender et al. [1] about the "power" of
singular matrices. They prove that the set of singular integer matrices is complete for
the complexity class C L . The result extends to the problem of verifying the rank of
integer matrices. Of course, our work is at a different level: we are essentially dealing
with presumably inherently sequential algorithms for problems that parallelize very well
(using different approaches). However, the coincidence suggests that nonsingular matrices
might not have enough power to map a general circuit. This is the major open problem
for the QR algorithms.
Also, for general matrices, it would be interesting to know the status of Householder's
algorithm with column pivoting, which is particularly suitable for the accurate rank determination
under floating point arithmetic.
For what concerns Givens' rotations, an obvious open problem is to determine the
status of other annihilation orderings, especially the ones that proved to be very effective
in limited parallelism environments [19, 15]. We suspect that these lead to inherently
sequential algorithms as well.
A M



--R

Inherently Seq.
NC5 A set is in C L provided that there is a nondeterministic logspace bounded Turing machine such that iff has the same number of accepting and rejecting computations on input

GEM NC GEMS NC Table 1: Parallel complexity of GE with different pivoting strategies and for different classes of input matrices.

As already mentioned
The complexity of matrix rank and feasible systems of linear equations


Fast parallel matrix inversion algorithms


Efficient Parallel Independent Subsets and Matrix Factorizations
Parallel Linear Algebra


Matrix Computations Limits to Parallel Computation Introduction to Parallel Algorithms and Architectures: Arrays Trees Hypercubes Journal of Computer and System Sciences
Theoretical Computer Science Computational Geometry
The Algebraic Eigenvalue Problem


On the Parallel Complexity of
Parallel complexity of Householder QR factorization

An alternative Givens ordering
Complexity of Parallel Matrix Computations


On Stable Parallel Linear System Solvers


GE computes the decomposition of (whenever it exists) by determining a sequence of 1 elementary transformations with the following properties:
pivoting strategy pivot row
GE with Minimal pivoting (GEM).
QR factorization via Householder's reflections (HQR).

GEP computes a decomposition of

HQR applies a sequence of 1 elementary orthogonal transformations to

a a a a a a
GQR applies to general real matri- ces
The rotation used to annihilate a selected entry of a matrix is the orthogonal matrix defined as follows:
A floating point system is characterized

In particular

Let denote the matrix corresponding to the circuit according to Vavasis' proof
Vavasis' proof is based on the observation that a NAND gate outputs unless one of its inputs is
Our matrix has order 3
Define the of the matrix to be the set of odd-numbered columns of the main submatrix
The proof of the theorem is now a consequence of the following two lemmas.

For the following facts hold of GEP on input
Figure 13: The structure of the matrix

Consider step 2 1 of GEP.





For the entries with row index larger than 2

By induction hypothesis we know that the first 2 2 elimination steps did not affect the auxiliary columns 2 2
In order to prove (d)


Suppose now that the pivot at step 2 1 is 4 0.

--TR
numerical stability;NC algorithms;inherently sequential algorithms;p-complete problems;parallel complexity;matrix factorization