This is the accepted version of the article: Blanes, Ian; Hernández-Cabronero, Miguel; Serra-Sagristà, Joan; [et al.]. «Re- dundancy and optimization of tANS entropy encoders». IEEE transactions on multimedia, 2021. DOI 10.1109/TMM.2020.3040547 This version is available at https://ddd.uab.cat/record/235940 under the terms of the license 1 Redundancy and Optimization of tANS Entropy Encoders Ian Blanes, Senior Member, IEEE, Miguel Hernández-Cabronero, Joan Serra-Sagristà, Senior Member, IEEE, and Michael W. Marcellin, Fellow, IEEE Abstract—Nowadays entropy encoders are part of almost all image compression [23], compressive sensing [24], or point data compression methods, with the Asymmetrical Numeral cloud data compression [25]. In [26], ANS is employed to Systems (ANS) family of entropy encoders having recently risen code large-alphabets sources. in popularity. Entropy encoders based on the tabled variant of ANS are known to provide varying performances depending on While most of the interest in ANS from the data compression their internal design. In this paper, we present a method that community has been focused on practical efforts, other rele- calculates encoder redundancies in almost linear time, which translates in practice to thousand-fold speedups in redundancy vant papers have been published recently, such as an approx- calculations for small automatons, and allows redundancy cal- imate method of calculating the efficiency of encoders [27], culations for automatons with tens of millions of states that [28] an improved probability approximation strategy [29], or would be otherwise prohibitive. We also address the problem their use as extremely parallel decoders [30]. of improving tabled ANS encoder designs, by employing the aforementioned redundancy calculation method in conjunction In this manuscript, we study Tabled ANS (tANS), which is with a stochastic hill climbing strategy. The proposed approach an ANS variant that is well suited for hardware implemen- consistently outperforms state-of-the-art methods in tabled ANS tations [31], [32], [33]. It employs a finite state machine encoder design. For automatons of twice the alphabet size, experimental results show redundancy reductions around 10% that transitions from state to state as symbols are encoded, over the default initialization method and over 30% for random which is certainly not a new idea [34]. However, the ANS initialization. theory enables the creation of effective encoding tables for such encoders. The underlying idea of this ANS variant is that Index Terms—Tabled Asymmetrical Numeral Systems, Entropy for non-integer amounts of information an encoder produces encoder redundancy, Optimization. integer-length codewords, while state transitions are employed to carry over fractional bits of information to subsequent I. I NTRODUCTION codewords. O F the diverse set of tools that are employed in data compression, entropy encoders are one of the most used tools, if not the most. They are the last stage in almost all In particular, we focus our study on efficiently measuring the redundancy produced by tANS encoders, and on optimizing tANS encoders to minimize their implementation footprint data compression methods, after a varying set of predictors or to improve their coding efficiency. We introduce a novel and transforms. One particular family of entropy encoders that redundancy measurement strategy that exploits the structure has gained popularity in recent years is the one based on the imposed by tANS encoders in state transition matrices. This Asymmetrical Numeral Systems (ANS) [1], [2], [3]. Some strategy achieves almost-linear time complexities for redun- ANS encoders and decoders have been shown to allow for dancy calculations when combined with an efficient calcula- very fast software implementations in modern CPUs [4], [5]. tion of average codeword length for all states. In addition, This has lead to its incorporation in recent data compression we employ the proposed redundancy calculation technique to standards and its use in many different cases [6], [7], [8], optimize tANS automatons. We do so by exploring automaton [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19]. In permutations through a stochastic hill climbing method. We addition, ANS-based encoders could be applicable to a very expect these efforts to translate into, for example, less area be- wide range of multimedia scenarios, such as an alternative ing allocated to a hardware implementation, to having smaller to the Rice-Golomb codes employed in the energy-efficient coding tables sent as side information, or to compression meth- scheme described in [20], as a high-throughput entropy en- ods that can dynamically determine adequate encoder sizes. coder in a high frame rate video format [21], or in general as Many compression techniques incorporate entropy encoders an entropy encoder in schemes for sparse coding [22], learned with suitably high efficiency. In this work, we aim at providing better tANS encoders at equivalent efficiency levels, but with This work was supported in part by the Centre National d’Etudes Spatiales, and by the Spanish Ministry of Economy and Competitiveness (MINECO), lower hardware costs. Particularly relevant to this research is by the European Regional Development Fund (FEDER) and by the European the optimization strategy for tANS encoders in [35]. Further Union under grant RTI2018-095287-B-I00, and by the Catalan Government comparison with this method is provided in the following under grants Beatriu de Pinós BP2018-00008 and 2017SGR-463. This paper has supplementary downloadable material available at sections. http://ieeexplore.ieee.org., provided by the author. The material contains the test frequency tables referenced in this document. This material is 304 KB in The rest of this paper is organized as follows. In what size." remains of this introduction we present some background 2 information on tANS. Afterward, we provide an efficient them as least-significant bits to x0 so that the resulting x is in method to calculate the redundancy of tANS encoders in I. Section II, and we show a method to optimize tANS encoders in Section III. In Section IV we discuss our experimental Algorithms 1 and 2 below describe the procedure to encode results, and conclusions are drawn in Section V. and decode one symbol, respectively. The encoding and de- coding functions are uniquely determined by the key f as discussed in more detail in what follows. A. Tabled Asymmetric Numeral Systems Algorithm 1 Algorithm for encoding one symbol. In this section we provide the necessary background on how Input: s, x, s, ms , Cs tANS encoders operate and establish some common notation. Output: y, s For the underlying rationale behind tANS see [2]. See Table V x0 ← x for a description of each symbol. while x0 > 2ms − 1 do Let an n-symbol discrete memoryless source have an under- push(x0 mod 2, s) lying alphabet A = {0, . . . , n − 1}. For symbol s ∈ A, let ps x0 ← bx0 /2c be the associated occurrence probability, which we assume to end while be finitely representable (otherwise, see [29]). y ← Cs (x0 ) A tANS encoder (or decoder) of size m is a deterministic finite-state automaton that transitions from state to state as Algorithm 2 Algorithm for decoding one symbol. symbols are being encoded (or decoded). Such automaton is Input: y, s, m, D defined by its key vector, f = f0 f1 . . . fm−1 ∈ Am , which Output: s, x, s uniquely describes the encoding and decoding functions that (s, x0 ) ← D(y) the automaton employs (i.e., how are symbols mapped to while x0 < m do codewords and which state transitions occur).PLet ms be the x0 ← 2x0 + pop(s) number of occurrences of symbol s in f , with s∈A ms = m. end while For example, given A = {0, 1, 2}, the key f = 00210 defines x ← x0 an automaton of size m = 5, with m0 = 3 and m1 = m2 = 1. Given the current state of the automaton, represented as an As an example, the encoding and decoding functions for the integer x ∈ I = {m, . . . , 2m − 1}, the symbol s is encoded automaton with key f = 10211011 are provided in Table I. as follows: For this key, m0 = 2, m1 = 5, m2 = 1 and m = 8. To encode symbol s = 0 in this example, suppose that the current state 1) First, state x is renormalized to x0 ∈ Is = of the automaton is x = 8. Renormalizing x into x0 ∈ I0 = {ms , . . . , 2ms −1} by removing as many least significant {m0 , . . . , 2m0 − 1} = {2, 3} requires taking the two least bits from x as necessary (i.e., a bitwise right shift significant bits from x, resulting in x0 = 2. Applying C0 to operation), and pushing them to a stack s (least-significant x0 yields y = 9, which is the new state of the automaton. The bit first). two least significant bits “00” are pushed to the stack. 2) Then, an encoding function Cs : Is → I is employed to The decoder starts with state y = 9 and applies D to obtain produce the resulting state Cs (x0 ) = y. s = 0 and x0 = 2. Renormalizing x0 by popping bits from Hence, the encoding process results in a state transition from the stack and appending until x ∈ I = {m, . . . , 2m − 1} = x ∈ I to y ∈ I and some bits pushed onto a stack. {8, . . . , 15} yields x = 8. The resulting compressed file contains only the contents of A key f unambiguously defines functions Cs and D as follows: the stack, the final state of the automaton, and, if it cannot • The first element of the ordered pair produced by the be deduced, the number of symbols encoded. The initial state decoding function is obtained directly from the symbols, of the automaton need not be included in the compressed file, in order, contained in f . That is, and can be chosen arbitrarily. However an initial state x = m is expected to produce less renormalization bits for the first D(x) = (fx−m , x0 ), m ≤ x < 2m. (1) encoded symbol (renormalizing any initial state larger than m The second element of the ordered pair, x0 , is given after produces at least as many least-significant bits and possibly the definition for Cs below. more). • The values in the coding table for Cs are, in order, the Decoding operations are performed in reverse order, with the states that have the symbol s as the first element of the peculiarity that the sequence of decoded symbols is produced ordered pair in the coding table for D. in reversed order. To decode a symbol, a decoding function D : I → (A, Is ) is employed to produce a decoded symbol s • Destination states in D (the second element of the ordered and a state x0 ∈ Is . Then, state x0 is renormalized by popping pair) are inverse to those in Cs . I.e., D(x) = (s, x0 ) where as many bits as necessary from the stack s and appending x = Cs (x0 ). 3 Table I: Encoding and decoding tables for the automaton is a computationally expensive procedure, which prevents its with key ‘10211011’. use in important applications. In particular, straightforward x C0 C1 C2 D calculation can become infeasible for large automatons, em- pirical redundancy studies, and most importantly, data-driven 1 10 key optimization procedures. 2 9 It is worth noting that a redundancy approximation to tANS 3 13 was provided by Duda [2] and later refined by Yokoo [28]. 4 However, the approximations may present notable divergences 5 8 from the true redundancy or may fail to distinguish the best 6 11 of several tANS encoders. See Fig. 1 for an example where Duda’s redundancy approximation yields the same result for 7 12 two different automatons, and where redundancy is signifi- 8 14 (1, 5) cantly underestimated around p1 ' 0.2. As seen in the figure, 9 15 (0, 2) the same statement is true for Yokoo’s approximation. 10 (2, 1) 100 11 (1, 6) 12 (1, 7) 13 (0, 3) 10−1 Redundancy (bits) 14 (1, 8) 15 (1, 9) 10−2 f1 = 10000 It is well known that low encoding redundancies are obtained 10−3 f2 = 00001 Duda’s approximation for f1 and f2 for keys where ms ≈ m ps [2]. However, as others have Yokoo’s approximation for f1 and f2 shown [35] and we further show, lower redundancies can be 10−4 achieved by taking into account the order of symbols in f . An 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 efective key construction method is described in [2], which we p1 use as a baseline. We reproduce that method in Algorithm 3 Figure 1: Example where Duda’s and Yokoo’s redundancy with slight modifications to ensure that all symbols appear approximations yield identical results for two different at least once in the key. In the algorithm, heap tuples are automatons having f1 and f2 , which in reality have compared lexicographically, and the ‘pop’ operation returns significantly different redundancy profiles. Additionally, the smallest element. both methods dramatically underestimate redundancy for Algorithm 3 Algorithm for baseline key creation. p1 ' 0.2. Input: m, ps ∀s ∈ A Output: f  In this section we describe a procedure to efficiently calculate g ← Set 0, . . . , n − 1 average codeword lengths and, in turn, the redundancy of h ← Heap (0.5/p0 , 0), . . . , (0.5/pn−1 , n − 1) tANS encoders. for i = 0 to m − 1 do repeat The average codeword length of a tANS encoder can be (v, s) ← pop(h) written as X until s ∈ g or |g| < m − i L= P(x) · `x (2) x∈I push((v + 1/ps , s), h) remove(s, g) where P(x) is the probability of the automaton being in the fi ← s state x in a given instant, and `x is the per-state average end for codeword length. State transitions can be modeled as a Markov process for II. E FFICIENT R EDUNDANCY C ALCULATION which P(x) is the stationary probability of being in state x. These probabilities can be obtained as the dominant unitary As for any entropy encoder, the redundancy of a tANS left eigenvector of the transition matrix P (i.e., the dominant encoder for a given source can be obtained as the difference unitary right eigenvector of PT ). Here we consider only between the entropy of the source and the average codeword irreducible aperiodic Markov chains, so that unique stationary length produced by the encoder. Thus, less redundant encoders probabilities are guaranteed to exist. In particular, stationary produce smaller compressed files. probabilities for automatons with reducible Markov chains While redundancy calculation for tANS is well understood, may not be unique, i.e., there can be more than one stationary the straightforward calculation of average codeword lengths distribution for a given automaton. On the other hand, a 4 stationary probability may fail to exist for a periodic Markov be obtained from (5) directly. This method is formalized in chain. While outside the scope of this article, if necessary, an Algorithm 4. approximate solution can be found for these cases by adding Regarding the complexity of the algorithm, given that |A| ≤ a small constant to all the elements of the stochastic matrix |I|, it can be seen that complexity is dominated by the first and thus ensuring irreducibly and aperiodicity [36]. and last loops. Each loop performs m individual operations, Regarding per-state average codeword lengths, these can be and thus the algorithm complexity is O(m). obtained as Algorithm 4 Per-state average codeword length calculation. X `x = ps C(s, b x), (3) s∈A Input: m, ms , ps ∀s ∈ A j k Output: `x ∀x ∈ I b x) = lg x is the number of bits pushed to where C(s, ms for x ∈ I do the stack by Algorithm 1 when encoding symbol s in state x. `0x ← 0 Here, and in what follows, note that lg denotes the base-two end for logarithm. It is well known that, for each symbol s, C(s, b x) for s ∈ A do can only take two values, which are consecutive integers, and k ← blg(m/ms )c + 1 that a threshold ts on x suffices to distinguish the correct t ← ms 2k value [2], [35]. Specifically, `0t ← `0t + ps `0m ← `0m + ps (k − 1) ( ks − 1 x < ts , C(s, x) = b (4) end for ks x ≥ ts , `m ← `0m where ks = blg(m/ms )c + 1 and ts = ms 2ks . for x = m + 1 to 2m − 1 do `x ← `x−1 + `0x . In what follows we show how to obtain `x and P (x) in a end for computationally efficient manner, to facilitate the calculation of Eq. 2 (as illustrated in Fig. 2). First we describe how to calculate per-state average codeword lengths directly from the B. Efficient calculation of the state transition matrix and its automaton size, symbol occurrence counts in the automaton dominant eigenvector key, and symbol probabilities. Afterwards, we describe how Having already seen how to obtain `x , in this section we to obtain state probabilities by exploiting the structure of the proceed to obtain P(x) by exploiting the particular structure state transition matrix of a tANS encoder. We obtain a compact of the transition matrix P of tANS automatons. representation this matrix, which we then use in a modified The re-normalization process in tANS creates structure and power method to obtain its dominant eigenvector. regularity in the automaton transition matrix which we can exploit. Due to re-normalization it can easily be seen that, A. Per-state average codeword length in Algorithm 1, multiple and consecutive input values of x produce the same output value (same state transition) due to While a straightforward approach to calculate all `x values the floor operation applied in the re-normalization step. through (3) and (4) requires O(mn) operations, we show a method that obtains per-state average codeword lengths in For example, we show in Fig. 3 the transition matrix P and O(m) by employing finite differences. Note that throughout key f for a three-symbol automaton assuming p0 = 0.2, the document we assume that arithmetic operations have a p1 = 0.3, and p2 = 0.5. It is readily apparent that row i of O(1) complexity, including exponentiation and logarithms. PT is formed by runs of value pfi , which occasionally wrap around. I.e., runs containing the probability of the symbol in Applying a forward difference operator (i.e., ∆n an = an+1 − the same row of f T , which occasionally reach the right-most an ) to `x yields column continue on the left-most column. We can employ this fact to efficiently create equivalent compact representations X ∆x `x = ps ∆x C(s, b x), (5) s∈A of the transition matrices of tANS automatons. In addition, it is possible to observe that runs for the same symbol are where ( successive and of equal length, up to the wrapping point, where 1 for x = ms 2ks , ∆x C(s, b x) = (6) its length is halved, but we do not exploit this for our purposes. 0 otherwise. We can represent row y−m of PT , associated with destination Given symbol s, the value of ∆x C(s, b x) is only 1 for a state y ∈ I as My = (α, β, p), where the interval [α, β) single x value, which implies that each term in the summation specifies a run of origin states x ∈ I that lead to y with in (5) is only non-zero for a single x value. Given this fact, probability p. For runs that wrap around, the end of the run we can obtain all values of ∆x `x at once by iterating over β is increased by m so that α < β, which simplifies notation alphabet symbols and only calculating non-zero terms in the in further steps. In the example of Fig. 3, there are 15 states summation. Values of `x can then be found through cumulative labeled 15 through 29. We then have M15 = (20, 22, 0.2) summation from ∆x `x values, except for `m , which needs to for the top row (row 0) of PT , and M27 = (28, 31, 0.3) 5 Figure 2: Diagram of the redundancy calculation method. 0     0.2 0.2  0.2 0.2   0    0.2 0.2     0     0.2 0.2     0     0.3 0.3 0.3 0.3     1    0.3 0.3 0.3 0.3   1  0      0.2 0.2   T P = 0.2 ,  fT =  0     0.2     0     0.2     0     0.2     0     0.3 0.3 0.3 0.3     1   0.3 0.3 0.3  1  0     0.2 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 2 Figure 3: Transition matrix (transposed) and key for an automaton with n = 3, p0 = 0.2, p1 = 0.3, p2 = 0.5, and m = 15. for the third row from the bottom (row 12). We can find Now we see that we can employ the power method on a all states x ∈ I with destination state y ∈ I starting from compact transition matrix to produce its dominant eigenvector the re-normalized state x0 ∈ Is such that Cs (x0 ) = y as (i.e., the vector with elements P(x) ∀x ∈ I). It is well known follows. We look for all x ∈ I that re-normalize to x0 , i.e., that the power method can be employed by iterating over x = x0 2k + z, with 0 ≤ z ≤ 2k − 1 and k, z ∈ Z. For PT vi T z = 0 we obtain a run start by solving x0 2k + 0 ≥ m, which vi+1 T = (7) ||PT vi T || yields k = dlg m/x0 e, as only one value for k is possible [2]. Similarly, for z = 2k − 1 we obtain a run end by solving   to produce a dominant eigenvector in O − lgR|λ2 | iterations, x0 2k + (2k − 1) + 1 ≥ m, which yields k = dlg m/(x0 + 1)e. 0 0 where R is the required precision and λ2 is the second- Hence, α = x0 2dlg m/x e , and either β = (x0 +1)2dlg m/(x +1)e 0 largest eigenvalue [37] (see [38] for further estimation of this or β = (x0 + 1)2dlg m/(x +1)e + m, depending whether there eigenvalue). In what follows, we show that, by exploiting the is a wrap around or not. This yields My = (α, β, p). transition matrix  structure,we can obtain the dominant eigen- Hereafter we refer to this representation, M , as a compact vector in a O m − lgR|λ2 | complexity. To do so, we show transition matrix, whichSAlgorithm 5 obtains in O(m) time how to obtain the result of the matrix-vector multiplication complexity, given that | s∈A Is | = m. wT = PT vT in the power method in O(m) instead of O(m2 ). Algorithm 5 Compact transition matrix creation. For non-wrapping runs, the ith element of wT can be formu- Input: ps , ms , Cs ∀s ∈ A lated as Output: M X wi = pfi vj−m , (8) for s ∈ A do α≤j<β for x0 ∈ Is do 0 α ← x0 · 2dlg(m/x )e with α and β being the closed and open boundaries of the run β ← (x 0 ( + 1) · 2 dlg(m/(x0 +1))e in row i of the transition matrix, as described previously. For β if α < β wrapping runs we can consider β← β + m otherwise   X X MCs (x) = (α, β, ps ) wi = pfi  vj−m + vj−m  . (9) end for α≤j<2m m≤j<β−m end for Due to the particular definition of β, both summations in (9) 6 are equivalent to 3.0   m = 25 m = 100 2.5 X m = 1000 wi = pfi  uj−m  , (10) m = 10000 Relative redundancy (%) α≤j<β m = 100000 2.0 where u = (v0 , . . . , vm−1 , v0 , . . . , vm−1 ). Thanks to this, arbi- trary summations of continuous ui elements can be obtained 1.5 in O(1) by subtracting two elements from a pre-calculated cumulative sum of the elements in u. 1.0 The resulting method to obtain stationary state probabilities for 0.5 tANS automatons is presented in Algorithm 6. As previously discussed, Algorithm 6 converges in O − lgR|λ2 | iterations, 0.0 0.0 0.2 0.4 0.6 0.8 1.0 with a O(m) cost per iteration.  Thus the complexity of Algo- p0 rithm 6 is O m − lgR|λ2 | . This complexity dominates those from the previous steps and results in the total complexity for Figure 4: Redundancy achieved with automatons of a given the tANS redundancy calculation. size where keys are restricted to be in order.1 Algorithm 6 Compact power method. Input: m, M call ordered automatons. There are m+1 such automatons. The Output: v 26 dotted black curves in the figure correspond to all ordered 1 v ← ||1|| automatons with m = 25. The blue curve is the minimum of repeat 2 these. The other colored curves show the minimum curve for v0 ← v other values of m. For comparison, the black curve represents u00 ← 0 the redundancy obtained by non-ordered automatons produced for i ← 1 to 2m − 1 do by Algorithm 3 for m = 100. It can be seen that, as m in- u0i ← u0i−1 + v(i−1 mod m)+1 creases, the redundancies obtained by even very large ordered end for automatons seem to stagnate, still having relative redundancies for i ← 0 to m − 1 do significantly larger than the automatons with non-ordered (α, β, p) ← Mi+m symbols. Note how the curves for m = 1000, m = 10000, and vi ← p · (u0β − u0α ); m = 100000 are nearly identical. This might have interesting end for implications for the redundancies achievable by the rANS v ← ||v||v variant of ANS, which for redundancy purposes are equivalent 2 until ||v − v0 ||∞ < 2−R to ordered tANS automatons, as it seems to suggest that rANS encoders may not converge to zero redundancy as automatons grow, or that they do so at a much slower rate than non-ordered At this point we have presented a method to obtain the average tANS automatons. codeword length of a tANS automaton and thus its redundancy in a very efficient manner, which allows us to explore the In the remainder of this section we describe how to obtain optimization of automatons. well-performing orders for the elements of tANS automaton keys. The approach, which is based on hill climbing, is as fol- III. R EDUNDANCY REDUCTION lows: given an initial automaton key, a succession of pseudo- The basic principle behind ANS dictates that, when encoding random pairwise permutations of key elements are considered. a symbol, an encoder should transition to a next state that is Whenever an element permutation reduces the automaton approximately 1/ps times larger than the current state (i.e., redundancy, the permutation is applied to the key, otherwise y ≈ x/ps ), in order to accommodate the extra lg(1/ps ) bits the permutation is discarded. Algorithm 3 can be employed of information [2]. When this is translated to tANS, it suggests to create initial automaton keys that have appropriate symbol that the number of occurrences of each symbol in tANS frequencies. keys should be ms ≈ m ps . However, in addition to symbol Traditionally, measuring the effects of any replacement or frequencies, the order of these symbols in the key also plays permutation of the key elements has been computationally an important role in tANS performance, as larger states emit challenging. In particular, obtaining the dominant eigenvector more symbols to reach the same Is during re-normalization, after a change in a state transition matrix is a difficult and well- and not all automaton states are equally probable [2]. studied problem [37], [39]. Existing solutions are significantly For example, Fig. 4 shows the redundancy of all two-symbol more computationally expensive than the straightforward ap- automatons that have keys sorted in ascending order, which we plication of our redundancy calculation method, or are just 1 This figure required more than 50 million redundancy evaluations for not very effective. For example, employing an approximate automatons of 105 states. Calculation was possible thanks to the redundancy eigenvector as the starting point for the power method may calculation method presented in the previous section. seem promising, but is not as effective as intuition would 7 suggest [37]. Hence, in what follows, we directly employ the proposed redundancy calculation method to evaluate the effect 20 of a key change. Redundancy reduction (%) Informed decisions could also be employed to explore the 15 space of automaton keys. For example, the estimator E = −P(x)`x − P(y)`y + P(x)`y + P(y)`x 10 (11) = (P(x) − P(y))(`y − `x ) 5 could be employed to guide the selection of keys to be evaluated. This estimator is the redundancy variation after 0 the elements of a key associated with states x and y are 101 102 103 104 105 106 107 swapped, under the assumption that state probabilities will Iterations remain unaltered. Recently, Dubé and Yokoo proposed a similar approach in [35], where key elements are rearranged proba02 book1_bzip2 aviris_123 as per the probability of their associated automaton state in descending order. In their state-of-the-art proposal this Figure 5: Redundancy reduction versus number of itera- approach is repeated multiple times, and the key producing tions for the proposed method (solid curves) and for the the least redundancy is finally selected. method by Dubé and Yokoo (dashed curves). Results are for the first frequency table of each dataset. The motivation of such approach is to reduce complexity by reducing the number of keys to be evaluated. Given the low cost of evaluating keys via the proposed method, such representative of the expected input to the entropy encoder approaches may not be necessary. Indeed, directly exploring of a general-purpose compression method. Three additional pseudo-random permutations reaches more key variations than synthetic distributions are included in the corpus. those suggested by heuristics and approximations, and ulti- mately achieve lower redundancies. Regarding the performance of the method presented in Sec- Fig. 5 shows redundancy reductions versus iterations for tion II, in addition to the previous theoretical time complexity three different frequency tables for our proposed method as discussion, several wall-clock measurements are reported here. compared to the informed approach proposed by Dubé and Table III reports on performance of the proposed redundancy Yokoo. It is worth emphasizing that for both approaches calculation method. Results are produced on an AMD Thread- the plots indicate actual redundancy reduction calculated via ripper 1950X with 3200 Mhz DDR4 RAM, by employing algorithms proposed earlier in this paper, without employing a hybrid implementation of higher level Python 3 code and any approximations or heuristics. In each case, the key size lower-level elements either in NumPy or custom C code. m is set to approximately 5n. It can be observed that the All results are for single-thread execution (including BLAS informed approach yields significant redundancy reductions routines). We have also developed a ‘reference’ implemen- with very few iterations, with little or no further improvement tation for comparison purposes where the power method provided by subsequent iterations. In contrast, for pseudo- is employed, but none of the other complexity reduction random permutations, significant redundancy reductions are strategies described in Section II are used. As compared with only obtained after a moderate number of iterations, but a this reference implementation, the proposed method yields larger redundancy reduction is eventually obtained. substantial execution time improvements. Albeit the reference implementation is not as optimized as the proposed one, the difference in wall-clock performance is of around three or four IV. E XPERIMENTAL R ESULTS orders of magnitude. For m ≥ 105 , results for the reference In this section we first present the corpus employed to test the method have not been produced due to its quadratic memory described methods, followed by experimental results regard- requirements. In addition, the proposed technique has been ing the performance of the method proposed in Section II, checked against this reference implementation for correctness, and results regarding the optimization method described in with additional checks against LAPACK routines for the power Section III. method. In order for the experimental results to be applicable in Regarding the performance of the optimization method pre- practical scenarios, a data corpus has been carefully curated. sented in Section III, results are presented in Fig. 6 and Ta- The corpus, detailed in Table II, consists of multiple frequency ble IV. Results in Fig. 6 show the relation between automaton tables obtained in various real data compression experiments.2 size and its redundancy for two different frequency tables. Two of the corpus datasets (aviris_123 and deer_iwt) represent Results for a ‘default’ automaton (Algorithm 3) are compared what an entropy encoder may expect during image compres- with the proposed ‘optimized’ method, and with three random sion, whereas another two (book1_bzip2, enwik8_gzip) are permutations of the ‘default’ automaton. It can be observed 2 The corpus is available online as supplementary material through IEEE that random permutations tend to yield poor results, whereas Xplore. the optimized method provides improvements over the default 8 Table II: Details of the corpus of frequency tables employed in the experimental results. The minimum, average, and maximum number of symbols is reported for each set of frequency tables. Name Tables Min. n Avg. n Max. n Source aviris_123 11 355 1,272.1 4,891 CCSDS 123.0-B-2 applied to the AVIRIS Yel- lowstone hyperspectral image (scene 00). book1_bzip2 24 92 92.0 92 Bzip2 (-9) on the ‘BOOK1’ file from the Calgary corpus. deer_iwt 16 3,412 9,395.2 14,447 Subbands of a 5-level 5/3 LeGall IWT transform applied to the ‘deer’ image from the 16-bit (log) Rawzor corpus. enwik8_gzip 1,985 13 82.1 219 Gzip (-9) on M. Mahoney’s 100 MB English Wikipedia test file. proba02 1 256 256.0 256 Synthetic distribution (02) from Y. Collet’s ‘probaGenerator’. proba14 1 53 53.0 53 Synthetic distribution (14) from Y. Collet’s ‘probaGenerator’. proba80 1 7 7.0 7 Synthetic distribution (80) from Y. Collet’s ‘probaGenerator’. Table III: Wall-clock measurements of the proposed method (measured in seconds). The first table of each dataset is employed for the measurements. Invalid automaton sizes (m < n) are denoted by ‘-’. Reference Proposal Speedup m book1_bzip2 deer_iwt proba02 book1_bzip2 deer_iwt proba02 book1_bzip2 deer_iwt proba02 10 2 0.152 - - 0.001 - - 152.0 - - 103 1.723 - 0.496 0.002 - 0.002 861.5 - 248.0 104 32.940 116.316 5.479 0.011 0.014 0.010 2994.5 8308.3 547.9 105 0.130 0.155 0.099 106 Insufficient memory 3.614 9.409 1.891 Insufficient memory 107 39.588 191.075 23.173 automaton. V. C ONCLUSIONS This paper addresses the problem of efficiently obtaining tANS For some frequency tables, such as that in Fig. 6(a), larger automaton redundancies and the optimization of said automa- automatons consistently yield smaller redundancies. However, tons. We describe a method to efficiently obtain automaton for the frequency table presented in Fig. 6(b) and others, the redundancies, and we do so by efficiently obtaining per- redundancy obtained by default automatons may significantly state average codeword lengths through a forward difference increase as their size is increased. Thus, even for cases operator, and by exploiting the structure of the state transition where automaton optimization may not be worth the extra matrix to obtain the stationary probabilities of the automa- computational effort, carefully evaluating automaton size can ton states. Our proposed method yields important benefits.   operates in almost-linear time complexity O m − lg |λ2 | , and three to four orders of R magnitude faster in practice. In addition, we show that we In Table IV, comprehensive results are presented for all can employ the proposed redundancy calculation method in frequency tables of all data sets. It is worth mentioning a stochastic hill climbing optimization process to improve that automaton keys that are close in size to the num- tANS encoder designs. Experimental results with novel corpus ber of alphabet symbols may present less opportunities to of frequency tables obtained from realistic data compression improve performance by permuting key elements (e.g., the processes indicate that automatons are consistently improved improvement obtained for the proba02 table is negligible for over other informed strategies of automaton optimization. In m ' 1.1n). However, these keys may also fail to mimic the particular, results show that improvements increase as automa- frequency table. This is particularly relevant for the highly- tons grow larger. biased table in the proba80 data set, where an automaton with m ' 1.1n = 7.7 cannot effectively represent the required occurrence probabilities. However, as automaton sizes are VI. ACKNOWLEDGMENTS increased, redundancy reductions become substantial, and even more so as compared to randomly selected permutations of an We would like to thank Charles Bloom for his well-thought- automaton key. out “rants” blog. 9 Table IV: Redundancy reductions, in percentage, over default and random key construction methods. Automaton sizes are set as a multiple of the alphabet size. Results are averaged over all frequency tables in each data set. For each result 50 thousand iterations are used. m ' 1.1n m ' 2n m ' 5n Corpus Default Random Default Random Default Random aviris_123 0.35 17.05 5.08 29.43 7.63 31.87 book1_bzip2 0.52 23.70 10.11 35.17 11.60 39.12 deer_iwt 2.17 41.33 17.40 72.66 17.04 90.83 enwik8_gzip 5.42 49.32 28.96 76.26 37.04 89.59 proba02 0.00 22.59 10.96 38.88 21.80 66.41 proba14 2.79 26.56 9.19 40.12 11.45 51.11 proba80 0.00 0.00 0.00 45.17 5.36 28.01 Table V: Symbol list. Symbol Description A An alphabet α, β Beginning and end of a run in My Cs Coding function for symbol s C b Bits pushed into the stack (function) D Decoding function ∆n Forward difference operator f Automaton key fi i’th element in f g A set (a) First table in ‘enwik8_gzip’. h A heap I Interval {m, . . . , 2m − 1} Is Interval {ms , . . . , 2ms − 1} ks Threshold value in C b for symbol s L Automaton average codeword length `x Average codeword length for state x λ2 Second largest eigenvalue M Compact transition matrix My Run of origin states that lead to state y in a compact transition matrix m Automaton size ms Number of occurrences of s in f (b) First table in ‘deer_iwt’. n Alphabet size Figure 6: Redundancy as a function of automaton size. P Probability transition matrix P(x) Probability of automaton being in state x A PPENDIX ps Occurrence probability of symbol s R Precision value A. Notation s A symbol in A A summary of the symbols employed through this manuscript s A stack is included in Table V. u, v, w Vectors employed in power method ui , vi , wi i’th elements in u, v, w, respectively R EFERENCES x, x , y 0 Automaton states [1] J. Duda, “Asymmetric numeral systems,” arXiv preprint arXiv:0902.0271, 2009. 10 [2] ——, “Asymmetric numeral systems: entropy coding combining speed [22] M. Kalluri, M. Jiang, N. Ling, J. Zheng, and P. Zhang, “Adaptive RD of Huffman coding with compression rate of arithmetic coding,” arXiv Optimal Sparse Coding With Quantization for Image Compression,” preprint arXiv:1311.2540, 2013. IEEE Trans. Multimedia, vol. 21, no. 1, pp. 39–50, Jan. 2019. [3] J. Duda, K. Tahboub, N. J. Gadgil, and E. J. Delp, “The use of [23] Z. Cheng, H. Sun, M. Takeuchi, and J. Katto, “Energy compaction- asymmetric numeral systems as an accurate replacement for Huffman based image compression using convolutional autoencoder,” IEEE Trans. coding,” in 2015 Picture Coding Symposium (PCS), May 2015, pp. 65– Multimedia, pp. 1–14, Aug. 2019 (early access). 69. [24] X. Yuan and R. Haimi-Cohen, “Image Compression Based on Com- [4] Y. Collet, “Finite State Entropy,” http://github.com/Cyan4973/ pressive Sensing: End-to-End Comparison with JPEG,” IEEE Trans. FiniteStateEntropy, 2013. Multimedia, pp. 1–16, Jan. 2020 (early access). [5] F. Giesen, “Interleaved Entropy Coders,” arXiv preprint [25] P. de Oliveira Rente, C. Brites, J. Ascenso, and F. Pereira, “Graph-Based arXiv:1402.3392, 2014. Static 3D Point Clouds Geometry Coding,” IEEE Trans. Multimedia, vol. 21, no. 2, pp. 284–299, Feb. 2019. [6] Y. Collet and C. Turner, “Smaller and faster data compression with Zstandard,” Facebook Code, 2016. [On- [26] A. Moffat and M. Petri, “Large-Alphabet Semi-Static Entropy Coding line]. Available: https://code.facebook.com/posts/1658392934479273/ Via Asymmetric Numeral Systems,” ACM Transactions on Information smaller-and-faster-data-compression-with-zstandard/ Systems, vol. 38, no. 4, Jul. 2020. [7] J. Alakuijala, A. Farruggia, P. Ferragina, E. Kliuchnikov, R. Obryk, [27] H. Yokoo, “On the Stationary Distribution of Asymmetric Binary Z. Szabadka, and L. Vandevenne, “Brotli: A general-purpose data Systems,” in 2016 IEEE International Symposium on Information Theory compressor,” ACM Transactions on Information Systems (TOIS), vol. 37, (ISIT), July 2016, pp. 11–15. no. 1, pp. 4:1–4:30, 2018. [28] ——, “On the Stationary Distribution of Asymmetric Numeral Systems,” [8] J. Alakuijala and Z. Szabadka, “Brotli Compressed Data Format,” in 2016 International Symposium on Information Theory and Its Appli- Internet Engineering Task Force, RFC 7932, July 2016. cations (ISITA), Oct. 2016, pp. 631–635. [29] H. Yokoo and T. Shimizu, “Probability Approximation in Asymmetric [9] E. Bainville, “LZFSE,” https://github.com/lzfse/lzfse, 2017. Numeral Systems,” in 2018 International Symposium on Information [10] A. Moffat and M. Petri, “ANS-based Index Compression,” in Proceed- Theory and Its Applications (ISITA), Oct. 2018, pp. 638–642. ings of the 2017 ACM on Conference on Information and Knowledge [30] A. Weiundefinedenberger and B. Schmidt, “Massively Parallel ANS Management. ACM, Feb. 2017, pp. 677–686. Decoding on GPUs,” in Proceedings of the 48th International [11] ——, “Index Compression Using Byte-Aligned ANS Coding and Two- Conference on Parallel Processing, ser. ICPP 2019. New York, NY, Dimensional Contexts,” in Proceedings of the Eleventh ACM Interna- USA: Association for Computing Machinery, 2019. [Online]. Available: tional Conference on Web Search and Data Mining. ACM, Feb. 2018, https://doi.org/10.1145/3337821.3337888 pp. 405–413. [31] S. M. Najmabadi, Z. Wang, Y. Baroud, and S. Simon, “High Throughput [12] F. Konstantinov, G. Gryzov, and K. Bystrov, “The Use of Asymmetric Hardware Architectures for Asymmetric Numeral Systems Entropy Numeral Systems Entropy Encoding in Video Compression,” in Dis- Coding,” in 2015 9th international symposium on image and signal tributed Computer and Communication Networks, V. M. Vishnevskiy, processing and analysis (ISPA), Sept. 2015, pp. 256–259. K. E. Samouylov, and D. V. Kozyrev, Eds. Cham: Springer International [32] S. M. Najmabadi, H. S. Tungal, T.-H. Tran, and S. Simon, “Hardware- Publishing, 2019, pp. 125–139. based architecture for asymmetric numeral systems entropy decoder,” [13] J. Townsend, T. Bird, J. Kunze, and D. Barber, “HiLLoC: Lossless Image in 2017 Conference on Design and Architectures for Signal and Image Compression with Hierarchical Latent Variable Models,” in International Processing (DASIP), Sept. 2017, pp. 1–6. Conference on Learning Representations, May 2020, pp. 1–14. [33] S. M. Najmabadi, T.-H. Tran, S. Eissa, H. S. Tungal, and S. Simon, “An Architecture for Asymmetric Numeral Systems Entropy Decoder- [14] J. Townsend, T. Bird, and D. Barber, “Practical Lossless Compression A Comparison with a Canonical Huffman Decoder,” Journal of Signal with Latent Variables using Bits Back Coding,” in International Con- Processing Systems, vol. 91, no. 7, pp. 805–817, 2019. ference on Learning Representations, Apr. 2019, pp. 1–13. [34] M. J. Gormish and E. L. Schwartz, “Apparatus and method for perform- [15] D. Marwood, P. Massimino, M. Covell, and S. Baluja, “Representing ing m-ary finite state machine entropy coding,” Jun 1999, US Patent images in 200 bytes: Compression via triangulation,” in ICIP 2018. 5,912,636. Proceedings of 2018 International Conference on Image Processing, 2018, pp. 405–409. [35] D. Dubé and H. Yokoo, “Fast Construction of Almost Optimal Symbol Distributions for Asymmetric Numeral Systems,” in 2019 IEEE Interna- [16] P. A. Chou, M. Koroteev, and M. Krivokuća, “A volumetric approach tional Symposium on Information Theory (ISIT), July 2019, pp. 1682– to point cloud compression, part i: Attribute compression,” IEEE Trans. 1686. Image Process., pp. 2203–2216, 2019. [36] T. Haveliwala, S. Kamvar, D. Klein, C. Manning, and G. Golub, [17] M. Krivokuća, P. A. Chou, and M. Koroteev, “A volumetric approach to “Computing pagerank using power extrapolation,” Stanford InfoLab, point cloud compression, Part II: Geometry compression,” IEEE Trans. Technical Report 2003–45, 2003. Image Process., pp. 2217–2229, 2019. [37] A. N. Langville and C. D. Meyer, “Updating Markov chains with an [18] R. Long, M. Hernaez, I. Ochoa, and T. Weissman, “Genecomp, a new eye on Google’s PageRank,” SIAM Journal on Matrix Analysis and reference-based compressor for sam files,” in 2017 Data Compression Applications, vol. 27, no. 4, pp. 968–987, 2006. Conference (DCC), 2017, pp. 330–339. [38] S. T. Garren and R. L. Smith, “Estimating the second largest eigenvalue [19] P. Peter, “Fast inpainting-based compression: Combining shepard inter- of a Markov transition matrix,” Bernoulli, vol. 6, no. 2, pp. 215–242, polation with joint inpainting and prediction,” in ICIP 2019. Proceedings 2000. of 2019 International Conference on Image Processing, 2019, pp. 3557– 3561. [39] M. J. O’Hara, “On low-rank updates to the singular value and Tucker decompositions,” in Proceedings of the 2010 SIAM International Con- [20] Q. Chen, H. Sun, and N. Zheng, “Worst Case Driven Display Frame ference on Data Mining, 2010, pp. 713–719. Compression for Energy-Efficient Ultra-HD Display Processing,” IEEE Trans. Multimedia, vol. 20, no. 5, pp. 1113–1125, May. 2018. [21] A. Mackin, F. Zhang, and D. R. Bull, “A Study of High Frame Rate Video Formats,” IEEE Trans. Multimedia, vol. 21, no. 6, pp. 1499–1512, June 2019.