{'fraction_non_alphanumeric': 0.031419982822788434, 'fraction_numerical': 0.016425708559977097, 'mean_word_length': 4.3968713789107765, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Compressed Matrix Computations

Frugal computing is becoming an important topic for environmental reasons. In this context, several techniques have been proposed to reduce the storage of scientific data by dedicated compression methods specially tailored for arrays of floating-point numbers. While these techniques are quite efficient to save memory, they introduce additional computations to compress and decompress the data before processing them. In this article, we introduce a new lossy, fixed-rate compression technique for 2D-arrays of floating-point numbers which allows one to compute directly on the compressed data, without decompressing them. We obtain important speedups since less operations are needed to compute among the compressed data and since no decompression and re-compression is needed. More precisely, our technique makes it possible to perform basic linear algebra operations such as addition, multiplication by a constant among compressed matrices and dot product and matrix multiplication among partly uncompressed matrices. This work has been implemented into a tool named blaz and a comparison with the well-known compressor zfp in terms of execution-time and accuracy is presented.


I. INTRODUCTION
At the beginning of this new decade, more than ten percent of the world total electricity production is consumed by information and communication technology systems [14]. This huge consumption can be roughly decomposed into three equivalent parts: storage (30%), computing (30%) and networking (40%). This distribution imposes to who wants to propose frugal computing solutions, for ecological or maybe economical reasons, to attack these three points all together: reducing storage, computing and networking simultaneously. In this article, we propose a solution to this problem for the special case of matrix computations [6] which is ubiquitous in computer science.
Clearly, at the beginning of this new decade, scientific data remain the most important part of the world stored data, before ERP or entertainment data for example [8]. Then it urgent to seek ways of processing scientific data which offer simultaneously gains in terms of storage, computing and networking. In this article, the scientific data that we consider are large matrices of floating-point numbers [1] used to represent smooth curves (as in other inspiring work [4], [18], we will assume that there exist correlations between adjacent elements of matrices).
Because usual compression techniques [9], [11] are not efficient enough on scientific data, specific (lossy) compressors have been developed such as zfp [18] or sz [4]. The performance of these tools is impressive in terms of data size reduction and accuracy, allowing to reduce storage by 10 to 20 without significant loss of relevant information. Doing so, they help to reduce the footprint of scientific data storage and communications. However, they increase the computing resources needed to solve a problem, since data must be uncompressed before processing and re-compressed after.
In this article, we introduce a new lossy compression technique for scientific data which makes it possible to compute directly among the compressed data, without decompression. While our method is not as efficient as zfp or sz in terms of compression rate and accuracy, it allows one to compute directly with the compressed data. This avoids a huge amount of operations since i) data do not need to be uncompressed and re-compressed and ii) the compressed data being smaller, less operations are needed to compute with them. More precisely, our technique operates upon two-dimensional arrays storing matrices and allows one to perform basic linear algebra operations such as addition and multiplication by a constant on the compressed data and dot product and matrix multiplication on partly decompressed matrices. To our knowledge, this is the first attempt to define a compression technique enabling one to compute directly on the compressed data.
We can give a bird's eye view of our technique as follows. Our compressor is lossy, fixed-rate (11.37 compression rate) and operates on 8×8 blocks of IEEE754 binary64 floatingpoint numbers [1]. After a block-splitting stage, we perform successively stages of normalization, prediction, transformation and, finally, quantization, all described in details in Section II.
A key-point is that all these stages act like linear maps (at the notable exception of the quantization done at the end of the scheme), i.e. the transformation t applied at any stage (but quantization) of our scheme satisfies t(A + B) = t(A) + t(B) and t(cA) = ct(A) where A and B are matrices represented in the adequate data structures and where c is a scalar constant. These properties are used in Section III to design algorithms to add compressed matrices and to multiply them by constants without decompression. Additionally, we will also see in Section III that the dot product among lines and columns of compressed matrices as well as the product A×B of matrices can be computed on partially uncompressed matrices.
Last but not least, our compression scheme has been implemented into an open-source library named blaz and experimental results are given in Section IV. They show that the loss of accuracy introduced by our technique remains small, compared to the gains obtained in terms of execution time and storage. For instance, the addition of matrices of size 2000×2000 with blaz is more than 50 times faster than with standard, uncompressed matrices and more than 5000 times faster than with zfp (see Section IV-A). On the other hand, The relative errors introduced by blaz are less than 10 times greater than with zfp for our basis of examples (see Section IV-B). This corresponds to the loss of one more decimal digit than with zfp. We strongly believe that this compromise is quite acceptable in many contexts.
As already mentioned, this article is organized as follows. The compression scheme used by blaz is detailed in Section II and the algorithms used to compute among the compressed structures are introduced in Section III. The experimental results are given in Section IV and Section VI concludes.

