LETTERS

Journal of the Physical Society of Japan 80 (2011) 093801
DOI: 10.1143/JPSJ.80.093801

Linear Complexity Lossy Compressor for Binary Redundant Memoryless Sources
Kazushi MIMURA
Graduate School of Information Sciences, Hiroshima City University, Hiroshima 731-3194, Japan(Received June 29, 2011; accepted July 7, 2011; published online August 17, 2011)

A lossy compression algorithm for binary redundant memoryless sources is presented. The proposed scheme is basedon sparse graph codes. By introducing a nonlinear function, redundant memoryless sequences can be compressed. We
propose a linear complexity compressor based on the extended belief propagation, into which an inertia term isheuristically introduced, and show that it has near-optimal performance for moderate block lengths.
KEYWORDS: nonlinear function, LDGM codes, lossy compression, replica method

Channel coding can be considered as a dual problem oflossy source coding.1,2) Recent research on error correcting
codes and lossy source coding has shown that the statisticalmechanical approach can be used to explain such problems.3)
Lossy compression for memoryless sources has been widely
investigated since Matsunaga and Yamamoto showed that it
is possible to attain the rate-distortion bound asymptoticallyusing low-density parity-check (LDPC) codes.4) The upper5)
and lower6,7) bounds on its rate-distortion performance of
low-density generator-matrix (LDGM) codes for lossy
compression with a given check degree are evaluated. Someother lossy compression schemes that have asymptotic
optimality have been proposed so far.8–19)
Eﬃcient compressors, on the other hand, are still in thestage of development. Some eﬃcient encoding algorithms,
which have near optimal performance, have been proposed,e.g., the nested binary linear codes,20) the inertia-termintroduced belief propagation,10,21) the survey-propagationbased message passing algorithm,22) the bit-ﬂipping-basedalgorithm,14) the exhaustive search of small words into what
an original message is divided,23) the linear-programmingbased algorithm,24) and the Markov–Chain Monte–Carlo(MCMC)-based algorithm.25)
For redundant memoryless sources, some low complexity
compressors, e.g., the near-linear complexity compressor
based on the exhaustive search of small words,23) thequadratic complexity compressor based on the bit-ﬂippingbased algorithm,14) and the MCMC-based compressor25)
whose complexity is independent of the sequence length,
have been proposed so far. One of other approaches to obtainnear-linear complexity compressors for redundant sources
is to apply the inertia-term-introduced belief propagation.
Hosaka and Kabashima have proposed an algorithm for
redundant sources, whose complexity is OðN 2 Þ.10) In thisstudy, we propose a linear complexity lossy compression
algorithm based on an inertia-term-introduced belief propagation by using a nonlinear function and a sparse matrix
such as low-density generator matrix (LDGM) codes forbinary redundant memoryless sources. This proposed
algorithm can be regarded as the perceptron-based one10)
whose edges are extremely deleted to have a ﬁnite
connectivity and has asymptotic optimality under someconstraints. We also show that it has near optimal
performance for moderate block lengths.
E-mail: mimura@hiroshima-cu.ac.jp

Let us ﬁrst provide the concepts of the rate-distortion
theory1) and some notations. Let x be a discrete randomvariable with an alphabet X . A source alphabet, a codeword
alphabet, and a reproduced alphabet are X , S, and X^ ,respectively. The compressor F encodes the M bit source
sequence x ¼ t ðx1 ; . . . ; xM Þ 2 X M into the N ð< MÞ bitcodeword  ¼ t ð1 ; . . . ; N Þ ¼ F ðxÞ 2 S N . The decompressor
G generates the M bit reproduced sequence x^ ¼t
ðx^1 ; . . . ; x^M Þ ¼ GðÞ 2 X^ M from the codeword . The coderate then becomes R ¼ N=Mð< 1Þ.
A distortion measure is a map d : X  X^ ! ½0; 1Þ. A
distortion between the sequences x ¼ t ðx1 ; . . . ; xM Þ 2 X M
and x^ ¼ t ðx^1 ; . . . ; x^M Þ 2 X^ M is measured P
by the averaged
single-letter distortion as dðx; x^ Þ ¼ ð1=MÞ M
k¼1 dðxk ; x^k Þ. A
rate distortion pair ðR; DÞ is said to be achievable if there
exists a sequence of rate distortion codes ðF ; GÞ withEx ½dðx; x^ Þ  D in the limit M ! 1. The rate distortionfunction RðDÞ is the inﬁmum of the rate R such that ðR; DÞis in the rate distortion region of the source for a given
distortion D.
We hereafter consider the binary alphabets X ¼ S ¼
X^ ¼ f1; 1g and a redundant binary memoryless source
whose distribution is given by ðX ¼ 1Þ ¼ 1  p; ðX ¼
1Þ ¼ p. The parameter p is a source bias. We use theHamming distortion^0; if x ¼ x,
^ ¼dðx; xÞ
ð1Þ
^1; if x 6¼ x,
as a distortion measure. The rate-distortion function of aBernoulli(p) i.i.d. source then becomes RðDÞ ¼ h2 ð pÞ h2 ðDÞ, where h2 denotes the binary entropy function which
is deﬁned by h2 ðxÞ ¼ x log2 ðxÞ  ð1  xÞ log2 ð1  xÞ.
Using a parameter w ¼ ðw1 ; w2 Þ 2 ðNnf0gÞ2 , we ﬁrst
introduce the nonlinear function g : N ! f1; 1g deﬁned as
1if w1 < jzj < w2gw ðzÞ ¼ gðw1 ;w2 Þ ðzÞ ¼;
ð2Þ
1 otherwise
where the operator t A denotes the transpose of A. For avector, this function acts componentwise. We here consider
the following decompressor:GðÞ ¼ gw ðAÞ;

