A completely parallel algorithm for the symmetric eigenproblem AX = Lambda BX is outlined. The algorithm is parallel in the sense that the numerical operations do not occur in a fixed sequence. Therefore, a large number of operations can be programmed to be performed concurrently on a computer with multiple central processing units. The standard symmetric eigenvalue problem AX = Lambda X has the property that the n eigenvalues of the principal submatrix of A of order n are separated by the (n-1) eignvalues of the principal submatrix of order (n-1). The separation property delineated n intervals containing one eigenvalue. Each eigenvalue and corresponding eigenvector can be computed independently. The n eigenproblem calculations can be divided among multiple processing units

Thurston, Gaylen A.

English

NASA Technical Reports Server

I 
- i . 
NASA Technical Memorandum 89082 
A PARALLEL SOLUTION FOR THE SYMMETRIC EIGENPROBLEM 
GAYLEN A ,  THURSTON 
(hASA-TM-890d2) 
S Y M M f i 3 R I C  E I G E N € & C E L E E  (NASA) 
A P A i i P L L E I  S C I U I I C N  FCH TBE 
14  p CSCL 2 O K  
Nd7-170Et j  
JANUARY 1987 
National Aeronautics and 
Space Administration 
Langley Research Center 
Hampton. Virginia 23665 
https://ntrs.nasa.gov/search.jsp?R=19870007653 2020-03-20T12:47:51+00:00Z
INTRODUCTION 
? 
This memorandum outlines a completely parallel algorithm for the symmetric 
eigenproblem AX - XBX. 
numerical operations do not occur in a fixed sequence. 
computer with multiple central processing units. 
The algorithm is parallel in the sense that the 
Therefore, a large 
-..-La- -C  r r r r - r t . i n - o  c=p- ).e r \ r r r n r n m m n A  t- be p-,rfcr-~,A ~ c ~ ~ ~ ~ r r ~ , p - t l _ v  22 2 
L L U I y Y b A  "I "pr&ursrr.- r'..b-".c.Y-- 
The standard symmetric eigenvalue problem AX = AX has the property that 
the n eigenvalues of the principal submatrix of A of order n are 
separated by the (n - 1) eigenvalues of the principal submatrix of order 
(n - 1). The separation property delineates n intervals containing one 
eigenvalue. Each eigenvalue and corresponding eigenvector can be computed 
independently. The n eigenproblem calculations can be divided among multiple 
processing units. 
A sufficient condition for the separation property to apply to the general 
symmetric problem AX - XBX is presented in this memorandum. The sufficient 
condition is slightly less restrictive than a condition which ensures that all 
the eigenvalues are real, namely, that B be positive definite. 
A parallel algorithm is readily derived for problems that obey the 
separation property. Once the eigenproblem of order n is solved the algorithm 
solves the eigenproblem of order (n + 1) and continues until n - N, the 
order of the square matrices A and B. Only one pass through the principal 
submatrices is required to solve for all the eigenvalues and eigenvectors of 
AX - XBX. 
Intermediate numerical results for eigenvectors must be stored. These 
results can be shared by the central processing units, or the results can be 
partitioned into separate blocks of data accessed by a single processor. 
The algorithm is numerically stable, avoids factoring matrix B, and 
automatically arranges the eigenvalues in nondecreasing order. Multiple 
eigenvalues and closely spaced eigenvalues do not present any numerical 
difficulties. 
The separation property is regarded in the literature as a theoretical 
property rather than the basis for a practical algorithm. The reason is the 
large number of matrix-vector multiplications required by the algorithm. 
total mumber of scalar multiplications is on the order of N /4 compared to 
the comparable number of 5N /3  for the Householder method for the standard 
problem. However, the input-output operations are straightforward for the 
parallel method so that the total computer cost is not necessarily as 
prohibitive as indicated by the count of scalar multiplications. 
The 
4 
3 
The body of the memorandum outlines the derivation of the separation 
property for the general symmetric problem AX = XBX. The steps in the 
algorithm are next written explicitly. 
for avoiding round-off errors in enforcing the separation property in a code 
for parallel processing is discussed. 
multiple eigenvalues are included in the discussion. 
Finally, some of the numerical analysis 
Special cases of uncoupled modes and 
2 
SEPARATION PROPERTY 
The separation property follows from examining the characteristic equation 
for the bordered eigenvalue problem of order n in the form 
. . .  ... ... ... . . .  
0. 0. ... (Pi4) ... 0. ( ci - Adi 1 
The pi are eigenvalues of the problem A (n-l) X - AB("-')X where the 
superscripts in parentheses denote the order of a principal submatrix. 
bordered form in equation (1) follows from a similarity transformation using 
U(n-l) , the matrix of orthonormal eigenvectors corresponding to the eigenvalues 
The 
Pi. 
A complete proof of the separation property for the standard eigenvalue 
problem is derived in ref. (1). 
vanish and b is unity. The effect of these terms on the proof of the 
separation property is shown here. 
For the standard problem, the coefficients di 
nn 
The characteristic equation for the matrix of coefficients in equation (1) 
can be written as 
The separation property is that the roots of equation (2) satisfy the 
conditions 
where the p are arranged in nondecreasing order. i 
A sufficient condition for the separation property to hold is derived by 
showing that f(X) 
infinity to negative infinity, 
slope of f(A) never changes sign and is always negative. 
has only one real root between consecutive values of the 
as indicated schematically in Fig. 1. Since f(A) goes from positive pi 
f(X) has only one root in each interval if the 
3 
n-1 -2di(ci - Adi) (ci - Adi)’ 
1 2 + 
(Pi - (Pi - 
- E  [ 
df 
d X  ‘bnn i-l 
or; ranrrangjnt tp-mn, 
A i-1 i-1 ( p  
(4) 
Therefore, a sufficient condition for the separation property to hold is that 
the right-hand side of equation (5) is negative. 
quadratic terms inside the square brackets, a more restrictive condition is 
determined that holds independent of the value of 
By dropping the positive 
A .  
- b  - C d f > O  
nn i-1 An 
However, the term An 
matrix 
sufficient condition for all of the eigenvalues of AX = XBX to be real. 
is a factor in a product equal to the determinant of 
B(n) and condition ( 6 )  holds when matrix B is positive definite, a 
The result in equation (5) is not a complete proof of the separation 
property, conditions ( 3 ) .  The special cases of multiplicity of equal 
eigenvalues p 
but these cases are included in the discussion in ref. (1) for the standard 
symmetric eigenvalue problem and are omitted here. The special cases are 
included in the algorithm presented in the next section. 
and of uncoupled modes where ci - di - 0 must be considered, i 
The significant result for parallel processing is that the separation 
property allows concurrent computation of the roots of the characteristic 
equation, equation (2), and of the corresponding eigenvectors from equations 
(1). 
complete parallel algorithm for the problem AX - XBX. The concurrent computation will be examined further after outlining the 
COMPLETE PARALLEL ALGORITHM 
The concurrent solution of the characteristic equation f(X) - 0, equation 
(2), is one step in the parallel algorithm. The steps in solving the nth order 
problem A (n)X - XB(n)X are the following: 
1. Given p the eigenvalues of A (n-l)X = XB(n-l)X, and the corresponding i’ 
matrix of eigenvectors 
(l), of order n. 
U(n-l), compute the bordered form, equations 
(n-1) n- 1 
k-1 uki bkn 
di - X 
i - 1,2,3, ..., (n-1) 
i - 1,2,3, ..., (n-1) 
2. Solve the characteristic equation for the bordered form, equation (l), 
for the eigenvalues X 
j' 
3. Compute the elements of the matrix of orthonormal eigenvectors for the 
bordered problem, V . (n) 
vij "ijVnj 
where 
sij - - (ci- X j  di)/(pi - Xj) i - 1,2,3,.  . ,n-1 
j - 1,2,3 , . . . ,  n 
The nth element of each eigenvector, 
the orthogonality relation containing B as a weighting matrix, 
v is determined by normalizing 
nj 
(9b) j - 1,2,3 , . . . ,  n - 1/2 n- 1 v - [ b n n + X  Sij(Sij + 2 dill 
nj i-1 
4. Compute the matrix of orthonormal eigenvectors, 
A(")X = X B(n)X. 
U(n), for the problem 
i - 1,2,3, . . . ,(  n-1) 
j - 1,2,3 ,..., n 
(n+l> (n) The complete algorithm updates the order n to (n + l), pj - A j  , 
and repeats the four steps above until 
problem AX - XBX. This completes the algorithm. The problem is solved for 
all the eigenvalues and corresponding eigenvectors with the eigenvalues 
n - N, the order of the original 
5 
automatically arranged in 
orthonormal with B as a 
I 
All the computations 
and can be programmed for 
doubly parallel since the 
nondecreasing order. The eigenvectors are 
weighting matrix, 
in the four steps with subscript j are independent 
computation in parallel. 
row operations with subscript i are also 
The algorithm is almost 
independent for each value of i. Any parallel processing with i subscripts 
depends on the data storage of the matrices and on their 
accessability by the different central processing units. The possibility of 
parallel processing with the row operations on subscript is merely noted 
here; the rest of this memorandum examines the parallel operations with the 
subscript . 
U (n-l) and lJ(n) 
i 
j 
NUMERICAL STABILITY 
Separation Property 
The algorithm summarized by the four steps in the preceding section is 
parallel in theory. In practice, the linear independence of the computed 
intermediate eigenvectors must be preserved without too much round-off error. 
The linear independence is achieved by enforcing the separation property, 
computing accurate roots of the characteristic equation, equations (2 ) ,  and 
avoiding zero or small divisors in the eigenvector computation, equations ( 9 ) .  
In programming the parallel algorithm, the separation property reduces to 
programming each central processing unit to solve equation (8) for a distinct 
X such that 3 
The characteristic equation, f(Xj) - 0, is nonlinear and must be solved 
iteratively. 
condition (12) is based on finding roots of another function, G ( X ) ,  
An iteration sequence that singles out the root that satisfies 
The 
The 
The 
function G(X)  can be written as the difference of two functions, 
first function is a cubic in A ,  
second function is 
6 
If the inequality (6) is satisfied, the cubic 
three real roots. At least one of the three roots, X = i, satisfies condition 
(12), as shown schematically in Fig. 2. The figure also shows the function 
g2(X) and X the exact root of G(X) that satisfies conditions (12). The 
zeroth approximation i for X can be improved by solving the cubic equation 
in AX generated by replacing X with 
gl(X) is guaranteed to have 
j' 
j 
in G(X) and retaining up through linear terms in AX in the Taylor's series 
for the terms in brackets in The iteration can be continued by 
updating i with X from equation (17) until AX approaches zero. 
g2(X). 
Checking The Separation Property 
At some point before the end of step 4, the results from the different 
central processing units for the eigenvalues X 
cases. There are two types of special cases. 
must be checked for special 
j 
The first special case is uncoupled modes where 
c = d I - O  I 
For this case, the complete ordering of the X in conditions ( 3 )  must be 
j 
but not X I+1 = 5- PI. checked so that XI - pI or XI+1 = pI, 
The second special case for the numerical analysis is multiple values of the 
recognize the correct order of the roots, and also avoid the apparent division 
by zero in equation (16) when the number of multiple values of 
two. Both special cases must be flagged at some point to prevent division by 
zero in computing the eigenvectors of the bordered problem, equations ( 9 ) .  
= p ,  then X = p .  The iterative solution for G(X) must 
j Pi. When Pj -1  = Pj 
p exceeds i 
Further Check on the Separation Property 
Before examining the eigenvector calculation for the special cases, two 
overall checks on the solutions of the characteristic equations are listed 
here. 
n. 
The sum and the product of the roots can be monitored at each value of 
Given the sum and the product of the eigenvalues for the problem of order 
(n-1) 9 
7 
n- 1 
i-1 'i 
S(n-1) - C 
an independent check is available for the sum and product of the next order 
n 
j -1 j 
S(n) - C X 
n 
j -1 j 
P(n) - II X 
The checks on equations (19) f o l l o w  from clearing the characteristic equation, 
f(X) 9 0, of fractions and writing it as an nth degree polynomial. The 
independent checks are then 
n- 1 
s(n) - S(n-1) + [ann + C di(pidi - 2 ci)l/An 
i-1 
is defined in equation (6). 'n where 
Linear Independence of Eigenvectors 
Two special cases were identified above that must be flagged to prevent 
division by zero in computing eigenvectors of the bordered problem. To avoid 
round-off errors in the eigenvector computation, a second subscript I(j) can 
be defined that is related to every eigenvalue by the condition that the 
absolute value of ( p  - A  ) is a minimum for the set of values (pi-X.). The 
quantity (pI- X . )  appears in the iterative solution of G ( X j )  - 0 and can be 
computed at that time rather than performing the indicated subtraction when X 
(cI- X.d ) in 
J I  
is nearly equal to 
absolute value, an alternate form of equations (9) is preferable for the 
eignvector computation, 
X 
j 
I j  J 
J 
j 
Then, if ( p  - A  ) is less than 
PI * I j  
vij - sij VIj /SIj i - I i - n  
8 
The special case of an uncoupled mode with cI = dI = 0 is a limiting case of 
equations (21) with V = 1 and the remaining v = 0. The second special 
i j 
is also a limiting case of 
1 j 
case of multiple eigenvalues p j - 1  - p j  = PI 
equations (21) where equation (21c) takes the form 
Therefore, if the I ( j )  subscripts are set properly and the quantities 
independent. The matrix of eigenvectors V (n) will be orthonormal without the 
need for checking linear independence. The accurate computation of the 
quantities 
- A  ) are computed precisely, the eigenvectors will be linearly (PI j 
- X . )  is discussed further in the appendix. ('1 J 
CONCLUDING REMARKS 
A completely parallel algorithm has been presented for the symmetric 
eigenproblem AX - XBX. The algorithm is straightforward and can serve as a 
benchmark problem for testing the efficiency of concurrent central processing 
units. 
The operation count for the algorithm is relatively high, but the storage 
reqirements and input-output operations are well-ordered. The total computer 
cost for the algorithm may prove to be competitive with more conventional 
algorithms. 
This memorandum outlines the basic algorithm for dense matrices. The 
algorithm can be adapted to banded matrices. It can also be modified for 
constrained problems where the B matrix is singular. An example of the 
latter problem arises when Lagrangian multipliers are introduced to constrain 
vibration modes. The formulation with multipliers retains symmetry, but results 
in some zero elements on the diagonal of the B matrix. The zero elements act 
to alter the separation property, but the algorithm remains parallel. 
APPENDIX 
Computing Small Roots of Cubic Equations 
The body of this memorandum outlines the solution of G ( X ) ,  equation (14), 
by an iterative sequence of cubic equations. 
root of the cubic equations approaches zero. 
affecting the final result of the iteration, the standard algorithm for real 
As the iteration converges, one 
To prevent round-off errors 
9 
roots of cubic equations should be slightly modified. 
has a change of variables that can be the source of round-off error. 
The standard algorithm 
Setting up the iterative sequence also involves a change of variables 
indicated by X - + AX in equation (17). For later use in computing 
eigenvectors, the change in variables should also explicitly include two other 
expressions 
As is updated by AX to approach X at each iteration step, the 
quantities in equations (Al) are updated in a similar fashion. 
j 
The final accuracy of the iteration depends on 
sequence of cubic equations of the form 
etx3 + efx2 + e'x + e - o 3 2 1 0 
where x = A A .  At each iteration step, an accurate 
is required for the smallest root in absolute value 
computing small roots of a 
solution of equation (A2) 
or else for the root that 
gives a current value of X - that satisfies conditon (12). 
j 
The standard algorithm, ref. 2, divides equation (A2) by the coefficient 
e' to reduce it to the form 3 
x 3 + e x 2 + e x + e  0 - 0  2 1 
Next, a change of variables is introduced, 
Y = x + (e2/3) (A41 
The change of variables transforms the cubic in x to the canonical form 
3 
y + 3py + 2q = 0 
where 
2 
P - -(e2/3) + (e1/3) 
The standard algorithm for solving equation (A5) computes a discriminant 
whose sign determines the number of real roots, 
D 
(A6 1 3 2  D - P  + q  
However, when eo is 
computing D from an 
10 
small compared to e less round-off error results from 2’ 
equivalent expression 
When the discriminant D is negative, the cubic equation has real roots. The 
three real roots of the cubic in x are 
x1 - -(e2/3) - 2r cos a (A8a) 
where 
and the sign of each square root is positive. 
When 
roots and 
x1 = 
x2 = 
x -  3 
where 
u = [-q + 
the discriminent D is positive, the cubic equation has two complex 
one real root, 
-(e2/3) + U + V (A9a) 
-(e2/3) - (U+V)/2 + i (31/2/2)(U-V) (A9b) 
-(e2/3) - (U+V)/2 - i (3’I2/2)(U-V) (A9c) 
(D)1’2]1’3, V = [-q - (D)1’2]1’3, and the signs on the cube roots 
are chosen so that W - -p. 
Theoretically, the iteration sequence for computing the X looks for 
.i 
real roots, x - A X ,  from the cubic equations. However, round-off error can 
return a small positive value for D when D is near zero. When D is 
exactly zero, the cubic has at least two equal roots. The iteration sequence 
searching for real roots can be continued for random small positive 
simply dropping the small imaginary parts of the complex roots. 
D by 
The iteration sequence converges when e - 0 in equation (A3). If the 
0 
assumption is made that the computation for 
than the round-off error introduced by the change in variables in equation 
(A4), then the smallest root of the cubic in x can also be computed from 
e contains less round-off error 
0 
11 
x1 - AA - -eo/(x2x3) 
where the x 
value. 
from equations (A8)  have been reorderd in nondecreasing absolute i 
REFERENCES 
1. Wilkinson, J. H.: The Aleebraic Eieenvalue Problem. Clarendon Press, 
(Oxford) , 1965. 
2. ADDlicable Mathematics. Rektorys, K., Editor, M.1.T Press (Cambridge, Mass.) 
1969. 
12 
\ 
: 
Figure 1 .  Schematic diagram i l lustrat ing  separation 
property for roots of f ( X ) ,  equation ( 2 ) .  
13 
Figure 2 .  
of g,(A), zeros of g2(A), and A .  J where g,(A) = g,(A). 
Schematic diagram showing roots  of the  cubic function 
1. Report No. I 2. Government Accuion No. I 3. Recipimt's GUIoQ No. 1 
NASA TM-89082 I 
4 Title wid Subtitle 
A P a r a l l e l  S o l u t i o n  for t h e  Symmetric Eigenproblem 
5. R.pon Date 
6. Performing Organization C o d  
January 1987 
505-63-11-03 
7 AuthorW 
Gaylen A. Thurston 
0 Pertwiriing Orydoizalion Name dnJ Addrcrr 
NASA Langley Research Center 
Hampton, VA 23665-5225 
8. Porforming Organization Report No. 
I 
10. Work Unit No. 
f 
11. Contract or Grant No. 
This memorandum o u t l i n e s  a completely p a r a l l e l  a lgori thm f o r  t h e  symmetric 
eigenproblem AX = ABX. 
numerical ope ra t ions  do not occur  i n  a f i x e d  sequence. Therefore,  a l a r g e  
number of ope ra t ions  can be programmed t o  be performed concur ren t ly  on a 
computer w i th  mul t ip l e  c e n t r a l  processing u n i t s .  The s t anda rd  symmetric 
eigenvalue problem AX = AX has  t h e  p rope r ty  t h a t  t h e  n eigenvalues  of 
The algori thm is  p a r a l l e l  i n  t h e  sense t h a t  t h e  
13. Type of Rapon and Period Covered 1 - 
1 7  K e y  Words lsuggertctd by AuthocIr) 1 
p a r a l l e l  s o l u t i o n  
I ?  SpwiroiIrlg Ayency Name and Address 
18. Distribution Statement 
U n c l a s s i f i e d  - Unlimited 
Nat ional  Aeronautics and Space Adminis t ra t ion 
Washington, DC 20546-0001 
I 
Unclas s i f i ed  14 A02 
symmetric eigenproblem Subject  Category - 39 


A parallel solution for the symmetric Eigenproblem

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server