This is the post-print version of the following article: Kenji Nakagawa, ``On the Optimal Markov Chain of IS Simulation,'' in IEEE Transactions on Information Theory, 2001. The original publication is available athttps://doi.org/10.1109/18.904558 (DOI: 10.1109/18.904558). (c) 2001 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other users, including reprinting/republishing this material for advertising or
promotional purposes, creating new collective works for resale or redistribution to servers or lists, or reuse of any copyrighted components of this work in other works.

442

IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 1, JANUARY 2001[21] L. Gurvits, “Stability of discrete linear inclusion,” Linear Algebra Applications, vol. 231, pp. 47–85, 1995.
[22] N. E. Barabanov, “Lyapunov indicators of discrete inclusions. I.,” Automat. Remote Contr., vol. 49, pp. 152–157, Feb. 1988.
[23] J. N. Tsitsiklis and V. D. Blondel, “The Lyapunov exponent and joint
spectral radius of pairs of matrices are hard—when not impossible—tocompute and to approximate,” Math. Control, Signals, Syst., vol. 10, pp.
31–40, 1997.
[24] V. D. Blondel and J. N. Tsitsiklis, “Boundedness of finitely generated
matrix semigroups is undecidable,” preprint, Jan. 1999.
[25] M. Maesumi, “An efficient lower bound for the generalized spectral radius of a set of matrices,” Linear Algebra Applications, vol. 240, pp.
1–7, 1996.
, “Spectral radius of sets of matrices,” in Proc. SPIE, vol. 2303, San[26]
Diego, CA, July 1994, pp. 27–29.
[27] D. Colella and C. Heil, “The characterization of continuous, four-coefficient scaling functions and wavelets,” IEEE Trans. Inform. Theory, vol.
38, pp. 876–881, Mar. 1992.
[28] G. Gripenberg, “Computing the joint spectral radius,” Linear Algebra
Applications, vol. 234, pp. 43–60, 1996.
[29] D. Lind and B. Marcus, An Introduction to Symbolic Dynamics andCoding. New York: Cambridge Univ. Press, 1995.
[30] R. Karabed and P. H. Siegel, “Matched spectral-null codes for partialresponse channels,” IEEE Trans. Inform. Theory, vol. 37, pp. 818–855,
May 1991.

On the Optimal Markov Chain of IS Simulation
Kenji Nakagawa, Member, IEEE

Abstract—We investigate the importance sampling (IS) simulation forthe sample average of an output sequence from an irreducible Markov
chain. The optimal Markov chain used in simulation is known to be a
twisted Markov chain, however, the proofs in [2], [3] are very complicatedand do not give us a good perspective. We give a simple and natural
proof for the optimality of the simulation Markov chain in terms of the
Kullback–Leibler (KL) divergence of Markov chains. The performancedegradation of the IS simulation by using a not optimal simulation Markov
chain, i.e., the difference between the obtained variance and the minimum
variance is shown to be represented by the KL divergence. Moreover, weshow a geometric relationship between a simulation Markov chain andthe optimal one.
Index Terms—Importance sampling simulation, information geometry,
Kullback–Leibler (KL) divergence, Markov chain.

I. INTRODUCTION
The importance sampling (IS) simulation technique has been used
to obtain quickly an accurate estimate for a very small probability thatis not tractable by the ordinary Monte Carlo (MC) simulation. The IS
technique is widely used for various types of engineering problems,e.g., the estimation of a blocking probability in queuing system [2],[4], [7], an error rate in communications system [2], [6], etc. See [2]for an overview of the application of IS simulation.
If the target event (blocking in queuing systems, or error in communications systems) is a rare event with small probability of less thanabout 1006 , it is impossible to obtain an estimate by the ordinary MC