II. COMPRESSION SCHEME
In this section, we introduce the five stages of our compression scheme which perform block-splitting, normalization, prediction, transformation and quantization. The originality of our compression scheme is to be compatible with the basic linear algebra operations that can be performed directly on the compressed data structures using the algorithms introduced in Section III. A summary of this scheme is given in Figure 1 and we detail the different steps hereafter. a) Block Splitting: First of all, as in most fixed-rate compressors [18], we split the original matrix into blocks. In practice, our blocks have size 8 × 8.
Example II.1 All along this section, we are going to illustrate how our compressor works using the following 8 × 8 matrix We may assume that A is one of the blocks obtained by splitting a larger original matrix. The values of A correspond to the functions f (x, y) = x×y displayed in Figure 7 for x and y starting at 0 with a step of 0.1. The graphical representation of Block A is also given in Figure 3. b) Block Normalization: After block splitting, the next step of our scheme is to normalize the values. This stage is twofold. First, for each block element, we compute the difference between itself and its preceding elements. More precisely, taking as entry a 8 × 8 matrix M corresponding to a block, we compute the new block ∆ such that, for any Normalization relies on the assumption that block elements are correlated and form a smooth surface. The differences between consecutive values are then assumed to be small and this step aims at reducing the range of values occurring in the block. In addition to the matrix ∆, we have to store the value of the first element of the block M 00 into a binary64 number.
Example II.2 By normalizing the matrix A of Example II.1 using the formula displayed in Equation (1), we obtain the new matrix Let us observe that the values in A N are smaller than those in A. We also store the element A 00 = 0.0 into a separate binary64 number.
The second step of the normalization stage consists of dividing the values in ∆ by the mean slope between consecutive values of the block A (considering only the non zero elements.) Since ∆ already contains the differences between adjacent elements of A, we have and we compute Note that the mean slope s must be stored into a binary64 number by our compressor.  Let s max be the maximal value in ∆ in absolute value, as our slopes of reference (see Figure  2). Then each value in ∆ is approximated by its closest element in S. Assuming that, for all 0 ≤ k < 255, s k ≤ s k+1 , we store in a new matrix Ψ the indexes k of the slopes s k which are the best approximations of the values in ∆ . Formally, for 0 ≤ i, j < 8, we have Note that Ψ is a 8 × 8 matrix of integers encoded on 8 bits.
In addition to the mean slope s, we need to store s max into a binary64 floating-point number.
The binary64 values A 00 = 0.0, and s = 0.04 are also stored. d) Block Transform and Quantization: As in many other compressors (e.g. JPEG-2000 [9]), we use a Type II twodimensional discrete cosine transform (DCT, [22]) of the block resulting from the normalization stage. DCTs are used to aggregate large coefficients in the first lines and columns of a matrix, small values occurring in the other elements after transformation. The quantization then consists of considering that these small coefficients are equal to 0 and of avoiding to store them in the compressed matrix. Our compressor only keeps the elements of the first two lines and columns of the matrix, i.e. 28 values. Note that the values returned by the DCT may be larger than 127 in absolute value. Let m be the greatest value in absolute value of the matrix D resulting from the transform. We re-scale D by multiplying all its elements by In order to recover the original matrix during the decompression phase, the value ϕ must also be stored into a 8 bits integer.
Example II.5 After re-scaling, the DCT of matrix P of Example II.4 gives the new matrix For this example, all the non-zero coefficient but one in position (4, 4) already occur in the first two lines and columns of the matrix P (in bold in the equation above) and the quantization only sets P 44 to 0 introducing a slight information loss. Recall that, in addition to P , the values A 00 = 0.0, s = 0.04 and the scale factor ϕ = 20 are also stored.
Example II. 6 The difference between the original block of Example II.1 and the result of its compression and decompression by our method is displayed in Figure 3.
Following our scheme, a compressed block is stored using 29 8-bits integers and 2 binary64 floating-point numbers. The 29 8-bits integers correspond to the 28 values that we keep after quantization plus one 8-bits integer corresponding to the re-scaling factor ϕ. A block is then stored into 360 bits instead of the original 4096 bits needed to store 64 binary64 numbers, yielding a compression rate of 11.37.

