At Play with J: Jacobi’s Method

VECTOR Vol.11 No.3

At Work and Play with J by Eugene McDonnell

Parallel Jacobi Warning; this column contains material which may either put you to sleep or turn you against applied mathematics altogether. To take someof the sting away I have added a problem which may give you somepleasurein tryingto solve.If you completely distrust your ability to read descriptions of programs, no matter how well-written, I advise you to go at once to the section headed Problem and avoid thepreliminary exposition, or the material following, valuableasit is.

Background Recently I had need of a program to perform eigenanalyses of square symmetric matrices, and went to Vector 9 3 for January 1993, which had Donald MclIntyre s article Jacobi s Method for Eigenvalues: an Illustration of J . I refer you to that article for McIntyre s lucid explanation of what the methodis. In the course of transcribing his 11-line jacobi program,along with its sixteen subprograms and its seven utility verbs, ] thought 1 saw the possibility of speeding it up significantly by taking advantage of some ofthe parallelism inherent in the problem. I have communicated with Mcintyre concerning this, and he tells me that he has used this method for many years, beginning with a Fortran program which he obtained from someone many years ago,transcribing it into APL and recently, as his article shows,into J.

If you lookat his program, you will see that at the heartofit are the lines

re. ((cos,-sin),sin,cos} (ia R)} I Q=. q ip j:r [C Re. r ip R ip fir

Thefirst line amends an identity matrix conforming to the argument matrix by replacing two ofits diagonal elements and the two corresponding off-diagonal elements with a 2-by-2 rotation matrix. The elements amended are chosen by finding the off-diagonal element of maximum magnitude, say at row-column indices p,q, and inserting the 2-by-2 matrix itemsat locations (p,p), (p,q), (qp) and (q,q). This amended identity matrix r is then used with two matrix products involving R, the original argument, and Q, originally an identity matrix. Those involving R have the effect of zeroing out elements (p,q) and (q,p) of R, while leaving the eigenvalues of R unaltered. When this operation has been performed a sufficient numberoftimes, onefinds thatall of the off-diagonal elements are

111

VECTOR Vol.11 No.3

essentially zero, and that the diagonal elements are the eigenvalues of the argument matrix . Those involving Q producethe eigenvectors of the argument matrix.

The valuable book Matrix Computations by Golub and Van Loandescribes this method (section 8.5), but because the search for (p,q) is O(n%2), goes on to suggest that it might be moreefficient to select p and q in a more rigid way. For the case of a 4-by-4 argument, they suggest that p and q beselected in the following order:

NRROGCOU wWanoanrnro

and go back to the beginning, repeating until a sufficiently good solution appears. Golub and Van Loan go on to point outthat the rowsofthe (p,q) table can be arrangedin a disjoint, or non-conflicting fashion:

a b c o 1 Qo 2 0 3 2 3 1 3 1 2

andthat, in a parallel machine, separate processors can be assigned to perform the individual matrix product operations. For example, in the 4-by-4 case, two processorsare needed,so that in step A one processor could do the (0,1) case and the other processor could do the (2,3) case; in step B one processor could do the (0,2) case and the other processor could dothe(1,3) case; and similarly for step C. They point outthat this method works only for even-order matrices, but that the odd case can be handled by bordering the argument matrix on the right and at the bottom with zeros, and then dropping these excess columns at the end. Thus the rotation matrices needed would looklike this:

step A { step B { step ¢ coi sot 0 Oo { co2 o s02 Oo {| co3 0 0 803 -sOi col ° o|] 0 oO 9 o{ o 0 o 0 proci 0 o o O }-s02 0 cO2 oo} 0 Qo ° Q Qo oO ° of; oa Qo 0 QO {-s03 Qo Q cog Oo 0 a oj}; 0 o Qo oj; 0 o Oo Qo 0 Q Qo o{ 0 cis Qo B13} 0 ci2 si2 Qo proc2 0 0 c23 s23{ 0 Qo 0 oO | 0 ~s12 c1i2 o 0 @ -s23 c¢23| O -s13 oO c1i3{ 0 0 0 o

112

VECTOR Vol.11 No.3

Mycontribution enters here. I realized that one doesn t need a parallel machine to obtain the benefits of this parallel Jacobi method. One can combine the rotation matrices, since they are disjunct, as follows:

step A step B step c col soi o o co2 0 802 o cO3 oO 0 s03 -s01l col 9 Qo 9 c13 9 s1i3 Q ec12 s12 o Q oO e23 s23 -s02 0 c02 oO @ -s12 c12 ° 9 QO -823 c23 oO -s13 9 e13 -s03 Qo a 03

This technique reduces the numberof matrix products required for a matrix of size n by a factor of n%2. Thus the larger the matrix, the greaterthe savings. A 10- by-10 problem can be reduced bya factor of 5; a 100-by-100 problem bya factor of 50, and soforth.

The Problem Now wecometo the playful part. As you can see, the row-column pairs to be included at each step must somehow bederived. In the case of a 4-by-4 matrix, wesee that step A uses the pairs (0 1) and (2 3); step B uses (0 2) and (1 3); and step C uses (0 3) and (1 2). The problem is to determine a permutation z that produces the desired result. For example, for n=4 any of the following permutations will do:

@ONNBRROO POROWNON ONY BRN OR w NHHROWOUNE

If wesetz=.0 3 1 2, wecan experimentas follows:

Jas. (z8()a:(1, <:4z) i. #2 NB. all of the possible permutations 04123 0312 o23a4

Jb=.((2!#2),2)$,a NB. exhibit all the pairs of items od 23 03 12 02 a4

113

VECTOR Vol.11 No.3

Je=.(>/"1)b NB. mask shows where lead item is greater than trail oogo O14

Jd=.c |]."_1 b NB. pairs with leading smaller item on 23 o3 12 02 i3

Je=./:-d NB. pairs in ascending order o4 o2 a3 12 13 23

Problem 1: Define a verb which takes as argumenta positive even integer n and yields a permutation which, repeatedly applied to a conforming identity permutation, produces, in successive pairs of items, all possible choices of 2 items from n, with no duplications.

Problem 2: How many of the !n permutations of even order n are solutions to problem 1?

Solutions to this problem may be sent by email to eem@ipsaint.ipsa.reuter.com or by ordinary mail to Eugene McDonnell / 1509 Portola Ave. / Palo Alto, CA 94306 / U.S. A.

Principal verbs The verbs described below were written for J8. If you are using an earlier version of J you may wish to get your system upgraded. Here are the verbs making up my solution to the parallel Jacobi problem. The two verbs CEA and CEAT produce identical results, but CEA is written using the rhetorical control structures which have been addedto J recently (see mylastarticle) and CEAI uses the algebraic control structures which havebeenin J from the beginning.

Each main verb CEA and CEAI (Complete EigenAnalysis) takes as argument a square symmetric matrix A and returns two conforming matrices, the first with the eigenvalues along the diagonal, and zeros elsewhere, and the second whose columns are the eigenvectors for the corresponding eigenvalues. They each test the parity of the number of rows of A. If this is even they laminate to A a conforming identity matrix, using the utility verb IM, and then apply the subverb

114

VECTOR Vol.11 No.3

PJ to this initial argument.Ifit is odd, the action is to border A on the right and the bottom with a columnand rowofzeros, using the utility verb bz, and then to apply CEA (or CEAI) to this, and at the end removing the bottom row and rightmost columnof each matrix of the result with the utility verb ub.

CEA =. 1 Tif. (2{#y.) do. ub"2 CEA bz y- else. PJ y.,:IM y. end.? CEAI=. (PJ@{,:1M)}*(ub"2@(CEAI@bz) )@. (2: |#)}

The subverb PJ (parallel Jacobi) takes as argument an array of two square matrices. It prepares four global variables for use by hsjr: a quantity eps as the product of a globally defined tolerance tol and the Frobenius norm ofthefirst matrix, yielded by the utility verb NF; a quantity s, the number of rowsin the first square matrix;a list k, the integers from 0 to s-1; and list p, a permutation which will be used to alter the arrangementof the atomsof k, using the utility verb mxp. It then employs the verb hsjr (half of s Jacobi rotations) to the limit. At the limit, it yields the desired complete eigenanalysis of the original argument,

PJ=. 3: 0 epss:tol«NF {. y. s=:# {. ¥. k=: i. § p=:mxp s hsjr a:_ y. )

The subverb hsjr (half of s Jacobi rotations) takes as argumentan array of two square matrices. It begins by makinga rotation matrix rm, using the verb RM. This rotation matrix is used with the first matrix of the argumentto develop PJo, the next stage of the eigenvalue matrix, one which has a smaller off-diagonal norm than the previous one, and setting to zero any of its elements whichareless than or equal to the quantity eps, using the utility verb clean. Next, it uses the same rotation matrix rm with the last matrix of the argument, to develop PJ1, the next stage in the eigenvector matrix. The two matrices are laminated to give the result array.

hsjr=.3 : 0 rm=.(k=:p{k) RM {.y. Pd CCfirm)+/ .*({.y.)+/ .srm) clean eps Poi=.({iy.)+/ .*rm PJO,:Pd1 )

The subverb RM (rotation matrix) builds a parallel Jacobi rotation matrix.

115

VECTOR Vol.11 No.3

It takes as left argumenta particular permutation ofthe integers from 0 through s?1. It fashions this into a two-column table t, then reverses those rows of t in which the first atom is greater than the second atom. An array ¢s of 2-by-2 cosine-sine matrices, one for each row oft, is formed, using the verb csm. These will be used to amend a matrix of zeros in locations specified by a conforming array of 2-by-2 boxes ix, whose atoms are each a 2-atom list derived from the corresponding row oft, formed usingtheutility verb ¢P (Cartesian product). For example, if a row of t is 23, the 2-by-2 boxes corresponding to it will be:

t---t---+ [2 2}2 3] torneo nt $3 213 3] tence o nt

Finally, a matrix of zeros is formed, conforming to the right argumenty., and the positions in this corresponding to positions given by the matrices of 1x will be amended with the corresponding matrices of cs, yielding the desired parallel Jacobi rotation matrix. RM=.3 : 0

t=.((-:8),2)$x. t=.(>/"1 t)]."0 1 t cs=.y. csm"2 1 ixs.cP t cs ix}0:"0 y. )

The subverb csm (cosine-sine matrix) takes as left argument a square matrix and as right argument a 2-elementlist of indices for that matrix, the first element giving a row number and the second element giving a column number, with the row numberless than the column number.If the entry in the matrix at that row- columnposition is zero, the result will be a 2-by-2 identity matrix.If it is nonzero the result will be a 2-by-2 Jacobi rotation matrix, using the verb makecs.

cesm=.makecs* (=@(1.@2:))@.(0:=<@]{[)

The subverb makecs (make cosine-sine table) takes as left argument a square matrix andasright argumenta 2-elementlist of indices for that matrix, the first element giving a row numberandthe second elementgiving a column number, with the row number less than the column number.It yields a 2-by-2 Jacobi rotation matrix. makecs=. 3: 0

116

VECTOR Vol.11 No.3

tau=.(((<2#}. yo) (x. )-(<24{. yo) (x. Hts (<7. IX. t=. (etaupxe({taujy+4 o. tau C=.%4 0. & s=.tee (c,s),:(-s),c ) The subverb mxp (make index permutation) takes a positive even integer as argumentandyieldsa list which is a permutation of the integers from 0 through one less than the argument. The permutation is such that when applied tepeatedly to a conforminglist, none of the successivepairsin thelists are equal. mxp=.E: C. O: 5 <: , (,~ >:@].)@>:@+:@1.@<rv

Utility verbs The utility verb CP takes a list as argument and returns the Cartesian product of the itemsofthelist. CPs. {@;"1~

The utility verb IM takes as argument a matrix and yields an identity matrix having the same numberof rows. IM=. C: = C: i. #

Theutility verb NF takes a matrix argument and yields its Frobenius norm as result. NF=. C: %: [i +/ Cr , *:

The utility verb clean takes a numeric array as left argument and positive atom as right argument. It yields a conforming array as result, wherein each elementof the left argument with magnitudeless than the right argumentis replaced by zero. clean=. C* J < [: | £

Theutility verb bz takes a matrix argument andyields a similar matrix bordered on the right and below by a new column and rowof zeros. bzs. >:@$ {. ]

The utility verb ub takes a matrix argumentandyields a similar matrix with the rightmost column and bottom row removed. ub=,. 1 _18}.

117

VECTOR Vol.11 No.3

Test Information Alter the following value as desired to control accuracy and speed:

tol=,ie_6 NB. value should be in the range ie_2 to 1e_17

NB. Test matrices

JAs.2 141,12 34,4 3 6 10,:1 4 10 20 tia 2 12 3 4 13 6 10 14 10 20 Jm=,1.5 _t 0.5, 1 2 _1,:_0.5 _1 4.5 1.5 _1 _0.5 i 2 4 0.5 1 4.5 Tr=.1 1 0.5,1 1 0.25,:0.5 0.25 2 1 i 0.5 1 20.25 0.5 0.25 2

NB. test results, using tol as specified above (executed on a Macintosh)

1 CEA A 0.453835 o 9 Oo O 0.038016 °o 9 o o 2.20305 9 0 Qo 9 26.3047 0.787275 _0.308686 0.530366 0,0601868 0.163234 0.723091 0.640331 0.201173 0.532167 _0.59455 0.391833 0.458082 0.265358 0.168411 0.393897 0.863752

CEA hm 2 0 0 0 3 0 0 0 0 0.707107 _0.408248 0.57735 9.8829e_10 ©,816497 0.57735 0.707107 _0.408248 0.57735 CEA r _0.0166473 o 0 Oo 1.480142 0 oO QO 2.53653 0.721208 0.44428 0.531483 0.686348 0.562411 0.461473 0.093729 _0.697601 0.710329

118