simulation because of the limit of simulation time and the precision
limit of pseudo-random numbers. To overcome these difficulties, a different probability distribution from the underlying probability distribution is used for simulation in order to generate more samples in thetarget event. Then the obtained value is modified by the likelihood ratio
to obtain an unbiased estimate. This estimate is called an IS estimate.
The probability distribution used for simulation is called a simulation
distribution. If the simulation distribution is appropriately chosen, thevariance of the IS estimate can be smaller than that of the MC estimate.
The simulation distribution that yields the IS estimate of the minimum
variance is referred to as the optimal simulation distribution. When weapply the IS technique to some simulation problem, it is critical to findthe optimal simulation distribution.
In the case of a Markov chain, the Markov chain which is used in theIS simulation is called a simulation Markov chain, and the optimal one
is called the optimal simulation Markov chain. It has been known [2],[3] that the optimal simulation Markov chain is unique and belongs to
the class of twisted Markov chains (TMC), but the proofs in [2], [3] arecomplicated and do not give a good perspective. In [2], the perturbation
technique is used to prove that the optimal simulation Markov chain is
a TMC. In [3], Jensen’s inequality is used in the proof and the line ofargument is elementary but complicated. Both of the proofs in [2] and[3] are lengthy.
In [10], we studied geometric properties of IS simulation and showedthat the Kullback–Leibler (KL) divergence of Markov chains plays animportant role in this problem. In this correspondence, we give a geometric view to this problem and provide a simple and natural proof
for the optimality of a simulation Markov chain in terms of KL divergence. The performance degradation of the IS simulation by using a not
optimal simulation Markov chain, i.e., the difference between the obtained variance and the minimum variance is shown to be represented
by the KL divergence. Moreover, we show a geometric relationship between a simulation Markov chain and the optimal one.
II. IMPORTANCE SAMPLING SIMULATION FOR MARKOV CHAINS
We investigate a simulation for the sample average of an output sequence from an irreducible finite-state Markov chain.
Let P0 denote an irreducible Markov chain on the state space f0; 1; . . . ; K g; K > 0. The0 state transition probability matrix of P0
is denoted by P0 = (P0 (x jx))x; x 2 , and the initial distribution is
given by the stationary distribution p0 = (p0 (x))x2 of P0 . The joint
probability of x; x0 2 is denoted by P0 (x; x0 )  p0 (x)P0 (x0 jx).
Consider a mapping f :2 ! Z , where Z denotes the set of
integers. Denote by EP [f ] the expectation of f with respect to P0

EP [f ] 

x; x

(1)

Assume EP [f ] = 0 without loss of generality.
Let x n = (x1 ; x2 ; . . . ; xn ) 2 n be a sample sequence generated
by the Markov chain P0 . We consider the probability of the following
set An :An = xn 2

n

1

01

n

n 0 1 i=1

f (xi ; xi+1 ) > c ;c > 0:

Writen

Manuscript received December 28, 1998.
The author is with Nagaoka University of Technology, Nagaoka, Niigata
940-2188 Japan (e-mail: nakagawa@nagaokaut.ac.jp).
Communicated by R. Cruz, Associate Editor for Communication Networks.
Publisher Item Identifier S 0018-9448(01)00580-6.

2P0 (x; x0 )f (x; x0 ): P0 (A ) =nx

2

P0 (xn )A

where

P0 (xn ) = p0 (x1 )P0 (x2 jx1 ) . . . P0 (xn jxn01 ):0018–9448/01$10.00 © 2001 IEEE

(2)IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 1, JANUARY 2001

We see from large deviations theory [2], [6], that the asymptotic behavior of the n is given as follows. Let R denote the set of real numbers. For t 2 R, denote by t the largest eigenvalue of a matrix

P0 (x jx)etf (x; x )
0

x; x

2

:

(3)

From the Perron–Frobenius theorem [2], we can see that t is a positive
eigenvalue with multiplicity 1. For y 2 R, the rate function I (y ) isdefined by

I (y) = sup(yt 0 log t ):t2R
The rate of the exponential decay of
tion.

n

!1 n

n

is represented by the rate func-n= I (c)(5)

We know the asymptotic behavior of n by this theorem, however,
it is often necessary to have an accurate value of n for some finite n.
So, we need to use a stochastic simulation, but if n is very small, it is
not tractable by the ordinary MC simulation. We apply the IS techniqueto the simulation of n .
Consider the IS simulation with a simulation Markov chain Q. Let
X n = (X1 ; . . . ; Xn ), n = 1; 2; . . ., be a random sequence genernnated by Q, and x n
be k independent realizations of
1 ; . . . ; xk 2
X n . The IS estimate n (Q) for n with a simulation Markov chain Q
is given by(Q) = k1

k

j

=1

0 2
0
PQ; s (x0 jx) = P0 (x 0jx) esf (x; x ) vQ;s (x ) ;Q(x jx)Q; s vQ; s (x)P (x )
1A (xnj ) 0 njn

Q(xj )

x; x0 2 :(9)

