A Parallel Algorithm for the Reduction to Tridiagonal Form for Eigendecomposition.
One-sided orthogonal transformations which orthogonalize columns of a matrix are related to methods for finding singular values. They have the advantages of lending themselves to parallel and vector implementations and simplifying access to the data by not requiring access to both rows and columns. They can be used to find eigenvalues when the matrix is given in factored form. Here, a finite sequence of transformations leading to a partial orthogonalization of the columns is described. This permits a tridiagonal matrix whose eigenvalues are the squared singular values to be derived. The implementation on the Fujitsu VPP series is discussed and some timing results are presented.
Introduction
Symmetric eigenvalue problems appear in many applications ranging from computational
chemistry to structural engineering. Algorithms for symmetric eigenvalue
problems have been extensively discussed in the literature [11, 9] and implemented
in various software packages (e.g. LAPACK [1]). With the broader introduction of
parallel computers in scientific computing new parallel algorithms have been suggested
[7, 2]. In the following another new parallel algorithm is suggested which is
particularly well adapted to vector parallel computers and has low operation counts.
Eigenvalue problems can only be solved by iterative algorithms in general as they
are in an algebraic sense equivalent to finding the n zeros of a polynomial. There are,
however, two main classes of methods to solve the symmetric eigenvalue problem.
The first class only requires matrix vector products and does not inspect nor alter
the matrix elements of the matrix. This class includes the Lanczos method [9] and
Date: November 1995.
1991 Mathematics Subject Classification. 65Y05,65F30.
Key words and phrases. Parallel Computing, Reduction Algorithms, One-sided Reductions.
Computer Sciences Laboratory, Australian National University, Canberra ACT 0200, Australia.
ANU Supercomputer Facility, Australian National University.
Centre for Mathematics and its Applications, Australian National University.
M. HEGLAND   , M. H. KAHN    , M. R. OSBORNE
has particular advantages for sparse matrices. However, in general, the Lanczos
method has difficulties in finding all the eigenvalues and eigenvectors.
A second class of methods iteratively applies similarity transforms
with A A to the matrix to get a sequence of orthogonally similar matrices which
converge to a diagonal matrix. This second class of methods consists mainly of two
subclasses. The first subclass uses Givens matrices for the similarity transforms and
is Jacobi's method. It has been successfully implemented in parallel [2], [12]. A disadvantage
of this method is its high operation count. A second subclass of methods
first reduces the matrix with an orthogonal similarity transform to tridiagonal form
and then uses special methods for symmetric tridiagonal matrices. Both parts of
these algorithms pose major challenges to parallel implementation. For the second
stage of tridiagonal eigenvalue problem solvers the most popular methods include
divide and conquer [7] and multisectioning [9]. Here the reduction to tridiagonal
form is discussed. Earlier algorithms use block matrix algorithms, see [6, 5, 3].
However, these methods have not achieved optimal performance. One problem is
that similarity transforms require multiplications from both sides.
It was seen [12] that the Jacobi method based on one-sided transformations
allows better vectorization and requires less communication than the original Jacobi
algorithms. Assuming that A is positive semi-definite the intermediate matrices B i
can be defined as factors of the A i , that is,
As will be seen in the following, the one-sided idea can also be used for the reduction
algorithm.
The algorithm will form part of the subroutine library for a distributed memory
computer, the Fujitsu VPP 500. Often the application of subroutines from libraries
allow the user little freedom in the choice of the distribution of the data to the
local memories of the processors. The one-sided algorithms allow a large range of
distributions and perform equally well on all of them.
In the next section the one and two-sided reduction to tridiagonal form is de-
scribed. Section 3 reinterprets the reduction as an orthogonalisation procedure similar
to the Gram-Schmidt procedure. This reinterpretation is used to introduce the
new one-sided reduction algorithm. In Section 4 the computation of eigenvectors is
discussed and Section 5 contains timings and comparisons with other algorithms.
REDUCTION TO TRIDIAGONAL FORM 3
2. Reduction to tridiagonal form
A class of methods to solve the eigenvalue problem for symmetric matrices A 2
R n\Thetan first reduces them to tridiagonal form and in a second step solves the eigenvalue
problem for this tridiagonal matrix. The problem of finding the eigendecomposition
of the tridiagonal matrix will not be discussed here but that of accumulating
transformations to be used for finding the eigenvectors of the symmetric matrix is
The reduction to tridiagonal form produces a factorization
where Q is orthogonal and T is symmetric and tridiagonal. If the offdiagonal elements
of T are (nonzero) positive the factorization is uniquely determined by the
first column of Q [9]. The proof of this fact leads directly to the Lanczos algo-
rithm. While the Lanczos method has advantages especially for sparse matrices,
methods based on sequences of simple orthogonal similarity transforms [9, p.118]
are preferable for dense matrices.
2.1. Householder's reduction. A method attributed to Householder uses Householder
transformations or reflections
with In the following let ff ij denote the elements of A. That is,
matrix H(w)AH(w) has zeros in rows 3 to n in the first column and in columns 3 to
n in the first row. The computation of H(w), or equivalently of w; fl and fi 1 requires
multiplications and n additions (up to O(1)).
Using the matrix vector product v
the application of H(w) is a rank two update as
4 M. HEGLAND   , M. H. KAHN    , M. R. OSBORNE
where
This takes n 2 +n multiplications and additions for the computation of p, 2n+2
multiplications and additions for the computation of q and n 2 multiplications
and 2n 2 additions for the rank two update. This gives a total of 2n 2
multiplications and 3n 2
In a second step, a v is found such that H(v)H(w)AH(w)H(v) has additional
zeros in columns 4 to n in the second row and in rows 4 to n in the second column
and the procedure is repeated until the remaining matrix is tridiagonal. The sizes
of the remaining submatrices decrease and at step n \Gamma k a matrix of size k has to be
processed requiring multiplications and 3k 2 additions. This
gives a total of
multiplications and
additions. The tridiagonal matrix is not uniquely determined by the problem. For
example, different starting vectors for the Lanczos procedure lead to different tridiagonal
matrices. Also, different matrices can be obtained if different arithmetic
precision is used [9, p.123/124]. However, despite this apparent lack of definition,
the eigenvalues and eigenvectors of the original problem can still be determined with
an error proportional to machine precision.
In summary, the sequential Householder tridiagonalization algorithm is as follows:
For
Calculate (fl; w) from A(k
End
For vector and parallel processors the Householder algorithm has some disad-
vantages. Firstly, as the iterations progress, the length of the vectors used in the
calculations decreases but for efficient use of a vector processor we prefer long vector
lengths. Secondly, in a parallel environment, if the input matrix A is partitioned in
a banded manner across a one-dimensional array of processors, the algorithm will
REDUCTION TO TRIDIAGONAL FORM 5
be severely load imbalanced. To avoid this various authors have suggested using
cyclic [5] or torus-wrap mappings of the data [10, 3]. Also, for a parallel imple-
mentation, the rank two update of A requires copies of both vectors w and q on all
processors which leads to a heavy communication load.
2.2. One-sided reduction. A one-sided algorithm is developed to overcome the
difficulties inherent in a parallel version of the sequential Householder algorithm. In
addition, it is expected that the one-sided algorithm generates less fill-in for sparse
matrices than the two-sided algorithms if the matrix is given in finite element form.
Real symmetric matrices A are either given by or can be factored as
so A can be represented in factored form by D and B. If Cholesky factorization
is to be used then the spectrum of A might have to be shifted so that A \Gamma -I is
positive definite. The parameter - can be chosen using the Gershgorin shift. If
exact arithmetic is used for the reduction, the eigenvectors of A \Gamma -I are equal to
the eigenvectors of A and the eigenvalues are shifted by -. However, the precision
of the computations will be affected by the introduction of the shift.
A Householder similarity transform of A is done by applying H(w) to B as
For this, a rank one modification must be computed
with flBw. The computation of p requires multiplications and
additions, the rank one modification n 2 multiplications and n 2 additions giving
together 2n 2 +n multiplications and 2n additions. In contrast to the two-sided
algorithm, the number of additions is approximately the same as the number of
multiplications in this algorithm. This is advantageous for architectures which can
do addition and multiplication in parallel as it means better load balancing.
The computation of w is more costly for this method than for the original Householder
method. As the Householder vector w is computed from the first column a 1
of A, this column has to be reconstructed first from the factored representation by
This requires n 2 +n multiplications and additions. Thus the computation of
the Householder matrix H(v) requires (up to O(1) terms) multiplications
and n 2 additions.
6 M. HEGLAND   , M. H. KAHN    , M. R. OSBORNE
Adding these terms up, one reduction step needs multiplications and
additions and so the overall costs of this algorithm are
multiplications and
additions. The total number of operations has increased compared with the original
algorithm. But the time used on a computer which does additions and multiplications
in parallel and at the same speed is the same. If the matrix A is not already
factorized, however, the time to do this would have to be taken into account as well.
In summary, the sequential version of the one-sided algorithm is as follows.
For
Calculate (fl; w) from a k
End
3. The one-sided algorithm as an orthogonalization procedure
In order to develop the parallel version of the algorithm the one-sided reduction is
interpreted as an orthogonalization procedure like the Gram-Schmidt process. Let
b i denote the ith column of B, that is,
Then the matrix A can be interpreted as the Gramian of the b i as follows,
The one-sided tridiagonalization procedure constructs
such that As T
is the Gramian of the c i the tridiagonality is a condition on the orthogonality of
certain c i as
2: (3.4)
REDUCTION TO TRIDIAGONAL FORM 7
In a first step the orthogonality of c established by setting
and
such that
and
Here ffi kj denotes the Kronecker delta. In a subsequent step, linear combinations of
are formed such that c 4 are orthogonal to c 2 . As the c are
orthogonal to c 1 the linear combinations are as well and so the subsequent steps do
not destroy the earlier orthogonality relations. This is the basic observation used in
the proof of this method.
The algorithm for reduction to tridiagonal form is then as follows:
c (1)
for k :=
end for
where the fl (k)
ij are such that the new c (k)
are orthogonal to c (k\Gamma1)
. That is,
and the matrix
i;j=k;:::;n
is orthogonal with
That is, at each iteration k, find an orthogonal transformation of the c (k\Gamma1)
i such
that c (k)
j is orthogonal to c (k\Gamma1)
This is equivalent to making
the offdiagonal elements ff j;k\Gamma1 and ff k\Gamma1;j of A zero for
Proposition 3.1. Let c (k)
j be computed by the previous algorithm and
c (1)
8 M. HEGLAND   , M. H. KAHN    , M. R. OSBORNE
where c (n)
tridiagonal and there is an orthogonal
matrix Q such that
Proof. The proof is based on the fact mentioned earlier that some orthogonality
relations are invariant. It uses induction. The proof is very similar to the one given
in the next section for the corresponding parallel algorithm.
Remark. The original Householder algorithm can be formulated in a similar way
treating B as an inner product. Coordinate transformations change the matrix until
it is tridiagonal.
3.1. Parallel Algorithm. In the following the single program multiple data model
(SPMD) will be used. The basic assumption is that all the processors are programmed
in the same way although their actions might be slightly different. Thus
an SPMD algorithm is described by the pseudocode denoting what one processor
has to do. The data in the matrix B is distributed to the processors by columns in a
cyclic fashion. This means that processor 1 contains columns
contains and so on where p is the total number of processors. More
formally, processor q contains b i for
pg. In order
to simplify notation let N ng where is the processor
number. Furthermore mod denotes the mod function mapping positive and negative
integers to
c (1)
broadcast c (1)
1 else receive c (1)for k :=
gather c (k)
end for
As before, the coefficients fl (h)
ij and ~
ij are such that the modified columns are
orthogonal to the last unmodified one. Thus
ng " N q ;
and
REDUCTION TO TRIDIAGONAL FORM 9
The coefficients also form orthogonal matrices so it follows that,
and
~
The c (k)
are overwritten with the ~ c (k)
. Note that the last calculations of the c (k)
are
duplicated on all p processors which leaves c (k)
k on all processors ready for the next
step. The only communication required for each iteration is in the gathering of at
most p vectors c (k)
.
Proposition 3.2. Let c (k)
j be computed by the previous algorithm and
c (1)
where c (n)
tridiagonal and there is an orthogonal
matrix Q such that
Proof. Let
c (1)
First show that, for
1. there is an orthogonal matrix Q (k) such that
2. the linear hull of c (k)
n is orthogonal to the linear hull of the c (1)
and
3. T
Proposition 3.2 is a consequence of this, obtained by setting 1. The proof
uses induction over k. The statement is easily seen to be true in the case of
with I, the n dimensional identity matrix. In this case the orthogonality
conditions are empty and the matrix T (1) is a two by two matrix and thus tridiagonal.
The remainder of the proof consists of proving the induction step. Assume that
the three properties are valid for m. Then it has to be shown that it is
also valid for
?From the construction it follows that C for an orthogonal G (m) .
Thus, with Q one retrieves the first property. In particular, the
existence of G (m) follows from the existence of the constructed Householder matrices.
Then the linear hull of c (m+1)
n is a subspace of the linear hull of
c (m)
n and thus orthogonal to c (1)
. Furthermore they have been
constructed to be orthogonal to c (m)
. From this the second property follows.
M. HEGLAND   , M. H. KAHN    , M. R. OSBORNE
Finally, T which is tridiagonal
and it remains to show that the the m+ 2nd column of T
has zeros in the first m rows. But the first m elements of this column are just
c (m+1)
m. They are zero because of the second property.
This proof can be used for both the sequential and the parallel algorithm.
In practical implementations the orthogonalisation of the c (k)
j is achieved by forming
the products c (k\Gamma1)
which correspond to individual elements in the up-dated
version of the symmetric matrix A. Householder transformations are then
used to zero the relevant off-diagonal elements. At the first stage in each iteration
these transformations are carried out locally and after the gathering step the transformation
on the (at most) p formed elements of A is replicated on all processors.
Although Householder transformations are used here, it would also be possible to
use Givens transformations.
3.2. analysis. When using the one-sided reduction to tridiagonal as one step
in a complete eigensolver there are several possible components to the error. First
there is the Cholesky factorization. However a result in Wilkinson [11] (equation
4.44.3) shows that the quantity extremely small. Thus the error
incurred in working with the Cholesky factor B in the subsequent calculations is
minimal.
In the second step, the tridiagonal matrix C (n\Gamma2) is produced by successive reduction
of orthogonal similarity transforms where C (k) is the transformed
matrix calculated after the k'th stage. The error in the eigenvalues of C (n\Gamma2) is
bounded by the numerical error in the transform [11]. After the k'th step this is
given by
Y
where P i is the exact orthogonal matrix corresponding to the actual computed data
at stage i. A bound on this difference follows from Wilkinson [11] (equation (3.45.3)).
In fact,
Y
ij and t is the word length. So the error introduced by the
reduction to tridiagonal is small.
The final stage is the calculation of the eigenvalues of the tridiagonal matrix.
These are determined to an accuracy which is high relative to the largest element
of the tridiagonal matrix. This applies for example to eigenvalues computed using
REDUCTION TO TRIDIAGONAL FORM 11
the Sturm procedure. But this result does not guarantee high relative accuracy in
the determination of small eigenvalues so, if this is important, the Jacobi method
becomes the method of choice [4].
4. Calculating Eigenvectors
In order to calculate the eigenvectors of the symmetric matrix, the orthogonal
matrix Q defining the reduction is accumulated. This is achieved by starting with
the identity matrix (distributed cyclically over the processors), then updating it by
multiplying it by the same Householder transformations used to update C. Forming
Q explicitly is preferable to storing the details of the transformations and then
applying them to the eigenvectors of the tridiagonal matrix which is the usual procedure
for sequential implementations. The reason for this is that the matrix of Q
is distributed in a way which renders multiplication from the left inefficient. This is
discussed more fully in the following.
4.1. Calculation of Eigenvectors for One-Processor Version. The n \Theta n symmetric
matrix A is reduced to tridiagonal form by a sequence of Householder transformations
represented by the orthogonal n \Theta n matrix Q and
n\Lambdan is tridiagonal. For the one-sided algorithm, we have actually calculated
The eigendecomposition of T is given by
where   is a diagonal matrix containing the eigenvalues and V 2 R n\Lambdan is the matrix
of eigenvectors. Combining the above two equations gives
and so the eigenvectors of A are the columns of
The matrix Q obtained during the rediagonalization procedure is represented by
the Householder matrices H
The eigenvector matrix V of the tridiagonal matrix T is represented by its matrix
elements. So in order to get the matrix element representation of U form the product
M. HEGLAND   , M. H. KAHN    , M. R. OSBORNE
This product can be done in two different ways. The first method multiplies V
with the H k . Formally, define a sequence U such that
U
This method is often used, for example [8],
where it is called backward transformation.
A second method computes the element-wise representation of Q first by the
recursion
specifically computes the matrix-vector product
QV where . This is called forward accumulation in [8].
The difference between these two methods is that the first one applies the Householder
transforms H k from the left and the second applies them from the right. In
addition the second method requires the computation of a product of two matrices
in element form. For the sequential case when using the two-sided Householder re-
duction, the first method is preferred as it avoids matrix multiplication. However,
for the multi-processor version, multiplication from the left by the Householder
transformation requires extra communication when the columns of the matrix of
eigenvectors V are distributed over the processors. In fact, one purpose of the one-sided
reduction is to avoid this communication in the reduction to tridiagonal form.
There is certainly communication required for the matrix multiplication QV but
this is of fewer large blocks of the matrix so will be less demanding.
4.2. Multi-processor Version. In the multi-processor case there is an added complication
in that B is assumed to be cyclically distributed (although it is not explicitly
laid out as such). So use
~
where P is the permutation that transforms B to cyclic layout.
The matrix Q is a product of 2 \Theta (n \Gamma 2) matrices formed from the Householder
transformations.
Here, H (1)
j refers to the transformations carried out locally at each step of the reduction
and H (2)
j refers to the transformation using one column of B from each processor
REDUCTION TO TRIDIAGONAL FORM 13
which is carried out redundantly on all processors. To find Q start with an identity
matrix distributed cyclically across the processors to match the implicit layout
of B. This matrix is then updated by transforming it with the same Householder
transformations that are used to update B to obtain the tridiagonal matrix.
The matrix of eigenvectors V of the tridiagonal matrix is obtained in block layout
whilst the matrix Q is in cyclic layout. To find the eigenvectors of the original
symmetric matrix this must be taken into account so instead of
we need QPV . The cyclic ordering must then be reversed and finally the eigenvectors
are given by
Pre-multiplication by the permutation P involves a re-ordering of the coefficients of
the eigenvectors which is carried out locally on the processor and does not involve
any communication.
5. Timings
The parallel one-sided reduction to tridiagonal has been implemented and tested
on the Fujitsu VPP500. The VPP500 is a parallel supercomputer consisting of vector
processors connected by a full crossbar network. The theoretical peak performance
of each PE is 1.6 Gflops and the maximum size of memory for each processor is
Gbyte. The VPP500 is scalable from 4 to 222 processing elements but access
was available only to a 16 processor machine. Each processor can perform send
and receive operations simultaneously through the crossbar network at a peak data
transfer rate of 400 Mbytes/s each.
The one-sided algorithm is particularly suited to the architecture of the VPP500
because the calculation of the elements of the updated matrix A from the current
version of C are vectorisable with loops of length n, the size of the input matrix.
In the conventional two-sided Householder reduction the vector length decreases at
each iteration.

Table

shows some timings and speeds obtained on up to a 16 processor VPP500.
The two times and speed given are, first, for the reduction without accumulating
the transformations for later eigenvectors calculations and, second, for the reduction
with accumulation. The speed is given in Gflops. The code was written in the
parallel language VPP-Fortran which is basically FORTRAN77 with added compiler
directives to achieve parallel constructs such as data layout, interprocessor
communication and so on.
?From these performance figures it appears that the algorithm is scalable in so far
as its performance is maintained as the size of the problem is increased along with the
14 M. HEGLAND   , M. H. KAHN    , M. R. OSBORNE
1024 3.8 .846 5.3 1.018
2048 26 1.000 37 1.164
3072 28 3.135 43 3.371

Table

1. Timings and Speed for Matrices of Increasing Size
number of processors. The one-sided reduction to tridiagonal of a matrix of size 2048
using 2 processors achieves nearly 70% of peak performance and approximately 60%
for a matrix of size 4096 using 4 processors. A formal analysis of the communication
required by the one-sided algorithm gives the communication volume proportional
to is the number of processors. As the computational load is
proportional to n 3 , isoefficiency is obtained when p(p \Gamma 1)n 2 =n 3 is constant, that is,
The scalability of the algorithm is evident if speeds for matrices of size n
on p processors are compared with those of matrices of size 4n using 2p processors.
For example, compare the speeds for
the doubling of Gflop rate.
It is interesting to compare these performance figures with other published results
for alternative parallel two-sided Householder reductions to tridiagonal. The comparisons
can only be general as the algorithms and machine architectures are very
different. The most straightforward comparison is time taken to reduce a matrix of
fixed size measured on machines of similar peak Gflop rate. In practice, the first step
of the Cholesky factorization adds an overhead of about one tenth of the time taken
for the one-sided reduction to tridiagonal. All accessible published results refer to
algorithms implemented on Intel machines.
REDUCTION TO TRIDIAGONAL FORM 15
Dongarra and van de Geijn [5] give times for a parallel reduction to tridiagonal
using panel wrapped storage on 128 nodes of a 520 node Intel Touchstone Delta.
Equating peak Mflop rates suggests that this is comparable to a 6 processor VPP500.
For a matrix of size 4000 their implementation on the Intel took twice as long as
the one-sided reduction of the same size matrix on a 4 processor VPP500.
The ScaLAPACK implementation of a parallel reduction to tridiagonal is given
by Choi, Dongarra and Walker [3]. Extrapolating from their graphs of Gflop rates it
can be seen that their times taken for various sized problems are two or three times
that taken by the one-sided algorithm on the same size problems on a VPP500 of
comparable peak performance.
Smith, Hendrickson and Jessup [10] use a square torus-wrap mapping of matrix
elements to processors and tested their code on an Intel machine corresponding to
a 12 processor VPP500. Their two-sided algorithm can be inferred to have taken
about the same time as a slightly larger problem on a 16 processor VPP500 using
the one-sided algorithm. The Smith et al algorithm is more sophisticated than the
other two-sided algorithms as it uses the torus-wrap mapping and Level 3 BLAS.
6. Conclusion
A new algorithm for reduction of a symmetric matrix to tridiagonal as the first step
in finding the eigendecomposition has been developed. Starting with the Cholesky
factorization of the symmetric matrix, orthogonal transformations based on Householder
reductions are applied to the factor matrix until the tridiagonal form is
reached. This is referred to as a one-sided reduction and leads to the updating
of the factor matrix at each iteration being rank one rather than rank two as in
the conventional Householder reduction to tridiagonal. Transformations can also be
accumulated to allow for calculation of eigenvectors. This algorithm is suited to
parallel/vector architectures such as the Fujitsu VPP500 where it has been shown
to perform well. In a complete calculation of the eigenvalues and eigenvectors of
a symmetric matrix the extra time for the Cholesky factorization is observed to
be about one tenth of that required for the reduction to tridiagonal so it is not a
significant overhead.
7.

Acknowledgement

The authors would like to thank Prof. R. Brent for his helpful comments and Mr.
Makato Nakanishi of Fujitsu Japan for help with access to the VPP500.
M. HEGLAND   , M. H. KAHN    , M. R. OSBORNE



--R


The solution of singular-value and symmetric eigenvalue problems on multiprocessor arrays
The design of a parallel
Jacobi's method is more accurate than QR
Reduction to condensed form for the eigenvalue problem on distributed memory architectures
Block reduction of matrices to condensed forms for eigenvalue computations
A fully parallel algorithm for the symmetric eigenvalue problem
Matrix Computations
The symmetric eigenvalue problem
A parallel algorithm for Householder tridiagonal- ization
The algebraic eigenvalue problem
A Parallel Ordering Algorithm for Efficient One-Sided Jacobi SVD Computations
--TR
one-sided reductions;reduction algorithms;parallel computing