ð3Þ

where A ¼ ðai Þ 2 f1; 0; 1gMN denotes a sparse regularmatrix whose row weight, i.e., the number of nonzero
elements, is C. Each nonzero element of the sparse matrix A
takes 1 with equiprobability. The function gw ðzÞ is adjusted

093801-1#2011 The Physical Society of Japan

J. Phys. Soc. Jpn. 80 (2011) 093801

L ETTERS

so that the bias of the reproduced sequence is close to that of
the original message as much as possible. More complex
functions can be considered as the function gw ðzÞ. As will bediscussed later, the complexity of our proposed algorithm isOðNÞ but it is proportional to 2C . So a small C is preferable;
therefore, we here consider a simple nonlinear function thatcan easily be adjusted.
The compressor is then deﬁned byF ðxÞ ¼ argmin dðx; GðsÞÞ:ð4Þ

s2f1;1gN

This compressor is identical to F ðxÞ ¼ argmaxs2f1;1gN
pðs; xÞ, which is the maximization of the following
distributionM1 Mdðs;xÞ
1 Ye¼eGk ðs;xk Þ ;ð5ÞZðÞZðÞ k¼1
P
where Gk ðs; xk Þ ¼ ð1=2Þð1  xk x^k Þ, x^k ¼ gw ð i2LðkÞ aki si Þand LðkÞ ¼ fijaki 6¼ 0g with the parameter  > 0. Here,ZðÞ denotes an normalization
constant of pðs; xÞ, which isPdeﬁned by ZðÞ ¼ s eMdðx;GðsÞÞ . The function Gk ðs; xk Þ
represents a distortion with respect to the kth bit.
We here consider a large row weight limit, which is a casewhere C ¼ N holds, to allows us to infer compression
performance. In this limit our scheme can be regarded as theperceptron-based code.9,10,17,18,26) When w2 > N, these are
equivalent to each other. To make the parameters w1 andw2 be of order unity, we introduce a constant into thedecompressor as GðÞ ¼ gw ðN 1=2 AÞ. The achievable distortion D is evaluated as D ¼ lim!1 @½f ðÞ=@ via thefree energy density f ðÞ ¼ ðMÞ1 EA;x ½ln ZðÞ, where E
denotes an expectation operator.
Applying the so called replica method, the free energy
density can be evaluated as f ðÞ ¼ 1 ð p lnfe þð1  e ÞKw g þ ð1  pÞ ln fe þ ð1  e Þð1  Kw Þg þR ln 2Þ R within the replica symmetric treatment, where2Kw ¼ fz2Rjgw ðzÞ¼1g ð2Þ1=2 ez =2 dz. The parameter Kw is
identical to the expectation value Ez ½gw ðzÞ ¼ 1 with arandom variable z which obeys the standard normal
distribution N ð0; 1Þ, which originates from the distributionof each element of N 1=2 A. It can be considered that thecompression performance is given using the distribution ofN 1=2 A and gw in this scheme.
The entropy density of pðs; xÞ is then obtained as sðÞ ¼ð@½f ðÞ=@  f Þ. The entropy density sðÞ must be nonnegative owing to the deﬁnition of pðs; xÞ; however it takes
negative values in the large  region. We therefore evaluate
the achievable distortion D at c which gives a zero entropy
density [sðc Þ ¼ 0], so this analysis is equivalent to theKrauth-Mézard approach which is a kind of one-step replica
symmetric breaking (RSB) treatment.27)
Using the zero-entropy-density criterion, minimizing theachievable distortion D ¼ lim!c @½f ðÞ=@ with respect
to w, one can obtain RðDÞ ¼ h2 ð pÞ  h2 ðDÞ, which is
identical to the rate-distortion function. From these twoconditions, i.e., the zero entropy density and the minimization of the achievable distortion, the following relationships
are obtained:D;ð6Þe ¼
1Dpðs; xÞ ,K. MIMURA