We can seex2

PQ; s (x0 jx) = 1; x 2 :We next define

JQ (y) = sup(ys 0 log Q; s );s2R

y 2 R:(10)

Then the following theorem holds.

where c is the value in (2).

nand vQ; s = (vQ; s (x))x2 its associated right eigenvector. A TMC
PQ; s concerned with the simulation Markov chain Q, which we call a
Q-TMC for short, is defined by(4)

Theorem 1 (Large Deviations Theorem [2]): We havelim 0 1 log

443

(6)

Theorem 2 (see [2]): For the variance V [ n (Q)] ofsimulation Markov chain Q, we haven(Q) by a1lim!1 0 n log V [ (Q)] = J (c)
n

n

Q

(11)

where c is the value in (2).

The optimal simulation Markov chain is defined to be the Markov
chain Q that maximizes the JQ (c).
B. Lemmas for TMC and Q-TMC
For preparation, we give some lemmas.
Let wt = (wt (x))x2 denote a left eigenvector of the matrix (3)associated with t , and recall that wt = (wt (x))x2 is a right eigenvector. The wt is assumed to be normalized as

w (x)w (x) = 1:
2
Similarly, let v= (v (x)) 2 denote a left eigenvector of thematrix (8) associated with which is normalized asv (x)v (x) = 1:
2
t

t

x

where 1A is the indicator function of the set An .

Q; s

Q; s

x

Q; s

III. THE OPTIMAL SIMULATION MARKOV CHAIN

Q; s

It is readily seen that the IS estimate (6) is unbiased [2]. Among theIS estimates, we would like to have the Markov chain Q that minimizes
the variance V [ n (Q)] of n (Q). A simulation Markov chain Q is
assumed to be of the same type of irreducible Markov chain as theunderlying Markov chain P0 .
A. Twisted Markov Chain and Q-Twisted Markov Chain

For the underlying Markov chain P0 , we define a TMC Pt ; t 2 R.
Recalling that t is the largest eigenvalue of the matrix (3), let us denoteby wt = (wt (x))x2 a right eigenvector associated with t . Then, a
TMC Pt is defined by0Pt (x0 jx) = P0 (x0 jx)etf (x; x ) wt (x ) ;x; x0 2 ; t 2 R:
t wt (x)From the definition of t and wt , one can readily see thatPt (x0 jx) = 1; x 2
x 2

(7)

holds.
Next, we define another type of TMC which plays an important role
when we evaluate the variance of an IS estimate. Consider the IS simulation with simulation Markov chain Q. For s 2 R, denote by Q; s
the largest eigenvalue of a matrix

P0 (x0 jx)2 esf (x; x )Q(x0 jx)
x; x 2

(8)

Q; s

x

Lemma 1: The stationary distribution pt = (pt (x))x2 of a TMC

Pt is given bypt (x) = wt (x)wt (x); x 2 :Furthermore, the stationary distribution pQ; s = (pQ; s (x))x2 of a
Q-TMC PQ; s is given bypQ; s (x) = vQ; s (x)vQ; s (x); x 2 :Proof: For Pt , see [10, Lemma 1]. Here we give a proof for PQ; s .
By the definition of vQ; s and vQ; s , we have
PQ; s (x0 jx)vQ; s (x)vQ; s (x)
(12)
x2
0 2
01 vQ; s (x0 )vQ; s (x) (13)
= PQ0 ((xx0jjxx)) esf (x; x ) Q;sx2= vQ; s (x0 )vQ; s (x0 ); x0 2 :(14)Lemma 2: For t = t3 that attains

sup (ct 0 log  )2t

we have EP[f ] = c. Furthermore, for s = s3 (Q) that attains
sup(cs 0 log  )2[f ] = c.
swe have EP

t

R

R

Q; s

444

IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 1, JANUARY 2001

Proof: For Pt , see [10, Lemma 2]. Here we give a proof for PQ; s .
First, we have

( 0 log  )

d
cs
ds

Q; s

=s (Q)s= c 0  01
= 0:dQ; sds s=s (Q)Q; sNext, we characterize the rate function I (c) by the KL divergence.
Lemma 5 (see [10, Theorem 4]): For t = t3 that attains

sup(ct 0 log  )
2

t

(15)
we have

Since vQ; s is a right eigenvector associated with Q; s , we havex2

P0 x0 x 2 sf (x; x )0vQ; s xe
Q x0 x

