Modular exponentiation of matrices on FPGA-s

Acta Univ. Sapientiae, Informatica, 3, 2 (2011) 172–191 Modular exponentiation of matrices on FPGA-s Tam´as HERENDI University of Debrecen email: [email protected] Roland S´andor MAJOR University of Debrecen email: [email protected] Abstract. We describe an efficient FPGA implementation for the exponentiation of large matrices. The research is related to an algorithm for constructing uniformly distributed linear recurring sequences. The design utilizes the special properties of both the FPGA and the used matrices to achieve a very significant speedup compared to traditional architectures. 1 Introduction Field-programmable gate arrays (FPGA) offer a number of special options in computation. Utilizing the unique properties of an FPGA, some algorithms that are impractical to implement on a more traditional architecture can become both convenient to create and resource-efficient. The programmable array of look-up tables commonly found on an FPGA provide both flexibility in creating logic to suit specific needs and naturally lend themselves to great parallelism in computations. Fast operations on matrices are of great practical interest. Ways to speed up certain matrix calculations still find their way into numerous applications. Faster implementations of matrix algorithms can be achieved either from a “software” point of view, by improving upon the algorithm itself, or from a Computing Classification System 1998: B.2.4 Mathematics Subject Classification 2010: 65F60 11Y55 Key words and phrases: pseudo random number generators, linear recurring sequences, uniform distribution, matrix exponentiation, parallel arithmetic, FPGA design, hardware acceleration of computations, hardware implementation of computations 172 Modular exponentiation of matrices on FPGA-s 173 “hardware” point of view, by using faster or differently structured architectures. Theoretical improvements on matrix algorithms include Strassen’s algorithm [12] and the Coppersmith-Winograd algorithm [2]. The naive algorithm for matrix multiplication is a well-known Θ(n3)algorithm. Strassen’s algorithm uses an idea similar to the Karatsuba-multiplication. It has a time complexity of O(nlg 7)by dividing the matrices into sub-matrices. Then by multiplying them in a different arrangement, it manages an overall lower multiplication count compared to the classical algorithm. Research implementing it on the Cell Broadband Engine can be found in [5]. Strassen’s algorithm and its applicability to the project is briefly discussed in Section 7. The CoppersmithWinograd algorithm further improves the complexity to O(n2.376)by combining the idea of Strassen with the Salem-Spencer theorem. [9] discusses and compares the performance of implementations of these algorithms. Numerous research has been done on creating efficient realizations of different matrix operations on different architectures. [8] and [10] both use FPGAs to perform matrix inversion. The design presented here is an implementation of matrix multiplication on an FPGA. Works of similar nature can be found in [1] and [4], dealing with FPGA configurations used for floating point matrix multiplication. [11] uses an FPGA design for digital signal processing. [3] discusses another FPGA implementation for accelerating matrix multiplication. The research in this paper is related to an algorithm for the construction of pseudo random number generators. It requires the exponentiation of large matrices to an extremely high power. This allows for numerous optimizations to be made on the FPGA implementation, resulting in an extremely fast design. A speedup factor of ∼200 is achieved compared to a highly optimized program on a more traditional architecture. We give the details of a design implemented on a Virtex-5 XC5VLX110T FPGA that multiplies two 896 ×896 sized matrices. The matrices are defined over the mod 4 residue class ring. Using this property and the fact that the hardware uses 6-LUTs (Lookup Tables), we describe first a module that computes the dot product of vectors taken from Z28 4in a single clock cycle at 100MHz clock speed. With these modules we construct a matrix multiplier module that computes the C∈Z20×20 4product matrix of A∈Z20×28d 4and B∈Z28d×20 4in dclock cycles at 100MHz. The significance of the value 28 in the implementation and its experimental determination is also discussed. Finally, we describe how to use these modules for multiplying matrices taken from Z896×896 4. The proposed algorithm deals with the management of stored 174 T. Herendi, R. Major data in such a way that it can be accomplished completely in parallel with the computations. The resulting design completes the multiplication in 64800 clock cycles at 100MHz. Future work for increasing the size of the used matrices, and further optimizing the design’s performance using Strassen’s algorithm is also described. 2 Mathematical background The present work is initiated by a method for the construction of uniformly distributed pseudo random number generators. (See [7].) The generator uses recurring sequences modulo powers of 2 of the form un≡ad−1un−1+ad−2un−2+· · · +a0un−dmod 2s, ai∈{0, 1, 2, 3}, s ∈Z+ The theoretical background can be found in [6]. The construction assumes that the values a0, a1,...,ad−1are such that xd−ad−1xd−1−· · · −a0≡(x−1)2P(x)mod 2 holds for some P(x)irreducible polynomial. It is practical to choose P(x)to have maximal order, since the order of Pis closely related to the period length of the corresponding recurring sequence. The sequence unobtained this way does not necessarily have uniform distribution, however exactly one of the following four sequences does: u(0) n≡ad−1u(0) n−1+ad−2u(0) n−2+· · · +a1u(0) n−d+1+a0u(0) n−dmod 2s u(1) n≡ad−1u(1) n−1+ad−2u(1) n−2+· · · +a1u(1) n−d+1+ (a0+2)u(1) n−dmod 2s u(2) n≡ad−1u(2) n−1+ad−2u(2) n−2+· · · + (a1+2)u(2) n−d+1+a0u(2) n−dmod 2s u(3) n≡ad−1u(3) n−1+ad−2u(3) n−2+· · · + (a1+2)u(3) n−d+1+ (a0+2)u(3) n−dmod 2s. For the details see [7]. Finding the sequence with uniform distribution is of interest. Let M(u) =        0 1 . . . 0 0 . . .. . ..... . .. . . 0 0 . . . 1 0 0 0 . . . 0 1 a0a1. . . ad−2ad−1        Modular exponentiation of matrices on FPGA-s 175 be the companion matrix of sequence u. To find which of the above sequences has a uniform distribution, we have to compute M(u)2d+1−2mod 4. If M(u)2d+1−2mod 4equals the identity matrix, then the period length of un is 2d+1−2, which means it is not the sequence we are searching for. The exponentiation of matrices to high powers can quickly become time consuming on traditional computers. The aim of the project was to utilize the special properties of an FPGA to achieve a significant upgrade in speed compared to implementations on more traditional architectures. 3 Hardware used in the implementation The project was implemented on a Xilinx XUPV505-LX110T development platform. The board features a variety of ports for communication with the device. As a first approach the RS-232 serial port was used to send data between the board and a PC. A high-speed PCI Express connection is also available if the amount of data transferred would necessitate its use. The board’s most prominent feature is the Virtex-5 XC5VLX110T FPGA. The FPGA’s main tool for computation is the array of 6-input look-up tables, arranged into 17280 Slices, with four look-up tables found in each Slice, adding up to a total of 69120 LUTs. A single 6-input LUT can store 64 bits of data, where its six input bits are used as an address to identify the single bit of data that is to be outputted. By manipulating the 64 bit content of the look-up table, it can be configured to carry out arbitrary Boolean functions with at most six input bits. In our design they are used to create LUTs performing a multiply-accumulate function, which are hierarchically arranged into larger and more complex modules. One out of four LUTs on the device can also be used as a 32 bit deep shift register; these are the basis to implement containers storing the data, which is directly fed to the computational module. Attached to the board, there is a 256MB DDR2 SODIMM module, which is used for storing data exceeding the amount that can be practically stored on the FPGA. 4 Structure of modules used in the computation The basic elements of the design are the LUTs denoted by L(a, b, s) = c, where a, b, c and sare two-digit binary numbers. The function carried out by Lis a multiply-accumulate (for short: MA) function, i.e.: c≡(a·b) + smod 4 . 176 T. Herendi, R. Major Let a=2α1+α0,b=2β1+β0,s=2σ1+σ0,c=2γ1+γ0, where α0, α1, β0, β1, σ0, σ1, γ0, γ1∈{0, 1}, and L= (l1, l0)where l1and l0are two single bit LUTs, according to the following: •l0(α0, β0, σ0) = γ0 •l1(a, b, s) = γ1 Figure 1: The structure of L(a, b, s) We remark that while l0needs only three input bits to accomplish its function, l1requires all six bits of input. The LUTs l0and l1were configured to the values shown in Table 1 and Table 2 to perform the multiply-accumulate function. PPPPPPPP P (α0, β0) σ00 1 (0,0) 0 1 (0,1) 0 1 (1,0) 0 1 (1,1) 1 0 Table 1: Contents of l0 With the help of these basic units one can compute the dot product wof two vectors u= (u0, u1,...,un−1)and v= (v0, v1,...,vn−1). Let us define a module m= (L[0], L[1],...,L[n−1]) by cascading nMA units denoted by L[i]. In this module mwe use the output of a given MA unit as the sum input of the next unit, i.e. si+1=cifor i=0, 1, . . . , n −2, where siand ciare the s input and coutput of L[i]. Therefor mis a function that accepts a pair of vectors u, v of two-digit numbers of length nand outputs on cn−1the two-digit dot-product of the two vectors, i.e. m(u, v) = w. Modular exponentiation of matrices on FPGA-s 177 PPPPPPPP P (a, b) s0123 (0,0) 0 0 1 1 (0,1) 0 0 1 1 (0,2) 0 0 1 1 (0,3) 0 0 1 1 (1,0) 0 0 1 1 (1,1) 0 1 1 0 (1,2) 1 1 0 0 (1,3) 1 0 0 1 (2,0) 0 0 1 1 (2,1) 1 1 0 0 (2,2) 0 0 1 1 (2,3) 1 1 0 0 (3,0) 0 0 1 1 (3,1) 1 0 0 1 (3,2) 1 1 0 0 (3,3) 0 1 1 0 Table 2: Contents of l1 Figure 2: The structure of m(u, v) In total, the number of LUTs used in mis 2n. Note that vectors of arbitrary length can be used in the computation if we connect the output of module m to the sum input of L[0](cn−1=s0), and then iteratively shift uand vonto the module’s input by nelements at a time: 178 T. Herendi, R. Major Function iterated m(u, v)// k=length(u) = length(v) 1. Define κ=dk ne,v0, u0∈Zκ·n 4 2. for i=0to κ·n−1do // fill vand uwith 0’s 3. if i<kthen v0 i=vielse v0 i=0 4. if i<kthen u0 i=uielse u0 i=0 5. end for 6. Define vtemp, utemp, w, let w=0 7. for i=0to κ−1do // shift v0and u0to vtemp and utemp 8. vtemp = (v0 i·n, v0 1+(i·n),...,v0 n−1+(i·n)) 9. utemp = (u0 i·n, u0 1+(i·n),...,u0 n−1+(i·n)) 10. w=w+m(vtemp, utemp) 11. end for 12. return w end Function Here u0and v0are the extensions of uand vby 0’s. We shall see that the number chosen for nis critical in setting many characteristics of the entire project. The experiment used for determining nwill be discussed in the following chapter. Our aim is to obtain a module that performs the matrix multiplication of A, B ∈Zk×k 4, where Z4is the mod 4 residue class ring. In the following, let C∈Zk×k 4be the output matrix, such that C=A×B. Furthermore, let aibe the ith row of matrix Aand let bjbe the jth column of matrix B. The multiplier units denoted by mare used to create more complex modules in a hierarchical manner. First, by taking ten mmultiplier blocks we create a row of multipliers R= (m0, m1,...,m9). This is used to compute ten consecutive elements of a single row of the output matrix: R(ai, bj, bj+1,...,bj+9)=(ci,j, ci,j+1,...,ci,j+9), where ci,j =ai·bj. The input vector aiis used by all ten multiplier units of R. The length of these vectors, as mentioned above, can be arbitrary, but vectors of length greater than nwill need to be iteratively shifted to the input of R. By taking ten row multipliers we can create a unit M10×10 = (R0, R1,...,R9) which outputs a 10 ×10 sub-matrix of C: Modular exponentiation of matrices on FPGA-s 179 M10×10(ai, ai+1,...,ai+9, bj, bj+1,...,bj+9) =      ci,j ci,j+1· · · ci,j+9 ci+1,j ci+1,j+1· · · ci+1,j+9 . . .... ci+9,j ci+9,j+1· · · ci+9,j+9      . Finally, four such units are arranged so that a 20 ×20 sub-matrix of Ccould be obtained as output: M20×20(ai, ai+1,...,ai+19, bj, bj+1,...,bj+19) =      ci,j ci,j+1· · · ci,j+19 ci+1,j ci+1,j+1· · · ci+1,j+19 . . .... ci+19,j ci+19,j+1· · · ci+19,j+19      . The M20×20’s inputs are twenty vectors from both matrices Aand B. Because of hardware constraints — in particular the number of LUTs on the used device — a larger arrangement of multipliers would be impractical to implement. The module M20×20 is comprised of 400 mmultiplier units. Figure 3 shows the hierarchy of units used to build M20×20. Figure 3: The structure of M20×20 The M20×20 unit can be used iteratively to multiply matrices of arbitrary size, producing 20×20 sub-matrices of the output matrix Cwith each iteration. After inputting twenty rows from matrix Aand twenty columns from matrix Band obtaining the desired output, we can simply repeat the process for a 180 T. Herendi, R. Major set of rows and columns of Aand Brespectively, until we obtain the entire output matrix C: Function large matrix mult(A, B) 1. Define κ=dk 20 e,A0, B0, C0∈Zκ·n×κ·n 4 2. for i=0to 20κ −1do 3. for j=0to 20κ −1do 4. if i<kand j < k a0 ij =aij else a0 ij =0 5. if i<kand j < k b0 ij =bij else b0 ij =0 6. end for end for 7. for i=0to κ−1do 8. for j=0to κ−1do 9. C0[i,j i+19,j+19] = M20×20(ai, ai+1,...,ai+19, bj, bj+1,...,bj+19) 10. end for end for 11. return C0[0,0 k−1,k−1] end Function Here C0[i,j k,l] =      c0 i,j c0 i,j+1· · · c0 i,l c0 i+1,j c0 i+1,j+1· · · c0 i+1,l . . .... c0 k,j c0 k,j+1· · · c0 k,l      . Note that in the naive algorithm large matrix mult(A, B), during the main loop (lines 7-10), for each twenty rows read from A, the entire matrix Bis read. During the whole procedure, matrix Awill be read entirely exactly once, while matrix Bwill be read κtimes. Methods improving on this number are described in section 6. Since for almost all practical cases the size kof matrices A, B ∈Zk×k 4will be greater than the parameter n, the vectors taken from these matrices will need to be iteratively shifted onto the input of the multiplier M20×20,nelements at a time. Therefore, an efficient way to both store and then use the vectors taken from the matrices is the creation of FIFO type containers made of shift registers. Let td nbe a shift register of width nand depth d. It means that td ncan store at most dvectors of length n, or equivalently a single vector of length at most nd. We choose dsuch that nd ≥k, thus it can store one row or column from the input matrices Aor B. Let the vector filling td nbe f= (f0, f1, . . . , fd−1), Modular exponentiation of matrices on FPGA-s 187 New rows are loaded in at a slower pace than columns are. By the time all columns are read once, the contents of the row-stores have shifted exactly to the next segment of data needed, the next (z−1)·20 rows. After matrix Bis completely read once, the row-stores are filled with rows a(z−1)˙ 20 →a2(z−1)·20−1. Reading rows and columns proceeds in this manner until we’ve completely read matrix Aonce. For this reason, it is practical to choose zsuch that (z−1)|κ. All together we read matrix Bκ z−1times and matrix Aonce. During each z−1iterations shown in Figure 6, twenty new columns and (z−1)·20 κ new rows are loaded into the column-store and row-store currently unused by the computation. When the unused row-store is filled with twenty new rows, it becomes active, to be used in the following iterations. The row-store containing the rows with the least index becomes inactive in the computation and starts accepting the new rows read. Figure 7: Progression of computations through matrix C Function improved matrix mult(A, B) 1. Define z, κ =dk 20 e,A0, B0, C0∈Zκ·n×κ·n 4 2. for i=0to 20κ −1do 3. for j=0to 20κ −1do 4. if i<kand j<kthen a0 ij =aij else a0 ij =0 5. if i<kand j<kthen b0 ij =bij else b0 ij =0 6. end for end for 188 T. Herendi, R. Major 7. Fill the row-stores with rows a0→a(z−1)·20−1 8. Fill the column-stores with columns b0−b19 9. For i=1to κ2 z−1 Do in parallel: |perform z−1iterations of the computation |READ the next 20 columns mod κ·20 |READ the next (z−1)·20 κrows mod κ·20 |WRITE the result of the previous z−1iterations 11. return C0[0,0 k−1,k−1] end Function The possible values for the parameters used in this section depend on the used hardware. The size of the matrices used in the implementation are determined by parameters n=28 and d=32. The LUTs on the device that comprise the Td 20n containers can be configured as d=32 bit deep shift registers. For this reason the matrices are of size 896 ×896. Rows with length k=896 are the largest that can be stored in containers that are one LUT deep, making them any larger would double the number of LUTs needed for creating a Td 20n. Because of the limited number of LUTs which can be used for storage purposes, z=10 was chosen. This yields that twelve Td 20n containers are defined in the design. Dealing with matrices larger than k=896 is part of future work. For convenience, time quantities are measured in clock cycles at 100MHz, the clock speed of the M20×20 multiplier. The value of δdepends on the DDR2 RAM used. The device was used at 200MHz, and has a 64 bit wide physical data bus. From these values we determine the following parameters: •κ=d896 20 e=45, •K(improved matrix mult) = κ z−1κ+κ=45 9·45 +45 =270, •δ=140 clock cycles at 100MHz, •Φ(improved matrix mult) = K(improved matrix mult)δ+κδ =270· 140 +45 ·140 =44100 clock cycles at 100MHz, •Γ(improved matrix mult) = dκ2=32 ·452=64800 clock cycles at 100MHz. Modular exponentiation of matrices on FPGA-s 189 The goal of Γ(improved matrix mult)> Φ(improved matrix mult)is achieved, meaning that the running time of the design is equal to the time used by the computation. The speedup provided by the configuration can be shown by comparing its performance to a similar implementation created on a more traditional architecture. A highly optimized C++ program was created for a machine using an Intel E8400 3GHz Dual Core processor with 2GB RAM. The algorithm is strongly specialized for the task, making use of all available options for increasing performance. It uses 64 bit long variables to perform multiplication on 16 pairs of two-digit elements at once in parallel on both processor cores. The running time of the multiplication of matrices of the same size is over 100 ms. The FPGA implementation, as mentioned above, achieves a runtime of ∼0.6 ms. On average, a speedup factor of 200 is reached using the described FPGA design. 7 Future work The future course of research will focus on increasing the size of the used matrices. As mentioned in the previous section, simply increasing the depth dof the Td 20n containers would be impractical. Since a single LUT on the device can only be configured as a 32 bit deep shift register, setting d > 32 would double the number of LUTs needed for a Td 20n, and the design is already using well over half of the device’s LUTs that can be configured this way (13440 out of 17280, to be exact). Increasing the size of the matrices this way would require the restructuring of both the multiplier module and the algorithm used for memory management. Instead, the currently implemented module can be used as a basic unit for the multiplication of larger matrices. Then the entries of the large matrices are 896 ×896 blocks. This also allows for further optimization using Strassen’s algorithm. Suppose we double the matrix sizes, interpreting them as matrices with four blocks. Using the classical algorithm, multiplying two 1792×1792 sized matrices would take eight multiplication of the blocks. Using a divide-and-conquer strategy, we can exchange one multiplication for a few extra additions. A11 A12 A21 A22·B11 B12 B21 B22=−D2+D4+D5+D6D1+D2 D3+D4D1−D3+D5−D7, 190 T. Herendi, R. Major where D1=A11(B12 −B22) D2= (A11 +A12)B22 D3= (A21 +A22)B11 D4=A22(B21 −B11) D5= (A11 +A22)(B11 +B22) D6= (A12 −A22)(B21 +B22) D7= (A11 −A21)(B11 +B12). This algorithm, with its O(nlg 7)time complexity, could speed up the design on large matrices. We should note however, that the speed of the extra additions have to be carefully considered. Since the multiplication is already extremely fast, a similar improvement may also be necessary for additions if the overall performance upgrade is to remain significant. Acknowledgements Research supported by the T´ AMOP 4.2.1/B-09/1/KONV-2010-0007 project and TARIPAR3 project grant Nr. TECH 08-A2/2-2008-0086. References [1] F. Bensaali, A. Amira, R. Sotudeh, Floating-point matrix product on FPGA, IEEE/ACS International Conference on Computer Systems and Applications, Amman, Jordan, 2007, pp. 466–473. ⇒173 [2] D. Coppersmith, S. Winograd, Matrix multiplication via arithmetic progressions, J. Symbolic Comput. 9, 3 (1990) 251–280. ⇒173 [3] N. Dave,K.Fleming, M. King, M. Pellauer, M. Vijayaraghavan, Hardware acceleration of matrix multiplication on a Xilinx FPGA, MEMOCODE ’07 Proc. 5th IEEE/ACM International Conference on Formal Methods and Models for Co-Design, Nice, France, 2007, pp. 97–100. ⇒173 Modular exponentiation of matrices on FPGA-s 191 [4] Y. Dou, S. Vassiliadis, G. K. Kuzmanov, G. N. Gaydadjiev, 64-bit floating-point FPGA matrix multiplication, Proc. 2005 ACM/SIGDA 13th international symposium on Field-programmable gate arrays, Monterey, CA, USA, 2005, pp. 86–95. ⇒173 [5] T. J. Earnest, Strassen’s Algorithm on the Cell Broadband Engine, 2008, http://mc2.umbc.edu/docs/earnest.pdf ⇒173 [6] T. Herendi, Uniform distribution of linear recurrences modulo prime powers, J. Finite Fields Appl. 10, 1 (2004) 1–23. ⇒174 [7] T. Herendi, Construction of uniformly distributed linear recurring sequences modulo powers of 2 (to appear). ⇒174 [8] A. Irturk, S. Mirzaei, R. Kastner, An Efficient FPGA Implementation of Scalable Matrix Inversion Core using QR Decomposition, UCSD Technical Report, CS2009-0938, 2009. ⇒173 [9] B. Kakaradov, Ultra-fast matrix multiplication, An empirical analysis of highly optimized vector algorithms, Stanford Undergraduate Research Journal 3(2004) 33–36. ⇒173 [10] M. Karkooti, J. R. Cavallaro, C. Dick, FPGA implementation of matrix inversion using QRD-RLS algorithm, Proc. 39th Asilomar Conference on Signals, Systems, and Computers, Pacific Grove, CA, USA, 2005, pp. 1625–1629. ⇒173 [11] S. M. Qasim, A. A. Telba, A. Y. AlMazroo, FPGA design and implementation of matrix multiplier architectures for image and signal processing applications, IJCSNS International Journal of Computer Science and Network Security 10, 2 (2010) 168–176. ⇒173 [12] V. Strassen, Gaussian elimination is not optimal, Numer. Math. 13 (1969) 354–356. ⇒173 Received: May 17, 2011 •Revised: October 11, 2011