Kw ¼

pD:
1  2Dð7Þ

The deﬁnition of the compressor means that it has
exponential complexity. We then utilize a suboptimalalgorithm based on message passing to construct thecompressor.21) Instead of the maximization of pðs; xÞ weuse a symbol MAP encoding scheme, which is maximizationof marginal distribution,X
i ¼ argmax
pðs; xÞ:ð8Þ
si 2f1;1g

snsi 2f1;1gN1

To evaluate the marginal distribution we apply the belief
propagation. Since GðsÞ ¼ GðsÞ holds for any s, theexpectation value of si becomes zero. To avoid thisuncertainty we heuristically introduce an inertia term as aprior, which gives the following inertia-term-introduced
belief propagation:21)
YXeGk ðs;xk Þti0 k ðsi0 Þ;ð9Þ^tki ðsi Þ ¼
si0 2LðkÞni
ttþ1ik ðsi Þ ¼ ik ri ðsi Þi0 2LðkÞni

Y^tk0 i ðsi Þ:ð10ÞA pseudo marginal can be evaluated asYt
qtþ1
^tki ðsi Þ;i ðsi Þ ¼ i ri ðsi Þ

ð11Þ

k0 2MðiÞnk

k2MðiÞ

where ik and i denote normalization constants and MðiÞ ¼
fkjaki 6¼ 0 8g. Here, the function rit ðsi Þ is an introduced
prior as the inertia term and the superscript t represents aniteration step. We here consider that si is binary, so we cansafely put tik ðsi Þ ¼ ð1=2Þð1 þ mik ðtÞsi Þ, ^tki ðsi Þ ¼ ð1=2Þð1 þm^ ik ðtÞsi Þ, qti ðsi Þ ¼ ð1=2Þð1 þ mi ðtÞsi Þ. We here deﬁne a prior
1
as rit ðsi Þ ¼ esi tanh ½mi ðtÞ , where the parameter  ð0   < 1Þ
denotes the amplitude of the inertia term, which isheuristically chosen. When  ¼ 0 [rit ðsi Þ ¼ 1], the inertiaterm-introduced belief propagation recovers the conventional belief propagation. It should be noted that theperformance does not strongly depend on the detailed shapeof the function, if it is an increasing function. It has not yet
been investigated how the inertia term works in detail so far;however it is known that the inertia term chooses a single
peak in the calculation of the pseudo marginal.
We calculate the equations of the belief propagation,
which givesaki xk tanhVki ðtÞ2m^ ki ðtÞ ¼;
ð12Þ
1 þ xk tanh
Uki ðtÞ2X
mik ðt þ 1Þ ¼ tanhtanh1 m^ k0 i ðtÞ093801-2

k0 2MðiÞnk

þ tanhmi ðt þ 1Þ ¼ tanh

1mi ðtÞ ;X

þ tanhð13Þtanh1 m^ ki ðtÞk0 2MðiÞ1!

!

mi ðtÞ ;ð14Þ

#2011 The Physical Society of JapanL ETTERS

J. Phys. Soc. Jpn. 80 (2011) 093801

Uki ðtÞ ¼

X
si0 2LðkÞni

uw

X!
aki0 si0Y 1 þ mi0 k ðtÞsi0;2i0 2LðkÞni!
XX
Vki ðtÞ ¼
vw
aki0 si0

Shannon boundTime sharing bound
Proposed algorithm
Hosaka's algorithm

0.8

i0 2LðkÞni



si0 2LðkÞni

1

ð15Þ

0.6
RwhereK. MIMURA

0.4
p=0.6

i0 2LðkÞni

Y 1 þ mi0 k ðtÞsi0;
2i0 2LðkÞni

0.2

ð16Þ

p=0.80uw ðxÞ ¼ Iðw2 < x < w1 Þ þ Iðw1 < x < w2 Þ Iðx < w2 Þ  Iðw2 < xÞ

0.1

0.2

0.3

0.4

0.5

D

ð17Þ