( j)
( j)

( )=

Q; s

()

vQ; s x ;

x

2 :
(16)

( ) = D(P kP0 ) = : inf[ ]= D(P kP0 )

I c

x; x

Substituting s

=

PQ; s

2

Q;s

ds

:

( ) into (17), we see from (15) that

Lemma 3: For a TMC Pt , we have the following:i) P ; s = t s0t ,ii) vP ; s (x) = wt (x)ws0t (x); x 2

( j)
( j)

x; x

(18)

2

and wt (x)ws0t (x) is its associated left eigenvector. iii) is easily
checked by i) and ii).

P

= (P (x0 jx))

x; x

2

= (Q(x0 jx)) 2
we define the KL divergence D(P kQ). Let p = (p(x)) 2 denote thestationary distribution of P and P (x; x0 )  p(x)P (x0 jx) denote theQ I (c):x; xjoint probability of x; x0 20) =P x; xx

2. Note that( 0 ) = p(x);

P x;x

x

2 :

x; x

2

P x; x0() log( jx) ;
( jx)

P x

Q x

P

Q

t

(23)

t

(24)

( )  2I (c):(25)

The equality holds if and only if Q = Pt .
Proof: From (9), we have

00log PP (x(x0 jxjx) ) = log PQ0((xx0 jjxx)) + sf (x; x0 ) + log v (x0 )00 log v (x) 0 log  : (26)
Q; s

Q; s

Q; s

Q; s

x; x

2

(20)

otherwise.

Lemma 4 (see [9]): For any Markov chains P and Q, we have( k )  0. The equality holds if and only if P = Q.

0(x; x0 ) log PP (x(x0 jxjx) )
Q; s

Pt

0

=
x; x

2

Pt

0
(x; x0 ) log PQ0((xx0 jjxx)) + sc 0 log 

Q; s

sc

0 log 

Q; s

:

(27)

= 2D(P kP0 ) 0 D(P kQ) 0 D(P kP )(28)
= 2I (c) 0 D(P kQ) 0 D(P kP ):
t

t

t

t

t

Q; s

Q; s

From (28) and Lemma 4, we have( )  2I (c):JQ c

The KL divergence is a kind of distance measure introduced on thespace of Markov chains which plays a fundamental role in stochastic
problems concerning Markov chains. For example, the KL divergence
determines the power exponent of the most powerful test function of
hypothesis testing for Markov chains [9], [11], or the error exponent of
the source coding of Markov chains [5], [11], etc.
The following is a fundamental property of the KL divergence.
D P Q

(22)

Hence, from (27) and Lemma 5(19)

If P (x0 jx) = 0 holds for any x; x0 2
with Q(x0 jx) = 0, we
write P  Q and say that P is dominated by Q. The KL divergence
D(P kQ) is defined by

1;

t

t

The equality holds if and only if P = Pt .

x( k )=(

P x; x

x; x

and

D P Q

0 ) log P (x0 jx) P (x0 jx)
P (x0 jx) P0 (x0 jx)
2= D(P kP ) + ct3 0 log ( k )=D P P0Multiplying both sides of (26) by Pt (x; x0 ) and taking summation
with x; x0 2 , we have, from Lemma 2 and (19)For two Markov chains

2

(21)

c

JQ c

P0 x0 x 2 sf (x; x )e
Pt x0 xxf

Theorem 3: For an arbitrary simulation Markov chain Q, we have

Proof: We can confirm by a straightforward computation thatt s0t is an eigenvalue of the matrix(

E

We can now give a new proof for the following theorem with the KL
divergence.

,iii) PP ; s = Ps0t .

P

Proof: For any P with EP [f ]= c, we see from (7) and Lemma 4

(17)

s3 Q

[f ] = c holds.

EP(x; x0 )f (x; x0 ) = dti.e., Pt is the dominating point [2] of the set fP jEP [f ] = cg withrespect to P0 .

By differentiating the both sides of (16) with respect to s, multiplying
by vQ; s (x), and taking summation with respect to x 2 , we havefrom Lemma 1
Q; s

t

R

(29)

The equality in (29) holds only if Q = Pt . If Q = Pt , from Lemma 3inf
2 D(P kP ) = inf
2 D(P kP 0 )= 0 (by putting s = 2t3 ):

s

R

t

Q; s

s

R

t

s

t

(30)
(31)