III. BASIC LINEAR ALGEBRA ALGORITHMS
In this section, we introduce the basic linear algebra algorithms used to compute among matrices following the scheme of Section II. More precisely, in what follows, we are going to introduce successively the algorithms for addition, multiplication by a constant, dot product and matrix multiplication. These operations are detailed for elementary 8 × 8 blocks. Obviously, for larger matrices, this has to be repeated for each elementary block.
All along this section, we use the following notations. Let B denote a 8 × 8 block of a compressed matrix, f and s denote respectively the values of the first element B 00 and of the mean slope s of B as defined in Equation (2). Both f and s are binary64 numbers. Next, ϕ denotes the scale factor ϕ introduced in Equation (6) Equation (7) motivates the computation carried out at Line 2 of Algorithm 1. Let D 1 and D 2 denoted the blocks obtained after the DCT by passing B 1 and B 2 through the scheme of Figure 1. Intuitively, if no re-scaling were done in our scheme, we could simply add D 1 and D 2 to obtain the block D corresponding to B. But two re-scaling are done in our scheme, during the normalization stage and among the values resulting from the DCT. The coefficients α 1 and α 2 are used to transform blocks made for scale factors ϕ 1 and ϕ 2 to the new scale factor ϕ, accordingly to Equation (4). Similarly, β 1 and β 2 are used to adapt the scale of the blocks resulting from the DCT to the scale factor ϕ. It is then possible to re-scale and add the coefficients of the DCT contained in C 1 and C 2 as done at Line 6, in the for loop of Algorithm 1. Note that adding the coefficients is possible since the DCT defines a linear map among our blocks. Formally, for a 8 × 8 block Ψ as defined in Equation (5), the two-dimensional DCT of Ψ, denoted DCT(Ψ) = D is defined, for 0 ≤ i, j < 8, by Ψij cos (2u + 1)iπ 16 cos (2v + 1)jπ 16 (8) and it follows that for two blocks Ψ = Ψ + Ψ , Ψij cos where the coefficients α i are the usual DCT coefficients defined by Example III.1 Figure ?? displays the points of the result of the addition of two 8×8 blocks corresponding to the functions x × y and x 2 × y 2 for x and y starting at (0, 0) with a step of 0.1. We display the result of the addition between the uncompressed blocks and the result of the addition of the compressed blocks (this latter result being uncompressed after the operation).
Let us remark that the subtraction of blocks works similarly to the addition. We omit to detail it in the present article.
b) Multiplication by a Constant: In our framework, the multiplication by a constant c is the simplest operation to implement. The first element of the block and the slope are multiplied by c. The other elements remain unchanged. This is summarized in Algorithm 2. This algorithm is straightforward and we do not give more details about it. c) Dot product: Let B 1 and B 2 be two 8 × 8 matrix blocks and let B 1 , B 2 ij denote the dot product between the i th line of B 1 and the j th column of B 2 , 0 ≤ i, j < 8. Our dot product requires a partial decompression of the blocks. More precisely, we have to recover the values resulting from the prediction stage of our scheme (third step of Figure 1 and Section II-c) and this implies to compute the inverse discrete cosine transform (IDCT) of the values contained in C 1 and C 2 (obviously, the values discarded by the quantization are replaced by zeros.) Let P 1 = IDCT(C 1 ) and P 2 = IDCT(C 2 ) be the blocks obtained by applying the inverse discrete cosine transform to C 1 and C 2 . We have to compute Algorithm 1 Addition of two compressed matrix blocks B 1 and B 2 .
Note that we do not need to compute the entire blocks B 1 and B 2 . It is sufficient to compute the first i lines of B 1 and the first j columns of B 2 in order to compute the dot product B 1 , B 2 ij . Algorithm 3 summarizes how we perform the dot product following equations (11) and (12). In this algorithm, the function IPREDICT LINES(C, i) and IPREDICT COLS(C, j) use Equation (12) to compute respectively the first i lines and j columns of B 1 and B 2 . d) Matrix multiplication: Our algorithm for the multiplication B = B 1 ×B 2 of two blocks uses the same ideas than for the dot product and we do not detail it hereafter. The difference comes from the fact that the blocks resulting from the IDCT and from Equation (12) are needed many times since, for each 0 ≤ i, j ≤ 7, the element B ij corresponds to the dot product B ij = B 1 , B 2 ij . To avoid redundant computations, these blocks are computed only once.
For large matrices made of many 8 × 8 blocks, the algorithms introduced in this section need to generalized. Briefly speaking, this corresponds to apply block-wise the former algorithms.