þ Iðx ¼ w1 Þ  Iðx ¼ w2 Þ;
ð18Þ
and IðPÞ denotes an indicator function that takes 1 if theproposition P is true, and 0 otherwise. After tm iterations, theith bit of the codeword can be obtained as i ¼ sgn½mi ðtm Þusing the mean of the pseudo marginal mi ðtm Þ. To derivetheseP iterative equations, weP use the identity x^k ¼
uw ð i0 2LðkÞni aki0 si0 Þ þ aki si vw ð i0 2LðkÞni aki0 si0 Þ, which holds
for any i 2 LðkÞ.
The computational cost of the terms Uki ðtÞ and Vki ðtÞ isOð2C Þ, which depends only on the row weights C, namely, itis Oð1Þ with respect to N. The complexity of this algorithm
is therefore OðNÞ when the number of iterations tm is ﬁxed.
Utilizing eqs. (6) and (7) which are obtained in the largerow-weight-limit analysis, we can approximately set allparameters C, w1 , w2, and  of our scheme with ﬁnite row
weights except .
We ﬁrst consider a setting of the parameter . Usingeq. (6) and the rate-distortion function, we set  as ¼ c ð p; RÞ for the given the source bias p and the code
1rate R, where c ð p; RÞ ¼ lnð½h1
 1Þ. Here,2 ðh2 ð pÞ  RÞ
1h2 denotes the inverse function of the binary entropyfunction.
We next consider a setting of the parameters C, w1 , andw2 . Each element of the vector As is the summation ofC Bernoulli random variables 1  2Bernoullið0:5Þ, wheres 2 f1; 1gN denotes a candidate of a codeword. This is asimilar situation to the row weight limit. To keep the row
weight ﬁnite, we restrict the row weight as C  Cmax . Usingeq. (7) and the rate-distortion function, we set ðC; w1 ; w2 Þ
^ 0 0 w02 Þ 00asðC; w1 ; w2 Þ ¼ argminðC0 ;w
 
P1 ;w2 Þ2DðCmax Þ jKðC ; w1 ;C^Kð p;RÞj, where KðC;w1;w2 Þ ¼ n2f0;...;Cg:gw ðC2nÞ¼1 2 Cn ,
1Kð p; RÞ ¼½ p  h12 ðh2 ð pÞ  RÞ=½1  2h2 ðh2 ð pÞ  RÞ, andDðCmax Þ ¼ fðC; w1 ; w2 Þj2  C  Cmax ; 0 < w1  C  1; w1 <w2  C þ 1g for the given p and R. Note that the parameterthat gives second smallest value might provide better
performance.
Lastly,  is determined by trial and error. In this study, wechoose  only within f0:2; 0:3; 0:4; 0:5g.
The empirical compression performance is shown inFig. 1. In this ﬁgure, the distortion averaged over 10 runs isplotted as a function of the code rate R for the source bias
p 2 f0:6; 0:7; 0:8; 0:9g. The length of codewords is fexed at
N ¼ 420, and the length of original sequence is adjusted.
We here choose Cmax ¼ 8 and p > 0:5. In this ﬁgure, the1

Shannon boundTime sharing bound
Proposed algorithm
Hosaka's algorithm

0.80.6R Iðw1 < x < w1 Þ;vw ðxÞ ¼ Iðx ¼ w2 Þ  Iðx ¼ w1 Þ00.4
p=0.7

p=0.9

0.2

00

0.1

0.2

0.3

0.4

0.5

DFig. 1. Empirical compression performance against the code rate R for
typical source bias p. The proposed algorithm (squares) and Hosaka’salgorithm (circles) are shown. The length of the orginal sequences isN ¼ 420, and all the measurements are averaged over 10 runs. Theparameter  is chosen within f0:2; 0:3; 0:4; 0:5g. The row weight C is chosen
within C  8 (Cmax ¼ 8). Top: p 2 f0:6; 0:8g. Bottom: p 2 f0:7; 0:9g.

time sharing bound is also shown. The time sharing boundis given by RðDÞ ¼ ð1  D=pÞh2 ð pÞ, which denotes thecompression performance achieved by the time sharing
scheme of a lossless coding [Rð0Þ ¼ h2 ð pÞ] and a trivial
encoding that always outputs an all-one vector for any input[Rð pÞ ¼ 0]. It can be conﬁrmed that the proposed linear
complexity compressor [with OðNÞ complexity] has slightly
better performance than Hosaka’s algorithm [with OðN 2 Þcomplexity].
We observed a not-so-good performance for the small-p
region that the source bias is less than about 0.2. In this
^ 0 ; w01 ; w02 Þ  Kð p; RÞj is not much smaller
region minjKðCthan that of the large-p region. When we compress theoriginal sequence x of which bias is p < 0:5, we can ﬁrst ﬂip
it as x and then compress. The information for determining
whether the sequence ﬂips requires one bit. To reduce^ 0 ; w01 ; w02 Þ  Kð p; RÞj, it might be helpful to introminjKðC
duce more complex nonlinear functions.
In this study, we have proposed a scheme using anonlinear function and a sparse matrix, and as well as alinear complexity message passing compressor based on theinertia-term-introduced belief propagation. The proposed
method can treat redundant memoryless sources and has
near-optimal compression performance for moderate block093801-3#2011 The Physical Society of Japan