This completes the proof.
Corollary 1: The degradation of the variance by using a not necessarily optimal simulation Markov chain Q, i.e., the difference 2I (c) 0
JQ (c) is represented by the KL divergence

2I (c) 0 J (c) = D(P kQ) + inf
2 D(P kP ):
Q

t

s

R

t

Q; s

(32)

IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 1, JANUARY 2001

445

IV. GEOMETRIC RELATION OF P0 AND Q

1

Let us consider the set of irreducible Markov chains on the statespace
f ; ; ; K g. On the set , two mutually dual coordinate systems  and  are introduced [1], [8], [10], namely= 0 1 ...

1( ) = p(x)P (x jx)  P (x; x ); x; x 2 ; x 6= 0 (33)P (x jx)P (0jx )(x; x ) = log;
x; x 2 ; x 6= 0 (34)P (0jx)P (0j0)
where p = (p(x)) 2 is the stationary distribution of P and P (x; x0 )is the joint probability of x; x0 2 .
Differential geometric structures are defined on 1 and statisticalquantities are represented in terms of geometric quantities, see [1], [8]for a reference of geometry on 1.
A straight line in the  coordinate is called a geodesic, which is indeed a geodesic defined by an affine connection introduced on 1 [1].
A subset H of 1 of the form
H = fP 2 1j
P (x; x0 )a(x; x0 ) = b; a(x; x0 ) 2 R; b 2 Rg
 x; x0

0

0

0

0

0

0

0

0

0

x

x; x

2

is called a hyperplane. For a geodesic l and a hyperplane H , we saythat l intersects orthogonally to H at P 3 if P 3 2 l \ H and

(

x; x

2 ; x 6=0dP (t) x; x0
dt) ( (x; x0 ) 0  (x; x0 )) = 0 (35)
Q

Q

()

hold for any P t 2 l with the  coordinate P (t) and any Q1 ; Q2 2 Hwith  coordinates Q ; Q .
Then we have the following lemma.
Lemma 6 (Pythagoras, [1], [8]): Let l be a geodesic that intersects
orthogonally to a hyperplane H at P 3 . Then for any P 2 l and Q 2 H ,we have(

D QkP) = D(QkP 3 ) + D(P 3 kQ):

(36)

We will show some geometric properties of the IS simulation. First,
we notice that the set of P which satisfies EP f
c forms a hyperplane. Let us denote it by Hf  fP 2 jEP f
cg. Then we havethe following theorem.

1

[ ]=
[ ]=

Theorem 4: The TMC Pt is a geodesic and it intersects orthogonally to the hyperplane Hf at Pt . Furthermore, the Q-TMC PQ; s
is a geodesic and it intersects orthogonally to the hyperplane Hf at
PQ; s (Q) .
Proof: For Pt , see [10, Theorem 5]. Here we give a proof for
PQ; s . Let Q; s denote the  coordinate of PQ; s . Then by the definition(34) of  coordinate, it is easy to see that(

Q; s x; x0) =  0 (x; x0 ) + s( 1 (x; x0 ) 0  0 (x; x0 ));x; x0 2 ; x0 6= 0: (37)
Q;

Q;

Q;

From (37), we see that Q; s is a straight line in the  coordinate, hence
PQ; s is a geodesic. For any PQ; s with  coordinate Q; s and for any
Q1 ; Q2 2 Hf with  coordinates Q ; Q , we have from (33), (34),and (37) that

(

x; x

2 ; x 6=00dQ; s x; xds) ( (x; x0 ) 0  (x; x0 )) = 0: (38)Q

Q

Thus, from Lemma 2, we see that the geodesic PQ; s intersects orthogonally to the hyperplane Hf at PQ; s (Q) .
Corollary 2: The infimum in (32) is attained by s
we have

2I (c) 0 J (c) = D(P kQ) + D(P kPQtt= s3(Q). Thus,

Q; s(Q) ):(39)

Geometric relation of P ; Q; P ; P

Fig. 1.

in the  coordinate.

( )(See Lemma 2 for the definition of s3 Q .)Proof: By Lemmas 2, 6, and Theorem 4, for any s 2 R we have(

D Pt kPQ; s) = D(P kP
 D(P kP

(Q) ) + D(PQ; s (Q) kPQ; s ) (40)(41)
Q; s (Q) ):

t

Q; s

t

= ( )

