We propose to store several integers modulo a small prime into a single
machine word. Modular addition is performed by addition and possibly
subtraction of a word containing several times the modulo. Modular
Multiplication is not directly accessible but modular dot product can be
performed by an integer multiplication by the reverse integer. Modular
multiplication by a word containing a single residue is a also possible.
Therefore matrix multiplication can be performed on such a compressed storage.
We here give bounds on the sizes of primes and matrices for which such a
compression is possible. We also explicit the details of the required
compressed arithmetic routines.Comment: Published in: MICA'2008 : Milestones in Computer Algebra, Tobago :
  Trinit\'e-et-Tobago (2008

Dumas, Jean-Guillaume

Fousse, Laurent

Salvy, Bruno

English

arXiv

International audienceWe propose to store several integers modulo a small prime into a single machine word. Modular addition is performed by addition and possibly subtraction of a word containing several times the modulo. Modular Multiplication is not directly accessible but modular dot product can be performed by an integer multiplication by the reverse integer. Modular multiplication by a word containing a single residue is a also possible. Therefore matrix multiplication can be performed on such a compressed storage. We here give bounds on the sizes of primes and matrices for which such a compression is possible. We also explicit the details of the required compressed arithmetic routines

Hal - Université Grenoble Alpes

Compressed Modular Matrix Multiplication

INRIA a CCSD electronic archive server