IV. EXPERIMENTAL EVALUATION
Our compression scheme as well as the operations among compressed matrices have been implemented into an opensource library named blaz 1 , written in C, and, in this section, we present two kinds of experiments, concerning the time performances (Section IV-A) and the accuracy (IV-B) of blaz.
Altogether these experiments show that our technique makes it possible to save simultaneously storage and computation time. In both cases, comparisons with the zfp [18] library are reported. All these experiments have been carried out on a quad core Intel Core i7-1165G7 processor with 16 gigabytes of RAM and running Ubuntu 20.04.3 LTS. For the sake of simplicity, all our experiments are carried out on square matrices even though our techniques works for matrices of 1 https://github.com/mmartel66/blaz any dimensions. For the relevance of the comparison, the compression rate of zfp is set to a ratio equivalent to blaz.

A. Performance
In this section, the performances of the basic linear algebra operations performed by blaz are compared to the same operations done without any compression and to the case where matrices are compressed by zfp [18] (C++ version).
First, we compare the execution time of operations among matrices compressed with blaz to the execution time of the same operations performed among uncompressed matrices. We use square matrices of sizes ranging from 80 to 2000 and the measures are done for addition, multiplication by a constant and for the dot product. Our results are displayed in Figure 4. Let us mention that the time in the y-axis of Figure 4 is given in logarithmic scale and that the execution time for the uncompressed dot product is almost null and is not represented. The first noticeable point is that the addition and multiplication by a constant are significantly faster with our compressed matrices than without any compression. Indeed, for square matrices of size 2000, the execution time is divided by 62 for addition (0.15ms vs 9.35ms) and by 58 for the multiplication by a constant (0.14ms vs 8.14ms.) Our technique makes it possible to compute faster than without any compression when small losses of accuracy are acceptable (see Section IV-B.) Concerning the dot product, the compressed operation is slower with a compressed matrix than without compression (0.12ms vs 0.004ms which is 30 times slower.) This is due to the fact that a partial decompression of the matrix is done in this case (see Algorithm 3).
To complete our study, we have also measured the time needed by blaz to compress and decompress matrices. These times are given in Figure 5.
Second, we compare blaz to the well-known compressor zfp. For both tools, we start with compressed matrices. In the     case of zfp, when a matrix element is read, the corresponding block is uncompressed and the value is returned. Conversely, when some value is assigned to a certain element, zfp directly modifies and re-compress the corresponding block. In the case of blaz, no decompression/re-compression is needed, the operations are carried out directly on the compressed matrices. The execution times are given in Figure 6 for our basic linear algebra operations. As in Figure 4, the execution time represented on the y-axis is displayed in logarithmic scale.
Again, we take square matrices of size ranging from 80 to 2000 and we show on a same graph the execution time taken by blaz, zfp and without compression. The top two graphs of Figure 6 are for addition and multiplication by a constant. The main observation is that zfp introduces a huge overhead due to the decompression and re-compression of matrices. For addition, taking as reference the execution time without any compression, blaz and zfp introduce respectively a speedup of 3 and a speed down of 121 for matrices of size 80 and a speedup of 62 and a speed down 115 for matrices of size 2000. In other terms, the addition is more than 50 times faster with blaz and more than 100 times slower with zfp compared to the uncompressed operation. Note that similar results are observed for the multiplication by a constant.
The execution times for the dot product are displayed in the left bottom corner of Figure 6. As already mentioned, this operation carried out with blaz matrices is slower than without any compression. Nevertheless, it remains significantly faster than with zfp matrices. Again, taking as reference the execution time without compression, blaz is 31 times slower for matrices of size 2000 while zfp is 62 times slower. In this case, blaz remains 2 times faster than zfp.
Finally, the graph in the bottom right corner of Figure  6 show the performances of matrix multiplication. A first observation is that blaz is not much slower than the uncompressed operation. This is due to the fact that time spent in the matrix multiplication dominates the partial compression/decompression time. The second observation is that zfp is much slower than blaz (approximately 10 times.)