The equality in (41) holds if and only if s s3 Q .
Theorem 5: For the underlying Markov chain P0 and a simulation
Markov chain Q and s 2 R, denote by 0 ; s ; Q ; Q; s the  coordinates of P0 ; Ps ; Q; PQ; s , respectively. Then we have a relation

Q; s= 0 +  0 
s

(42)

Q

which can be depicted in Fig. 1.
Proof: (42) is easily checked by the definition (34) of the coordinate.
From Fig. 1, we can see by geometric consideration that the optimal
simulation Markov chain Q is Q Pt and s
t3 .

=

=2

V. CONCLUSION
We investigated the IS simulation for the sample average of an output
sequence from an irreducible Markov chain. We provided a new proof
for the optimality of a simulation Markov chain by means of the KL
divergence of Markov chains. Furthermore, we gave a geometric viewto the relation of the underlying Markov chain and a simulation Markov
chain.
In a more complicated simulation problem, usually the optimal simulation distribution is searched only from among the twisted distributions. This is because the set of twisted distributions is one-dimensional
and the optimization is easy. The introduction of the geometric view to
the IS simulation can yield a simple proof for the proposition that the
optimal simulation distribution among a wider class is necessarily atwisted distribution.
REFERENCES[1] S. Amari, Differential-Geometrical Methods in Statistics (Lecture Notes
in Statistics. New York: Springer-Verlarg, 1985.
[2] J. A. Bucklew, Large Deviation Techniques in Decision, Simulation, andEstimation. New York: Wiley, 1990.
[3] J. A. Bucklew, P. Ney, and J. S. Sadowsky, “Monte Carlo simulation andlarge deviations theory for uniformly recurrent Markov chains,” J. Appl.
Probab., vol. 27, pp. 44–59, 1990.
[4] M. Cottrell, J.-C. Fort, and G. Malgouyres, “Large deviations and rare
events in the study of stochastic algorithms,” IEEE Trans. Automat.
Contr., vol. AC-28, no. 9, pp. 907–920, Sept. 1983.
[5] L. D. Davisson, G. Longo, and A. Sgarro, “The error exponent for thenoiseless encoding of finite ergodic Markov sources,” IEEE Trans. Inform. Theory, vol. IT-27, pp. 431–438, July 1981.

446

IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 1, JANUARY 2001

[6] A. Dembo and O. Zeitouni, Large Deviations Techniques and Application. Boston, MA: Jones and Bartlett, 1993.
[7] M. Devetsikiotis and J. K. Townsend, “Statistical optimization of
dynamic importance sampling parameters for efficient simulation ofcommunication networks,” IEEE/ACM Trans. Networking, vol. 1, pp.
293–305, June 1993.
[8] H. Itoh and S. Amari, “Geometry of information sources,” in Proc. 11th
SITA, 1988, pp. 57–60.
[9] K. Nakagawa and F. Kanaya, “On the converse theorem in statistical
hypothesis testing for Markov chains,” IEEE Trans. Inform. Theory, vol.
39, pp. 629–633, Mar. 1993.
[10] K. Nakagawa, “On the twisted Markov chain of importance sampling
simulation,” IEICE Trans. Fund., vol. E79-A, no. 9, pp. 1423–1428,
Sept. 1996.
[11] S. Natarajan, “Large deviations, hypotheses testing, and source coding
for finite Markov chains,” IEEE Trans. Inform. Theory, vol. IT-31, pp.
360–365, May 1985.

Quasi-Convexity and Optimal Binary Fusion forDistributed Detection with Identical Sensors in Generalized
Gaussian NoiseWei Shi, Member, IEEE, Thomas W. Sun, Student Member, IEEE, andRichard D. Wesel, Member, IEEE
Abstract—In this correspondence, we present a technique to find the
optimal threshold for the binary hypothesis detection problem withidentical and independent sensors. The sensors all use an identical andsingle thresholdto make local decisions, and the fusion center makes
a global decision based on thelocal binary decisions. For generalized
Gaussian noises and some non-Gaussian noise distributions, we show thatfor any admissible fusion rule, the probability of error is a quasi-convex
function of threshold . Hence, the problem decomposes into a series ofquasi-convex optimization problems that may be solved using well-known
techniques.
Assuming equal a priori probability, we give a sufficient condition of the
non-Gaussian noise distribution ( ) for the probability of error to bequasi-convex. Furthermore, this technique is extended to Bayes risk andNeyman–Pearson criteria. We also demonstrate that, in practice, it takes
fewer than twice as many binary sensors to give the performance of infinite
precision sensors in our scenario.
Index Terms—Distributed detection, fusion rule, generalized Gaussian
noise, hard decision, non-Gaussian noise.

I. INTRODUCTION
Consider distributed detection of s 2 f0m; mg, where the ith of nlocal sensors observes xi = s + zi with independent and identically
distributed (i.i.d.) noise zi . The ith sensor compares xi to a threshold to compute a binary decision ui as

ui =

0;

1;

when xi < when xi   .

Manuscript received November 3, 1999; revised July 18, 2000. The material
in this correspondence was presented in part at the IEEE International Symposium on Information Theory, Sorrento, Italy, June 2000 and at the 33rd Asilomar
Conference on Signals, Systems and Computers, Pacific Grove, CA, October
1999.
The authors are with the Electrical Engineering Department, University ofCalifornia, Los Angeles, CA 90095-1594 USA (e-mail: wshi@ee.ucla.edu;
wsun@ee.ucla.edu; wessel@ee.ucla.edu).
Communicated by P. A. Chou, Associate Editor for Source Coding.
Publisher Item Identifier S 0018-9448(01)00582-X.

Each binary decision ui is transmitted to a fusion center, which apn u to produce the final decisionplies a fusion rule F to k =i=1 i
F (k). The case with n = 2 and Gaussian noise is considered as an
example in [1]. This correspondence extends this specific case to thegeneral case with any number of sensors n for the family of generalized Gaussian noise and describes a technique for identifying the pair(; F ) that minimizes the probability of error. For equal a priori probabilities, a sufficient condition of the non-Gaussian noise distribution
g(x) for the probability of error to be quasi-convex is given and some
non-Gaussian noise distributions which satisfy this condition are listed.
This technique extends to Bayes risk and Neyman–Pearson criterion.
The identical threshold  in local sensors generally does not resultin an optimum system. However, the identical threshold assumption
reduces the complexity dramatically. For binary hypothesis detection,
Irving and Tsitsiklis [2] showed that no optimality is lost with identical
local detectors in a two-sensor system. Tsitsiklis [3] as well as Chen
and Papamarcou [4] showed that identical local detectors are asymptotically optimum when the number of sensors n tends to infinity. These
results provide some justification for restricting attention to identicalquantizers.
Even with identical local thresholds, the problem is still complicatedby the existence of multiple local minima. Furthermore, minimizing
probability of error is difficult because the Bayesian error probability
is not a smooth function, i.e., its first derivative is a discontinuous function. Hashlamoun and Varshney [1], [5] have overcome this difficulty
by using a smooth bound on the Bayesian error probability to approximately determine the optimal pair (; F ). Avi-Itzhak and Diep [6]provided an even tighter bound leading to a very good approximation.
Even so, these techniques still produce a minimization problem that
may have local minima, since neither convexity nor quasi-convexity
was shown in these references.
In this correspondence, we show that for any admissible fusion rule
Fi (i.e., any Fi that is optimal for at least one  ), Pe (; i) defined as theprobability of error for (; Fi ) is a quasi-convex function of  . Quasiconvexity [7] means that every sublevel set S = f : Pe (; i)  g
is convex. The admissible functions Fi are simply threshold tests of theforms = Fi (k) =0m;m;

if k < iif k  i.

(1)

Hence, the problem decomposes into a series of n quasi-convex optimization problems (one for each Fi ) that may be solved exactly using
a variety of techniques including, for example, the ellipsoid algorithm[8].
Section II shows the proof of the quasi-convexity of the probabilityof error versus  for every admissible fusion rule Fi with generalized Gaussian noise. In Section III, with equal a priori probability,
a sufficient condition of the non-Gaussian noise distribution g (x) forthe probability of error to be quasi-convex is given and some nonGaussian noise distributions that satisfy this condition are listed. Section IV shows some illustrative examples. Section V extends this technique to Bayes risk and Neyman-Pearson criterion. Section VI concludes this correspondence.
II. GENERALIZED GAUSSIAN NOISE
A. Admissible Fusion Rules
This section identifies the admissible fusion rules for generalized
Gaussian noise and shows that the probability of error versus  is quasiconvex for these rules.

0018–9448/01$10.00 © 2001 IEEE