Various compression algorithms use linear transforms to represent data vectors in dierent coordinate systems, where they can be compressed better. The matrices usually have oat coecients, and the data vectors are integer, so some rounding is required. And to make the compression lossless it is required to make these transforms reversible, i.e. to be able to exactly restore the original vectors from the results of their transform. In this article we'll consider a straightforward algorithm for nding the decomposition of linear transform matrix and the approaches for estimation of the mean square approximation error and for nding the optimal decomposition, which minimizes the error.
When developing data formats for storing large volumes of multichannel data it is worth to take into account the correlation between the channels to avoid storing almost the same data again and achieve a better compression ratio. To measure the correlation and decorrelate the data we can use such methods as Karhunen-Loeve Transform (KLT) [1] to transfer the original values to a better coordinate system, where most part of coordinates would represent dierences between some original channels, which are generally small. In this article we’ll consider only integer input data. Because it is improbable, that the matrix of KLT transform would have all integer coecients, the result of the transform will be represented by real values. But it is highly ineective to immediately use ? work supported by RAS (projects: -17-117032210079-1, -19119111990037-0), RFBR (projects:18-07-00758-, 17-57-44006-ª-) and the project "Fundamentals, methods and technologies of digital monitoring and forecasting of the ecological situation in the Baikal natural territory" (2020-1902-01-071). The results were obtained using the core facilities center "Integrated Information
the oating point values to represent the real results of the transform, because the larger memory consumption of the oating point values would outweigh any advantages of the decorrelation.
So it is required to represent somehow the results of the transform by integer values to take advantage from the decorrelation. And it is not enough to just round the transform results, because it would cause information loss and may not allow to restore the original input values from the stored rounded values. Thus, we need to approximate the real-valued linear transform by some integer transform, that will be very close to the real-valued transform and will be reversible. I.e. it should be possible to restore the original values from the results of the integer approximation of the transform.
The problem of building reversible integer approximation of linear transform was thoroughly considered in the works of Hao and Shi [25]. The main idea of their approach is to represent the linear transform matrix A by a product of elementary reversible matrices (ERM), e.g. triangular matrices with diagonal elements of absolute value 1. The process of computation of the product of a vector by ERM with rounding resulting vector coordinates can be organized so, that it can be reversed. We’ll consider the idea in more details in the next section.
The article [4] is the most close to our work: here the authors prove, that matrix A with jdetAj = 1 of size n n can be decomposed into the product of permutation matrix and at most n + 1 single-row elementary reversible matrices (SERM). The algorithm of SERM decomposition by Hao and Shi is based upon high-level matrix operations. They also develop some error estimation approach, which is based upon the innity norm k k1, i.e. they estimate and try to minimize the maximum error of vector coordinate.
In this article we’ll consider a straightforward algorithm for nding the SERM decomposition of the transform matrix and the approaches for estimation of the mean square approximation error and for nding the optimal decomposition, which minimizes the error. 2
When developing some lossless data compression format, which performs a transform of integer vectors, say, the RGB color components, it is required to make the transform reversible. So, the developers of the formats use some reversible discrete approximations of the continuous transforms. 2.1
Yr = b(R + 2G + B)=4c; Ur = R
G; Vr = B
G which has the following reverse transform ("Inverse RCT"): G = Yr
b(Ur + Vr )=4c; R = Ur + G; B = Yr + G
This transform is rather rough approximation of the exact transform of the color space RGB to the color space YCrCb, which is dened by the CCIR 601 standard: After analyzing correlation between color components of particular Earth remote sensing images We have suggested the following new reversible integer transform of the color components, which allows to obtain better compression ratio for this kind of images:
Cr 1 = b(R + G + B)=3c; Cr 2 = b(R + B 2G)=2c; Cr 3 = B
R
(
C = 2Cr 2 + Cr 3 mod 2; G = Cr 1
bC=3c;
R = G + (C Cr 3)=2; B = G + (C + Cr 3)=2 (
Proof. Here we consider the operation b c as rounding down and the values like X mod D are always positive. For example 7 mod 5 = 2 + 5 = 3. So bX=Dc = (X
X mod D)=D C = 2Cr 2 + Cr 3 mod 2 = 2b(R + B = (R + B (R + B) mod 2) 2G)=2c + (B 2G + (B + R
R) mod 2 2R) mod 2 = R + B 2G bC=3c = b(R + G + B)=3c b(R + B = b(R + G + B)=3c 2G)=3c b(R + B + G)=3
Gc = G
Cr 1 = = G + (C
Cr 3)=2
G + (R + B 2G (B
R))=2 G + (C + Cr 3)=2
G + (R + B 2G + (B
R))=2 = =
R B tu
In fact, development of this kind of reversible transforms is more like a mathematical trick and it is impossible to generalize this kind of results to integer approximations of arbitrary continuous transforms. 2.3
To perform a reversible turn of a raster image by an angle a it is widely known the lifting scheme: cos sin
sin cos =
1 1 cos sin where ki 2 Z, ki 6= 0 and y0 is an auxiliary intermediate value, which will not be included into the resulting vector y. Usually the coecients ki 2 f 1; 1g, because larger absolute values will extend the range of yi, and it would make ineective the following compression steps.
Therefore we can dene reversible integer approximation of the transform
y0 = F 0(x), where the vectors x and y0 have all integer components:
y00 = k0xn + [f0(x1; :::; xn 1)]
y10 = k1x1 + [f1(y00; x2; :::; xn 1)]
yi0 = kixi + fi(y00; :::; yi0 1; xi+1; :::; xn 1)
yn0 = kny0 + fn(y10; :::; yn0 1)
(
Then we’ll be able to reverse this transform by expressing from the equations
y0 and xi and rewriting the steps in backward order:
(
When abs(ki) > 1 the integer transform becomes not perfectly reversible, because for some values of yi0 it will produce non-integer xi, but this fact makes no problem for our purposes of being able to restore the values represented by the results of the forward transform.
The error of approximation dy = y0 y results from the combination of the rounding error of fi and the errors induced by replacement of yi by yi0 in the arguments of fi. 3
The authors of [4] introduce the concepts of elementary reversible matrices (ERM), which can be computed in the reversible manner. They consider triangular elementary reversible matrices (TERM) and single-row elementary reversible matrices (ERM).
The SERM is dierent from identity matrix I in the non-diagonal cells of a single row, it can be considered as a linear shift along the corresponding to the row coordinate. They prove that the matrix A of the size n n with det A = 1 can be factorized into the product of a row permutation matrix P and n + 1 SERM matrices. The proof uses TERM factorization as an intermediate step and it is hard to use the proof as a guide for implementation of the SERM factorization algorithm. That’s why we are going to give in this section our own proof of an analogous theorem, which is straightforward and can be easily used for implementation. 3.1
Straightforward Constructive Proof of the Possibility of Matrix
Let us show constructively that for a matrix of the size n n with the absolute value of determinant 1 one can construct a factorization into n + 1 matrices B(i) of shifts along the axes of the coordinate system with unit elements on the diagonal and a single non-zero row:
M = PLAPR = B(n) ::: B(
2 1
(
Of the rows of the shift matrix B(i) only the i-th row (and the n-th for
the matrix B(0)) is dierent from the corresponding row of the identity matrix
I. So, when using this kind of decomposition each step changes a single vector
coordinate. Total the transform performs n+1 coordinate change, and it requires
n + 1 rounding operation to implement integer approximation of the transform.
The inversion of the integer approximation is performed the same way, as it was
considered before for the more general case in (
Theorem 2. To each result of permutations of rows and columns M = PLAPR
may correspond a particular factorization into n + 1 shift by coordinate matrices.
Proof. Let’s show how can we nd the single coordinate shift decomposition of
the matrix M . Consider the transform y = M x. Note, that the space of matrices
M is n2 1 - dimensional (because the matrix belongs to the subspace of the n2
- dimensional space of coecients, which is dened by the equation det M = 1
or det M = 1). The space of bij is also of the size (n + 1)(n 1) = n2 1
(where n + 1 is the number of equations (matrices) and n 1 is the number of
coecients to be found in the equations). So the dimensions of both spaces are
equal for any n. Because det B(i) = 1 for i > 0 and det B(0) = k = det M , we
have det M = det B(n) ::: B(
v = n 1 X b0j xj + kxn j=1 the n-th coordinate of R(0)x.
n 1 n 1
y1 = x1 + X b1j xj + b1nv = (1 + b1nb01) x1 + X (b1j + b1nb0j ) xj + b1nkxn (
n yi = X j=1 mij xj
xj in the expressions for y1 (
yi = i 1 X bij yj + xi + j=1 n 1 X bij xj + binv j=i+1 i 1 n = X bij X mjlxl + kxn + xi + bin @ j=0 l=1 0n 1 1
X b0j xj + kxnA
j=1
(7)
M (
i 1 i 1 yi = P P bilmlj + binb0j j=1 l=1 xj + i 1 P bilmli + binb0i + 1 xi l=1 n 1 + X i 1
X bilmlj + binb0j + bij j=i+1 l=1 ! xj + i 1 X bilmln + bin l=1 ! xn (8)
Let us suppose that as a result of analysis of the expressions for yj for j < i we will know the values of b0j for j < i, then using the rst and the last members of the expression for yi (8) we will have the system of i linear equations: 8i 1 > P bilmlj + binb0j = mij > <l=1
i 1 >> P bilmln + bin = min :l=1 for 1 j < i for i unknown variables: bil for 1 l < i and bin.
After solving the system using the second term in the expression for yi (8) we obtain: i 1 X bilmli + binb0i + 1 = mii l=1 i 1 Hence b0i = mii P bilmli 1 =bin.
l=1
Thus we have extended to the step i the inductive hypothesis that it is possible to nd b0j for j < i using the analysis of the expressions for yj for j < i. Note, that the induction base is also valid, because using the analysis of the expression for y1 we have already computed the value of b01.
i 1 X bilmlj + binb0j + bij = mij for i + 1 l=1 j < n (9) 1
After substitution to the equations (9) of the already computed values bil for l < i and bin we can nd the unknown variables bij : i 1 X bilmlj l=1 bij = mij binb0j for i + 1 j < n (10)
And after computation of b0i in (9) we can exchange i and j and also nd still unknown values: j 1 X bjlmli l=1 bji = mji bjnb0i for 1 j < i
Similarly to the already considered general case, we can group the tems with xj to obtain: As a result we have the following system of n linear equations: 8n 1 > P bnlmlj + b0j = mnj > < l=1
n 1 >> P bnlmln + k = mnn : l=1 for 1 j < n for n 1 unknown variable bnl , 1 l < n. I.e. one of the equations (14) will be redundant. The reason for the redundancy is that the condition det M = 1 = k is automatically satised due to the structure of the matrices involved in the decomposition.
Since the redundancy can consist in the fact that all the coecients in the equation with unknowns bnl, will be equal to 0 or that one of the equations will be expressed through others, we should solve the set of equations (14) by the version of Gauss algorithm for overdetermined sets of equations.
Some of the computations involved in the proof include division by values, which may be zero for some M . And solving of the systems of linear equations here may also fail. That‘s why we say " may correspond " in the theorem statement. But, once we have found a solution by the algorithm, we can be sure, that it is the only solution of the problem of nding decomposition for the matrix M . (11) (12) (13) (14)
Later on we will consider the process of nding the optimal decomposition of the matrix A by looking through all the combinations of row and column permutations PL and PR. So, when the process of nding the decomposition for a particular matrix M will fail, it will be just required to consider another version of M . 4
Let‘s recall, that M = PLAPR, i.e. it is a result of permutation of rows and columns of the matrix A. Thus, it is required to nd the best order of computation of the coecients bij . We will minimize the mean square deviation of the results of the integer approximation from the results of the original continuous transformation. For the small dimensions n we can aord to use the exhaustive search of all the (n!)2 variants of permutations. The exhaustive search will simply skip the variants, which cause division by zero during computation of the matrix M decomposition.
To select the best of the variants we should be able to estimate and compare the errors, caused by the variants of approximations. Let us consider a method for estimating the mean square dierence of the results of the integer approximation of the transform from the results of the original continuous transform. 4.1
fxg = [x] x denotes the dierence between the result of rounding and the original value.
When rounding a linear expression on integer variables with rational coefcients the size of the set of possible results of rounding depends on the least common multiple (LCM) of the denominators of the coecients. Subsequently we’ll call real rational coecients with large denominators. And the term rational will be used for the rational coecients with small denominators.
Assuming that rounding of an expression x having some real coecients is performed to the nearest integer for large ranges of x, we can consider the value fxg as a random variable, uniformly distributed over the interval [ 0:5; 0:5).
fxg2 =
Z 0:5 0:5 x2dx =
When we deal with the matrices A produced from some real-world data, say, as a result of Karhunen-Loeve Transform (KLT), the estimate (16) will always be in eect, because it is highly improbable to get rational elements here. (15) (16)
Once U and V - are independent random variables, and W = aU + bV , then
W 2 = (aU + bV )2 = a2U 2 + b2V 2 But if U = V then the rounding error will be higher:
W 2 = (aV + bV )2 = (a + b)2 V 2
Thus, to calculate the mean square rounding errors of linear expressions it is required to sum the squares of the coecients at the sub-expression errors multiplied by i2.
Let’s denote i (17) (18) (19) (20) i 1 n 1 yi = xi + X bijyj + X bijxj + binv = j=1 j=i+1
i 1 n 1 xi + X bijyj0 + X bijxj + binv0 j=1 j=i+1 i 1 X bijdyj j=1 bindv dyi = yi0 Let’s call
8i 1 n 1 9 i 1 yi = <X bijyj0 + X bijxj + binv0= + X bijdyj + bindv :j=1 j=i+1 ; j=1
8i 1 n 1 i = <X bijyj0 + X bijxj + binv0 :j=1 j=i+1 9 = ;
From (18) we can see that dyi depends on dyj for j from 0 to i 1, and also on i and 0 . We can substitute into (18) for dyj the corresponding expressions (18) to nally get: i 1 dyi = i + X Dij j j=0
dyj from (20) dyi = i + i 1 X bij j=1 j 1 j + X
l=0
X bij Dj0 + binA 0 + j=1
0i 1 0i 1 1 1
X bij Dj0 + binA 0 + bi;i 1 i 1 + j=1
i 1 i 2 X bij j + X j=1 l=1 l i 1 X bij Djl j=l+1 i 2 i 2 X bij j + X j=1 j=1 j i 1 X bilDlj l=j+1
X bij Dj0 + binA 0 + bi;i 1 i 1 + X j=1 j=1 i 2 0
1 i 1 X bilDlj A j (21) l=j+1 The formula (21) describes the process of computation of Dij using the values already computed for j from 0 to i 1.
dy1 = 1 + b1n 0
When computing the average square of the expressions (21) we will assume that i2 = 112 because we suppose here, that the coecients bij are real. 5
Let‘s consider the results of the algorithm usage for random rotation matrices. Table 1 shows some values, computed by the algorithms for a 3 3 matrix. The columns of the table correspond to the vector components. The rows of the table are:
NErr - the number of dierences of the coordinates of source vectors from the coordinates of the results of the forward transformation of the vectors by the obtained reversible integer approximation and back by its inverse detected in the loop over an area in the space of the vectors; sqD - root mean square error of the results of continuous transform from the results of its reversible integer approximation computed in the loop over the area in the space of the vectors; sqDNet - root mean square of deviation of the coordinates of the results of continuous transform from their rounded values (it allows us to estimate the lower bound for sqD); DCalc - root mean square error of the results of continuous transform from the results of its reversible integer approximation estimated by our algorithm; Top# - the permutation of columns, which gives the minimum of DCalc; Left# - the permutation of columns, which gives the minimum of DCalc;
Error = 0:650244800045056 and the following reversible integer approximation v = x2 + [0:578125087465738 x3 0:53671762651506 x1] y1 = x3 + [ 0:397744145953282 x1 0:862892315808962 v] y2 = x1 + [0:419964254965846 y1 + 0:797925922936165 v] y3 = v + [0:608264013914563 y1 0:502304008291562 y2]
Error = 0:768078871487727 and the following reversible integer approximation v = x4 + [0:2739172 x3 0:8386267 x2 0:00005832195 x1 + 0:1419508 x5] y4 = x3 + [ 0:5215370 x2 0:006188384 x1 0:04225347 x5 0:3529485 v] y1 = x2 + [0:1916207 y4
0:002794371 x1 + 0:28758131 x5 + 0:8573734 v] y5 = x1 + [0:006326542 y4 + 0:01397231 y1 0:005955142 x5 0:01004929 v] y3 = x5 + [0:05669100 y4
0:774736 y1 + 0:006511915 y5 + 0:5056881 v] y2 = v + [ 0:07052526 y4 1 2 3 4 5 6 7 NErr 0 0 0 0 0 0 0 sqD 0.3846963 0.3202919 0.3138482 0.3830634 0.394112 0.312817 0.4199523 sqDNet 0.2886750 0.2886750 0.2886689 0.2886647 0.2886753 0.2884511 0.2886757 DCalc 0.3846966 0.3202547 0.3138602 0.3830685 0.3941317 0.31305 0.5059632 Top# 7 3 5 4 2 6 1 Left# 7 6 4 3 2 1 5 M1 0.1839351 -0.3193916 0.1334853 0.113963 -0.2213811 -0.09188760 0.8808507 M2 -0.3417553 0.6397698 -0.1804105 -0.1702700 0.3857258 0.2043536 0.4709703 M3 -0.02518445 0.4264799 0.8528390 -0.08033195 -0.2374124 -0.1614647 -0.03546056 M4 -0.8741814 -0.439128 0.2010590 -0.01850904 0.04347172 -0.01722685 0.004371392 M5 0.2798122 -0.3364512 0.3273608 -0.566333 0.5675667 0.241562 0.01108168 M6 -0.02110229 -0.0127339 0.04954787 -0.03171776 -0.4384360 0.8962660 -0.02031064 M7 -0.07618058 0.02455606 -0.26880 -0.793401 -0.4791590 -0.2485626 0.02184048
and the following reversible integer approximation
The examples considered demonstrate, that the root mean square error grows steadily with the increase of the matrix size. The largest error we usually have for the vector component, which is computed last. And even for the matrix 7 7 it was slightly over 0:5, which is quite acceptable for most use cases. 6
In contrast to the previous works on the reversible integer transforms we give a straightforward proof of the possibility to construct a reversible integer approximation of a linear transform, represented by matrix with unit determinant.
We have developed a method to estimate the mean square of rounding errors for the transform matrices with real coecients and using the estimates we perform a brute-force search of the best computation order, which minimizes the error estimation.
For the matrices of rational coecients with small denominators it is still required to rectify the error estimates, but a thorough examination of the related questions would substantially increase the size of this article beyond the prescribed limits.
The algorithms suggested in the article are detailed enough to be easily implemented in software. And the tests performed show, that they allow us to compute high quality reversible approximations of linear transforms.
We have a successful experience of using the reversible integer transforms in conjunction with a specialized compression algorithm [6] for compression of multichannel images.