B. Accuracy Measurements
In this section, we introduce a second set of experiments concerning the accuracy of blaz. Again, we compare our tool to zfp. Contrarily to the execution time measurements introduced in Section IV-A, the values of the matrix elements now matter. We use matrices corresponding to the discretisation of the six non-linear functions f k : R 2 → R, 1 ≤ k ≤ 6, presented in Figure 7. We set M kij = f k (x j , y i ) where the points (x j , y i ) are taken in the range [−2, 2] × [−2, 2] with a constant step depending on the size of M k .
Our experimental results are displayed in Table I. They correspond to mean relative errors computed as follows. Let M , M 1 and M 2 denote square matrices of size n whose elements are floating-point numbers. Let C and U denote the compression and decompression functions and let M = U(C(M )), the mean relative error e on M is We proceed similarly for matrices resulting from basic linear algebra operations. For instance, we apply Equation (13) to M = M 1 + M 2 and M = U(C(M 1 ) + C(M 2 )).
The first two lines of Table I indicate the mean relative error on the compression/decompression of the matrices M 1 to M 6 using blaz and zfp. While zfp is more accurate than blaz, this latter remains accurate with relative errors of less than 0.5% in half of the cases and never greater than 2%.
The second part of Table I is dedicated to addition. We display the mean relative errors for all the additions M i + M j , 1 ≤ i, j ≤ 6, i = j performed on blaz matrices (upper right triangle) and on zfp matrices (lower left triangle). For example, the mean relative errors on M 1 + M 6 with blaz and zfp respectively are 0.78% and 0.31%. The errors introduced by blaz and zfp are comparable (same magnitude in general) even if zfp is better in most cases (at the price of the huge overhead in terms of execution time shown in Section IV-A.) In some cases, blaz is more accurate than zfp (e.g. for M 4 + M 6 blaz error is 1.70% while zfp error is 2.32%.) We may conclude that, for addition, there is not a clear advantage to use zfp instead of blaz.
The third part of Table I deals with the multiplication by a constant and, in our experiments, the matrices M 1 to M 6 are multiplied by 2. While being less accurate than zfp, the accuracy of the computations remain accurate, with relative errors around 1%. Note that changing the constant does not change the relative error.
The last part of Table I is for matrix multiplication. Let us mention that, in blaz, the value and accuracy of the dot product M, M ij is the same than these of the element (i, j) of the matrix M × M . Then we only discuss the accuracy of the matrix multiplication which encompasses the one of the dot product. The matrix product not being symmetric in general, we display in Table I the mean accuracies of all the products M i × M j , 1 ≤ i, j ≤ 6, computed following Equation (13). These accuracies are all for blaz. The errors obtained with zfp being almost zero in any case, we do not display them. We can see that the accuracy of blaz for matrix multiplication is very good (relative errors of order 10 −3 or less in more than half of the cases and never greater than 10 −1 in magnitude), even if they remain worse than zfp. We can also remark that blaz multiplication is more accurate than blaz addition. This is due to the fact that a partial decompression is done for the multiplication, contrarily to the case of addition. This better accuracy then results from more computations and, in future work, we would like to design even more frugal algorithms (yet less accurate) for the dot product and multiplication.