HAL Id: hal-00259950https://hal.archives-ouvertes.fr/hal-00259950v2Submitted on 13 Mar 2008HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestinée au dépôt et à la diffusion de documentsscientifiques de niveau recherche, publiés ou non,émanant des établissements d’enseignement et derecherche français ou étrangers, des laboratoirespublics ou privés.Compressed Modular Matrix MultiplicationJean-Guillaume Dumas, Laurent Fousse, Bruno SalvyTo cite this version:Jean-Guillaume Dumas, Laurent Fousse, Bruno Salvy. Compressed Modular Matrix Multiplication.MICA’2008 - Milestones in Computer Algebra, May 2008, Tobago, Trinidad and Tobago. pp.133-140.￿hal-00259950v2￿hal-00259950, version 2 - 13 Mar 2008Compressed Modular Matrix MultiplicationJean-Guillaume Dumas∗ Laurent Fousse∗ Bruno Salvy†March 13, 2008AbstractWe propose to store several integers modulo a small prime into a sin-gle machine word. Modular addition is performed by addition and possi-bly subtraction of a word containing several times the modulo. ModularMultiplication is not directly accessible but modular dot product can beperformed by an integer multiplication by the reverse integer. Modularmultiplication by a word containing a single residue is a also possible.Therefore matrix multiplication can be performed on such a compressedstorage. We here give bounds on the sizes of primes and matrices forwhich such a compression is possible. We also explicit the details of therequired compressed arithmetic routines.1 IntroductionCompression of matrices over fields of characteristic 2 is naturally made via thebinary representation of machine integers [1, 8].The FFLAS/FFPACK project has demonstrated the need of a wrapping ofcache-aware routines for efficient small finite field linear algebra [4, 5].Therefore, a conversion between a modular representation of prime fields ofany (small) characteristic and e.g. floating points can be performed via thehomomorphism to the integers [2]. In [3] it is proposed to transform polynomialover a prime field into a Q-adic representation where Q is an integer thanthe field characteristic. We call this transformation DQT for Discrete Q-adicTransform. With some care, in particular on the size of Q, it is possible to mapthe polynomial operations into the floating point arithmetic realization of thisQ-adic representation and convert back using an inverse DQT.Efficient matrix computations over very small finite fields of characteristic otherthan two are required e.g. to study strongly regular graphs [9], in order toprove/disprove and help in the comprehension of the conjectures of [10].In this note we propose to use this fast polynomial arithmetic within machinewords to compute dot products. We show in section 2 how to recover a dot∗Laboratoire J. Kuntzmann, Université de Grenoble, umr CNRS 5224. BP 53X, 51, rue desMathématiques, F38041 Grenoble, France. {Jean-Guillaume.Dumas,Laurent.Fousse}@imag.fr†Projet ALGO, INRIA Rocquencourt, 78153 Le Chesnay. France. Bruno.Salvy@inria.fr1product of size d+1 can be recovered from the single coefficient of degree d of apolynomial product. Whenever the prime modulus is small enough this enablesto compute several accumulations of binary products in a single machine oper-ation. Then we propose in section 3 an alternative matrix multiplication usingmultiplication of a compressed word by a single residue. The latter requiresalso a simultaneous modular reduction, called REDQ in [3]. In general, theprime field, the size of matrices and the available mantissa are given. This givessome constraints on the possible choices of Q and d. In both cases anyway, weshow that these compression techniques represent a speed-up factor of up to thenumber d + 1 of residues stored in the compressed format.2 Q-adic compression or Dot product via poly-nomial multiplicationSuppose that a(X) =∑di=0 aiXi and b(X) =∑di=0 biXi are two polynomialsin Z/pZ[X ]. One can perform the dot product∑di=0 aibd−i by extracting thecoefficient of degree d of a(X)b(X).2.1 Modular dot product via machine word multiplicationThe idea here, as in [3], is to replace X by an integer Q, usually a power of 2in order to speed up conversions. Thus the vectors of residues a = [a0 . . . ad]and b = [b0 . . . bd] are stored respectively as b̄ =∑di=0 biQi and the reverseā =∑di=0 ad−iQi.This is done e.g. over floating points via the following compressions:double& init3( double& r,const double u, const double v, const double w) {r=u; r*=_dBase; r+=v; r*=_dBase; return r+=w;}2.2 GainNow for matrix multiplication A × B one wishes to convert a whole row of theleft m× k matrix A and a whole column of the right k×n matrix B. Thus A istransformed into a m×⌈kd+1⌉CompressedRowMatrix, CA and B is transformedinto a⌈kd+1⌉× n CompressedRowMatrix, CB.Therefore the matrix multiply CA × CB can gain a factor of d + 1 over themultiplication of A × B, for classical multiplication as shown on the 2 × 2 ex-ample below where the matrix product[a bc d]×[e fg h]=[ae + bg af + bhce + dg cf + dh]isperformed via integer multiplications. The precise gain will be given in table 2.2[Qa + bQc + d]× [e + Qg f + Qh] =[∗ + (ae + bg)Q + ∗.Q2 ∗ + (af + bh)Q + ∗.Q2∗ + (ce + dg)Q + ∗.Q2 ∗ + (cf + dh)Q + ∗.Q2]The result matrix C = CA × CB is m × n. Thus in order to compare similarcomputations one has either to consider multiplication of compressed matriceswhich is then the procedureC = CA × CB; CC = ReduceAndCompress(C) (1)or to consider multiplication of normal matrices via compression and thus theprocedureCA = CompressRows(A); CB = CompressColumns(B); C = CA × CB (2)2.3 Partial compressionNote that the last column of CA and the last row of B might not have d + 1elements if kd+1 /∈ Z. Thus one has to artificially append some zeroes to theconverted values. On b̄ this means just do nothing. On the reversed ā thismeans multiplying by Q several times.2.4 Delayed reduction and lower bound on QFor the results to be correct the inner dot product must not exceed Q. With apositive modular representation mod p (i.e. integers from 0 to p−1), this meansthat (d + 1)(p − 1)2 < Q. Moreover, we would like to use delayed reductionson the intermediate results and thus accumulate the āb̄ before any modularreduction. It is thus possible to perform matrix multiplications of with commondimension k as long as:kd + 1(d + 1)(p − 1)2 = k(p − 1)2 < Q. (3)2.5 Available mantissa and upper bound on QIf the product āb̄ is performed with floating point arithmetic we just need thatthe coefficient of degree d remains fully in the mantissa β. Write āb̄ = cHQd+cL,the latter means that cH , and cH only, must remain lower that 2β. It couldthen be exactly recovered by multiplication of āb̄ by the correctly precomputedand rounded inverse of Qd and floored, as shown e.g. in [3, Lemma 2].With delayed reduction this means that∑di=0kd+1(i+1)(p− 1)2Qd−i < 2β . Wecan use equation 3 in order to show that∑di=0kd+1(i + 1)(p− 1)2Qd−i ≤ Qd+1.With this we just have to enforce thatQd+1 < 2β. (4)Thus a single reduction has to be made at the end of the dot product as follows:3Element& init( Element& rem, const double dp) const {double r = dp;// Multiply by the inverse of Q^d with correct roundingr *= _inverseQto_d;// Now we just need the part less than Q=2^tunsigned long rl( static_cast<unsigned long>(r) );rl &= _QMINUSONE;// And we finally perform a single modular reductionrl %= _modulus;return rem = static_cast<Element>(rl);}Note that one can avoid the multiplication by the inverse of Q when Q = 2t:by adding Q2d+1 to the final result one is guaranteed that the t(d+1) high bitsrepresent exactly the d+1 high coefficients. On the one hand, the floating pointmultiplication can be replaced by an addition. On the other hand, this doublesthe size of the dot product and thus reduces by a factor of d+1√2 the largestpossible dot product size k.2.6 ResultsOne can see on figure 1 that the compression (d + 1) is very useful for smallprimes since the gain over the double floating point routine is quite close to d.Indeed choosing a power of 2 for Q simplifies and speeds up conversions andthus gives the following compression factors modulo 3:Compression 2 3..4 5..8 8 7 6 5 4 3Degree d 1 9 7 6 5 4 3 2Q-adic 23 24 25 26 27 28 210 213 217Dimensions 2 ≤ 4 ≤ 8 ≤ 16 ≤ 32 ≤ 64 ≤ 256 ≤ 2048 ≤ 32768Table 1: Compression factors for different common matrix dimensions modulo 3,with 53 bits of mantissa and Q a power of 2.Before n = 256 the compression is at a factor of five and the time to perform amatrix multiplication is less than a hundredth of a second. Then from 257 to2048 one has a factor of 4 and the times are roughly 16 times the time of thefour times smaller matrix, as visually shown on figure 2, left. The same is trueafterwards with the respective factor 3 of compression.Remark that the curve of fgemm with underlying arithmetic on single floatsoscillates and drops. This is because the matrix begins to be too large and thatmodular reductions are required between the recursive matrix multiplicationsteps. Then the floating point BLAS1 routines are used only when the sub-matrices are small enough. One can see the subsequent increase in the numberof classical arithmetic steps on the drops around 2048, 4096 and 8192.1http://www.tacc.utexas.edu/resources/software/4 0 5 10 15 20 25 0  2000  4000  6000  8000  10000GIGA finite field operations per secondMatrix orderFinite field Winograd matrix multiplication with Goto BLAS on a XEON, 3.6 GHzdouble fgemm mod 11float fgemm mod 3Compressed double fgemm mod 3dgemmFigure 1: Compressed matrices multiplication of equation 1 compared withdgemm (the floating point double precision routine of GotoBlas) and fgemm(the exact routine of FFLAS) with double or single precision.3 Right Compressed matrix multiplicationAnother way of performing compressed matrix multiplication is to multiply anuncompressed matrix m × k to the right by a row-compressed k × nd+1 ma-trix. A dot product with this algorithm will be of the form a = [a0, . . . , an] ×[∑dj=0 b0jQj , . . . ,∑dj=0 bnjQj ]. Therefore, a single entry of the resulting matrixwill be∑ki=0 ai(∑dj=0 bijQj) =∑dj=0(∑ki=0 aibij)Qj as shown on the examplebelow.[a bc d]×[e + Qfg + Qh]=[(ae + bg) + Q(af + bh)(ce + dg) + Q(cf + dh)]Here also Q and d must satisfy equations (3) and (4).The major difference is in the reductions. Indeed now one needs to reducesimultaneously the d + 1 coefficients of the polynomial in Q in order to get theresults. This simultaneous reduction can be made by the REDQ algorithm of[3, Algorithm 2].Thus the whole right compressed matrix multiplication over two compressedmatrices CA and CB, is the following algorithm as shown also on figure 2,right:A = Uncompress(CA); CC = A × CB; REDQ(CC) (5)54 Full compressionOf course one would like to compress simultaneously two dimensions of thematrix product. The required dot products can there be achieved by polynomialmultiplication with two variables Q and Θ. Let dq be the degree in Q and dθbe the degree in Θ. The dot product is then: a = [∑dqi=0 ai0, . . . ,∑dqi=0 ain] ×[∑dθj=0 b0jΘj , . . . ,∑dθj=0 bnjΘj ]. The latter is∑kl=0(∑dqi=0 ail)(∑dθj=0 blj)QiΘj =∑dqi=0∑dθj=0(∑kl=0 ailblj)QiΘj as shown on the example below.[a + Qc b + Qd]×[e + Θfg + Θh]= [(ae + bg) + Q(ce + dg) + Θ(af + bh) + QΘ(cf + dh)]In order to guarantee that all the coefficients can be recovered independently,Q must still satisfy equation (3) but then we have this additional equation onΘ:Qdq+1 ≤ Θ (6)This gives thus upper bounds on dq and dθ:Q(dq+1)(dθ+1) < 2β (7)5 ComparisonWe summarize the differences of the presented algorithms on figure 2 and 3.C = A BUncompressedxA, Row CompressedB, Column CompressedxB, Row CompressedC = A B, Row CompressedA, UncompressedFigure 2: Algorithms of equations (1), left, and (5), right.We see on the one hand that the first algorithm compresses the common di-mension whereas the Right (or also Left) compressions compress an external6xC=AB, Column CompressedA, Column CompressedB UncompressedB, Row CompressedxA, Column Compressed C=ABFigure 3: Left Compression and Full Compressionmatrix dimension. Thus in the case of rectangular matrices one can choose be-tween those routines the fastest one. This will be the routine compressing thelargest dimension as shown on table 2. On the other hand the full compressionalgorithm compresses both external sizes but, as shown by equation (7) theavailable mantissa is is only shared by both compression. Therefore if we definethe compression factor to bee =⌊βlog2(Q)⌋then the degree of compression for the first three algorithms is just be d = e− 1where it becomes d =√e − 1 for the full compression with equal degrees forboth variables Q and Θ. For ω the exponent of matrix multiplication, the table2 shows that the gain in terms of arithmetic operations is eω−2 for the first threevariants and eω−12 for the full compression. When ω = 3 for classical matrixmultiplication the speed-up is the same. Now, when fast matrix multiplicationis used, as e.g. in [6, §3.2], full compression performs less operations. This isnot only of theoretical interest but also of practical value since the consideredmatrices are then less rectangular. This enables more locality for the matrixcomputations and usually better performance.Algorithm Operations Reductions Conversions(1) O(mn(ke)ω−2)m × n REDC 1emn INITe(5) O(mk(ne)ω−2)m × neREDQe1emn EXTRACTeLeft Comp. O(nk(me)ω−2)me× n REDQe1emn EXTRACTeFull Comp. O(k(mne)ω−12)m√e× n√eREDQe1emn INITeTable 2: Number of operations for the different algorithmsThe difference there will mainly be on the number and on the kind of modularreductions. Since the REDQe reduction is faster than e classical reductions,7see [3], and since INITe and EXTRACTe are roughly the same operations, thebest algorithm would then be one of the Left, Right or Full compression.For example, with algorithm (1) on matrices of sizes 10000×10000 it took 92.75seconds to perform the matrix multiplication modulo 3 and 0.25 seconds toconvert the resulting C matrix. This is less than 0.3%. For 250 × 250 matricesit takes less than 0.0028 seconds to perform the multiplication and roughly0.00008 seconds for the conversions. There, the conversions count for 3%.The full compression algorithm seems the better candidate for the locality andfast matrix multiplication reasons above ; howbeit the compression factor is aninteger, depending on the flooring of either βlog2(Q)or√βlog2(Q). Thus there arematrix dimensions for which the compression factor of e.g. the right compressionwill be larger than the square of the compression factor of the full compression.There the right compression will have some advantage over the full compression.Further work would thus include implementing the Right or Full compressionand comparing the conversions overhead with that of algorithm (1).References[1] Don Coppersmith. Solving linear equations over GF (2): block Lanczosalgorithm. Linear Algebra and its Applications, 192:33–60, October 1993.[2] Jean-Guillaume Dumas. Efficient dot product over finite fields. In Vic-tor G. Ganzha, Ernst W. Mayr, and Evgenii V. Vorozhtsov, editors, Pro-ceedings of the seventh International Workshop on Computer Algebra inScientific Computing, Yalta, Ukraine, pages 139–154. Technische Univer-sität München, Germany, July 2004.[3] Jean-Guillaume Dumas. Q-adic transform revisited. Technical Re-port 0710.0510 [cs.SC], ArXiv, October 2007. http://hal.archives-ouvertes.fr/hal-00173894.[4] Jean-Guillaume Dumas, Thierry Gautier, and Clément Pernet. Finite fieldlinear algebra subroutines. In Teo Mora, editor, Proceedings of the 2002International Symposium on Symbolic and Algebraic Computation, Lille,France, pages 63–74. ACM Press, New York, July 2002.[5] Jean-Guillaume Dumas, Pascal Giorgi, and Clément Pernet. FFPACK: Fi-nite field linear algebra package. In Jaime Gutierrez, editor, Proceedingsof the 2004 International Symposium on Symbolic and Algebraic Computa-tion, Santander, Spain, pages 119–126. ACM Press, New York, July 2004.[6] Jean-Guillaume Dumas, Pascal Giorgi, and Clément Pernet. Dense linearalgebra over prime fields. ACM Transactions on Mathematical Software,2008. to appear.[7] Chao H. Huang and Fred J. Taylor. A memory compression scheme formodular arithmetic. IEEE Transactions on Acoustics, Speech, and SignalProcessing, 27(6):608–611, December 1979.8[8] Erich Kaltofen and Austin Lobo. Distributed matrix-free solution of largesparse linear systems over finite fields. In A.M. Tentner, editor, Proceedingsof High Performance Computing 1996, San Diego, California. Society forComputer Simulation, Simulation Councils, Inc., April 1996.[9] John P. May, David Saunders, and Zhendong Wan. Efficient matrix rankcomputation with application to the study of strongly regular graphs. InChristopher W. Brown, editor, Proceedings of the 2007 International Sym-posium on Symbolic and Algebraic Computation, Waterloo, Canada, pages277–284. ACM Press, New York, July 29 – August 1.[10] Guobiao Weng, Weisheng Qiu, Zeying Wang, and Qing Xiang. Pseudo-paley graphs and skew hadamard difference sets from presemifields. De-signs, Codes and Cryptography, 44(1-3):49–62, 2007.9

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.218.6439

Compressed Modular Matrix Multiplication

Abstract

Similar works

Full text

Available Versions

Hal - Université Grenoble Alpes

INRIA a CCSD electronic archive server