J. Phys. Soc. Jpn. 80 (2011) 093801

L ETTERS

lengths. The adjustment of the column weight distribution of
the sparse matrix might enable us to improve the compression performance. The analysis of this scheme with ﬁnite
row weights is one of our future studies.
Acknowledgments The author would like to thank Tadaaki Hosaka,
Tatsuto Murayama, and Jun Muramatsu for their valuable comments and Ankit
Gupta, Sergio Verdú, and Tsachy Weissman for providing him valuable
preprints. This work was partially supported by a Grant-in-Aid for ScientiﬁcResearch (C) No. 22500136 from the Ministry of Education, Culture, Sports,
Science and Technology of Japan (MEXT).

1) T. M. Cover and J. A. Thomas: Elements of Information Theory(Wiley, New York, 2006) 2nd ed.
2) I. Csiszár and J. Körner: Information Theory: Coding Theorems for
Discrete Memoryless Systems (Academic Press, New York, 1981).
3) H. Nishimori: Statistical Physics of Spin Glasses and Information
Processing—An Introduction (Oxford University Press, Oxford, U.K.,
2001).
4) Y. Matsunaga and H. Yamamoto: IEEE Trans. Inf. Theory 49 (2003)
2225.
5) A. G. Dimakis, M. J. Wainwright, and K. Ramchandran: Proc.
Information Theory Workshop, 2005, p. 650.
6) E. Martinian and M. J. Wainwright: Proc. Workshop Information
Theory and Its Applications, 2006.
7) M. J. Wainwright: IEEE Signal Process. Mag. 24 [5] (2007) 47.
8) T. Murayama and M. Okada: J. Phys. A 36 (2003) 11123.
9) T. Hosaka, Y. Kabashima, and H. Nishimori: Phys. Rev. E 66 (2002)

K. MIMURA066126.
10) T. Hosaka and Y. Kabashima: Physica A 365 (2006) 113.
11) S. Ciliberti and M. Mézard: J. Stat. Mech. 2005 (2005) P10003.
12) S. Ciliberti, M. Mézard, and R. Zecchina: Complex Syst. Methods 3(2006) 58.
13) K. Mimura and M. Okada: Phys. Rev. E 74 (2006) 026108.
14) A. Gupta and S. Verdú: IEEE Trans. Inf. Theory 55 (2009) 1961.
15) S. Miyake and J. Muramatsu: IEICE Trans. Fundam. Electron.
Commun. Comput. Sci. E91-A (2008) 1488.
16) J. Muramatsu and S. Miyake: Proc. IEEE Int. Symp. Information
Theory (ISIT2008), 2008, p. 424.
17) F. Cousseau, K. Mimura, T. Omori, and M. Okada: Phys. Rev. E 78
(2008) 021124.
18) F. Cousseau, K. Mimura, and M. Okada: Proc. IEEE Int. Symp.
Information Theory (ISIT2008), 2008, p. 509.
19) K. Mimura: J. Phys. A 42 (2009) 135002.
20) T. Wadayama: Proc. 3rd Int. Symp. Turbo Codes and Related Topics,
2003, p. 231.
21) T. Murayama: Phys. Rev. E 69 (2004) 035105(R).
22) M. J. Wainwright and E. Maneva: Proc. IEEE Int. Symp. Information
Theory (ISIT2005), 2005, p. 1493.
23) A. Gupta, S. Verdú, and T. Weissman: Proc. IEEE Int. Symp.
Information Theory (ISIT2008), 2008, p. 847.
24) S. Miyake, J. Honda, and H. Yamamoto: Proc. IEEE Int. Symp.
Information Theory and Its Applications (ISITA2008), 2008, p. 589.
25) S. Jalali and T. Weissman: Proc. IEEE Int. Symp. Information Theory(ISIT2008), 2008, p. 852.
26) K. Mimura, F. Cousseau, and M. Okada: J. Phys. Soc. Jpn. 80 (2011)034802.
27) W. Krauth and M. Mézard: J. Phys. (Paris) 50 (1989) 3057.

093801-4#2011 The Physical Society of Japan