A SYSTOLIC ARCHITECTURE FOR THE SINGULAR VALUE DECOMPOSITION Richard P. Brent* Franklin T. Luk TR 82--H522 September 1982 Department of Computer Science Upson Hall Cornell University Ithaca, New York 14853 *Departmentof Computer Science, Australian National University, Canberra, A.C.T. 2600, Australia Supported in part by the U.S. Amy Office under grant DAAG 29-79-C0124 and in part by the Mathematical Sciences Research Centre and the Centre for Mathematical Analysis, Australian National University. A Systolic Architecture for the Singular Value Decomposition Richard P Brent Department of Computer Science Australian National University Canberra, A.C.T. 2600 Australia Franklin T. Luk* Department of Cpmputer Science Cornell University Ithaca, N.Y. 14853 U.S.A. TR--HCS-82--H09 August, 1982 * Supported in part by the U.S. Army Office under grant DAAG 29-79-C0124 and in part by the Mathematical Sciences Research Centre and the Centre for Mathematical Analysis, Australian National University. Abstract We propose a systolic architecture for computing a singular value decomposition of an m >C n matrix, where m >- n Our algorithm is stable and requires only 0 (mn) time on a linear array of 0 (n) processors. Extensions to algorithms for two--Hdimensional arrays are also discussed. KeyWordsandPhrases: Systolic arrays, singular value decomposition, Hestenes method, threshold Jacobi method, real-time computation. 1. 2 Introduction A singular value decomposition (svD) of an m X n (m > n) matrix A is its factorization into the product of three matrices: (1.1) A = u ? vT where U is an m X n matrix with orthonormal columns, E is an n x n nonnegative diagonal matrix and the n X n matrix V is orthogonal. This decomposition has many important scientific and engineering applications (cf. [1], [5] and [11]), for some of which real-time computation is desirable. With the advent of VLSI technology, there is interest in a systolic network for computing the decomposition. In [3] Finn et al. describe one such architecture and two related algorithms. Unfortunately, the convergence of their algorithms is only an empirical observation and has not been proved. Heller and Ipsen [7] consider the problem of computing only the singular values of a banded matrix with bandwidth w They present a systolic array of 0(w) processors and an O(wn2) algorithm. In this paper we present a systolic array of 0(n) linearly--H connected processors which computes the SVD in time O(mn) Alternatively, a two--Hdimensional systolic array of O(mn) processors with a more complicated interconnection pattern can compute the SVD in time O(nlogm) . Our systolic architectures implement a one--Hsided orthogonalization method due to Hestenes [8]. The Hestenes method is essentially the serial Jacobi procedure T for finding an eigenvalue decomposition of the matrix A A , and has been used successfully by the second author [9] on the ILLIAC IV computer. A complication arises here in that, for the sake of parallel computing, we discard the classical scheme of rotating column pairs in the order (l,2),(l,3), ..., (l,n),(2,3), ..., (n--Hl,n) To enforce convergence, we choose a threshold approach (see Wilkinson [13, pp. 277--H278]). 2 Hestenes method The idea is to generate an orthogonal matrix V such that the transformed matrix AV has orthogonal columns. Normalizing the euclidean length of each nonnull column to unity, we get the relation (2.1) AV = U ? where U is a matrix whose nonnull columns form an orthonormal set of vectors and Z is a nonnegative diagonal matrix. A SVD of A is thus given by (1.1') A = ? vT U As a null column of ? is always associated with a zero diagonal element of ? , there is no essential difference between (1.1) and (1.1') Hestenes [8] suggests that the orthogonal matrix V should be constructed as a sequence of plane rotations, say Q1,Q2Let (2.2) Ac A Q1 Q2 WitI?out loss of generality, we may assume that the next rotation acts on the i--Hth and j--Hth columns ?c? and ?c , ? < j) of the current matrix Ac . Let us consider (2.3) ?ac.,?c.?[_cosO sinO NN nO cosO, (?a???aj) Our problem is therefore choosing a rotation angle 6 to generate N N an orthogonal pair of column vectors a. and a . We use the Th. ?1'7?T' ?V Thir?s;hauser FlOl. Definin? 4 (2.4) -= II?II22 , 11Aj1122 and y we set 6 = 0 if Y = 0 , and otherwise compute (2.5) = 2y = sign(?) l?l + cosO = 1 and sin6 = t.cos6 The rotation angle always satisfies (2.6) 161 < --H? which would guarantee convergence had we rotated the columns in the classical order of (l,2),(l,3), ..., (l,n),(2,3), ..., (2,n),(3,4), (n-l,n) (Forsythe and Henrici [4]). However, in order to perform the rotations in parallel, we choose a scheme (to be described in Section 4) in which all column pairs are rotated just once in any sequence (called a "sweep") of n(n-l)/2 rotations, but not in the above order. Very little is known about the convergence properties of this alternate scheme (see Hansen [6]), although we note that it has an asymptotically optimal "preference factor" [6, eqn. (23)] To enforce convergence, we adopt a threshold approach (Wilkinson [13, pp. 277-278]). Let us associate with each sweep a threshold value, and when making the transformations of that sweep, we omit any rotation based on a normalized inner product IJ?II2II?jII2 5 which is below the threshold value. Our method enjoys ultimate quadratic convergence (Wilkinson [12]) and numerical experience suggests that only six to ten sweeps are required (Rutishauser [10]). We may compute the matrix V by accumulating the plane rotations. If we let (2.7) vc Q1Q2 ?k (?v?1, o+??, c then the (k+l)--Hst rotation will affect only the columns vc. and vc We can therefore determine the new columns using (2.8) - cosO sinO v. -sinO cosO v v c c i.e. we update both A and V simultaneously 3. Generation of all airs (i,?) In this section we show how 0(n) linearly-connected processors can generate all pairs (i,j), 1 < i < j -< ii , in 0(n) steps. The application to the computation of the SVD is described in Sections 4 and 5. First suppose n is even. We use n/2 processors Pl, ?n/2 ` where ?k and ?k+l communicate (k=l,2,...,n/2-l). Each k d , output lines outL? and processor ?k has registers L an out? , and input lines inL? and in? , except that outL1 , inL1' outRni2 and inRni2 are omitted. The output outP???is connected to the input inL?+1 and the output out?+1 is connected to the input? in? , as shown in Figure 1 6 L1 R1 Pi outR1 inL2 inR1 outL2 outR2 inL3 L2 R2 inR2 outL3 L3 R3 P2 P3 Fi ure 1: Inter--H rocessorconnectionsforn=8 Initially Lk = 2k - 1 and Rk = 2k executes the following program: outR3 inL4 inR3 outL4 L4 R4 P4 At each time step processor if Lk < then process (Lk,Rk) else process if k=l then out?:=? else if k < n/2 then if k > I then outLk:=Rk; fwait for outputs to propagate to inputs of adjacent processors1 if k < n/2 then ?:=in? else if k > 1 then Lk:=inL?; Here process (i,j) means perform whatever operations are required on the pair (i,j), 1 < i < j < n The operation of the systolic array is illustrated in Figure 2. We see that the index 1 sta?s in the register L of processor P1 . Indices 2, ..., n travel through a cycle of length n--Hl consisting of the registers L2,L3, ..., Ln/2?Rn12? Rn/2?l? ..., R1 During any n--Hl consecutive steps a pair (i,j) or (j,i) can appear in a register pair (Lk?Rk) at most once. A parity argument shows that (i,j) and (j,i) can not both occur (see Figure 2). Since there are n/2 register pairs at each of n--Hl time steps, each possible pair (i,j) , 1 < i < j < n , is processed exactly once during a cycle of n--Hl consecutive steps. 7 ste 0 ste 1 ste 2 ste 3 ste 4 ste 5 ste 6 1 Z\ 3?% Il4)? _____ I 6 ______ II 8:.' IIm? M [ffi3j 13 ______ 7) 8 N jm.?8t N[ ? o+?2A3 1$.N$?[ [6..?9?j _____ 82 t64 2 7 4 [8 I Fl ure. 2: Full c cle of the 5 stolic arra for n = 8 If n is odd, we use [n/2i processors but initialize Lk = 2k - 2 , = 2k - 1 for k=l, ..., ??n21 and omit any ??process1V calls from processor P1 4. A one--Hdimensional systolic array for SVD computation Assume that n is even (else we can add a zero column to A or modify the algorithm as described at the end of Section 3). We use n/2 processors ?l? o+` ?n/2, as described in Section 3, except that Lk and Rk are now local memories large enough to store a column of A (i.e. Lk and each have at least m floating-point 8 words). Shift registers or other sequential access memories are sufficient as we do not need random access to the elements of each row. c Suppose processor ?k contains column ? in L and column 1 k ac in . It is clear that ?k can implement the transformation Rk defined by (2.3)-(2.6) in time 0(m) by making one pass through c and c a. a to compute the inner products (2.4), and another pass to perform the transformations (2.3) or (2.8). Adjacent processors can then exchange columns in the same way that the processors of Section 3 exchange indices. This takes time 0(m) if the bandwidth between adjacent processors is one floating-point word. (Alternatively, exchanges can be combined with the transformations (2.3) or (2.8).) Consequently, we see that n/2 processors can perform a full sweep of the Hestenes method in n --H 1 steps of time 0(m) , i.e. in total time 0(mn) . Initialization requires that the (2k-l)-th and 2k--Hth columns of A be stored in the local memory of processor ?k for k=l, ..., n/2; clearly this can also be performed in time 0(mn) Although the process is iterative, in practice only 6--H10 sweeps are required to orthogonalize the columns to full machine accuracy (see Section 2). Hence, we have a systolic array of n/2 processors which computes the SVI) in time 0(mn) By comparison, the serial Hestenes algorithm takes time 0(mn2) After an integral number of sweeps the columns of the matrix AV = U? (see (2.1)) will be stored in the systolic array (two per processor). If V is required, it can be accumulated at the same time that AV is accumulated, at the expense of increasing each processor's local memory (but the computation time remains 0(mn) ) see (2.8). Two--Hdimensional 5 stolic arra 5 for SVI) corn utation 9 5 It is natural to ask if a rectangular array of m by n/2 processors can be used to obtain a greater speedup than a linear array of n/2 processors. The only difficulty (but a major one) is in the fast accumulation of the inner products (2.4). If the processors have a binary tree interconnection pattern superimposed on each column of the rectangular array, so that the inner products (2.4) can be accumulated in time O(log m) , then an array of mn/2 processors can compute the SvD of A In time O(nlog m) , by an obvious modification of the method of Section 4. As observed in Section 1, Hestenes method is theoretically T equivalent (with exact arithmetic) to the JacobI method applied to A A In [2j we show that Jacobi's method can be implemented in time 0(n) on a square array of n/2 by n/2 systolic processors with nearest- neighbour connections. Hence, the factor log m above (and the non- local connection pattern) can be avoided, at the expense of some loss of numerical accuracy. References 10 [1] H.C. Andrews and C.L. Patterson, "Singular value decomposition and digital image processing',, IEEE Trans. Acoustics, Speech and Signa? Processinj ASSP-24 (1976), 26-53 [2] R.P. Brent and F.T. Luk, 1'A systolic architecture for linear time solution of the symmetric eigenvalue problem", Tech. Report TR-CS-82-1O, Dept. of Computer Science, Aust. Nat. Univ., to appear. [3] A.M. Finn, F.T. Luk and C. Pottle, `1Systolic array computation of the singular value decomposition", Proc. SPIE Symp. East 1982, Vol. 341, Real Time Signal Processing V (1982), to appear. [4] G.E. Forsythe and P. Henrici, "The cyclic Jacobi method for computing the principal values of a complex matrix", Trans. Amer. Math. Soc. 94 (1960), 1-23. [5] G.H. Golub and F.T. Luk, "Singular value decomposition: applications and computations", ARO Report 77--H1, Trans. of 22nd Conf. of Army Mathematicians (1977), 577-605. [6] E.R. Hansen, "On cyclic Jacobi methods", J. Soc. Indust. Appi. Math. 11 (1963), 448--H459 [7] D.E. Heller and I.C.F. Ipsen, "Systolic networks for orthogonal equivalence transformations and their applications", Proc. 1982 Confi. on Advanced Research in VSEI, MIT (1982), 113-122. [8] M.R. Hestenes, `Inversion of matrices by biorthogonalization and related results", j. Soc. Indust. Appi. Math. 6 (1958), 51--H90. [9] F.T. Luk, "Computing the singular-value decomposition on the ILLIAC IV", ACM Trans. Math. Softw. 6 (1980), 524--H539. 11 [10] H. Rutishauser, "The Jacobi method for real symmetric matrices", in Handbook for Automatic Computation, vo?. 2 (Linear A?o?ebra) eds. J.H. Wilkinson and C.Reinsch, Springer--HVerlag, Berlin, (1971), 202--H211. [11] J.M. Speiser and H.J. Whitehouse, "Architectures for real time matrix operations," Proc. Covernment Mtcrocircuit Applicaftons Conf., Houston, Texas (1980). [12] J.H Wilkinson, "Note on the quadratic convergence of the cyclic Jacobi process", ivumer. Math. 4 (1962), 296--H300 [13] J.H. Wilkinson, The A1?ebraic Ei?enva?ue Prob?em, Clarendon Press, Oxford (1965).