V. RELATED WORK
Roughly speaking, compressors for scientific data (i.e. arrays of floating-point numbers) can be classified into two categories: Either the user sets the compression rate and this determines the accuracy of the encoding or the user sets the accuracy and this determines the compression rate. The compressors zfp [18] and fpzip [19] before it have been initially designed as compressors of the first category. Conversely, sz [4], [17] has been designed as a compressor of the second category. Elements of comparison between these two approaches can be found in [23]. Note that the current versions of zfp and sz offer many options which make them encompass the categories mentioned earlier. Let us also mention an extension of sz [26] which uses spline functions for more accuracy and another approach based on topological control [21].
Another related topic is about the use of GPUs for scientific data compression [13], [16], [24], [25]. Finally, much attention has also been paid recently to the compression of neural networks [12], [15], [20]. We strongly believe that the ideas developed in this article could be useful in the context of neural networks which perform mainly basic linear algebra operations.

VI. PERSPECTIVES AND CONCLUSION
The work presented in this article opens a new research direction and many perspectives remain to be explored. First of all, we would like to determine formal error bounds for our compressor (in the spirit of [5]) and, more interestingly, for the operations on compressed matrices. We strongly believe that the error on the result of some operation between compressed matrices can be expressed in function of the errors on the compressed matrices corresponding to the operands.
Second, we aim at making the compression rate tunable, more or less compression being permitted with the corresponding impact on accuracy. This can be achieved by skewing the quantization or by introduction an interpolation stage in our scheme. For example, we could add a Lorenzo predictor which would be compatible with our basic linear algebra operations [10].
Third, we believe that compressors allowing to compute on the compressed matrices could be useful in other contexts. For example a compressor for integers could be useful to perform fast and frugal basic linear algebra operations in the context of image processing. As well, a blaz-like compressor designed for embedded systems could be useful for many applications with memory constraints.
Next, we plan to extend our library, with new linear algebra operators but also, our compression scheme processing every 8×8 block of a matrix independently, we would like to develop a parallel implementation of our blaz library, typically based on the MPI library [7].
Yet another perspective is to use compressed matrices in precision tuning tools [3]. Precision tuning consists of reducing the precision of the variables of a program up to some accuracy requirement set by the user. It would be interesting to use the precision recommended by the precision tuning tools to compress arrays up to this requirement.
Finally, we aim at testing further the blaz library. We aim at evaluating it in the context of large numerical simulations to assess its performances in term of execution-time and accuracy on more realistic use-cases [2], [27].