JMLR: Workshop and Conference Proceedings 4: 5-20

New challenges for feature selection

Approximating Mutual Informationby Maximum Likelihood Density Ratio Estimation
Taiji Suzuki

s-taiji@stat.t.u-tokyo.ac.jp
Department of Mathematical Informatics, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan

Masashi Sugiyama

sugi@cs.titech.ac.jp

Department of Computer Science, Tokyo Institute of Technology,
2-12-1 O-okayama, Meguro-ku, Tokyo 152-8552, JapanJun Sesesesejun@is.ocha.ac.jp

Department of Information Science, Ochanomizu University,
2-1-1 Ohtsuka, Bunkyo-ku, Tokyo 112-8610, Japan

Takafumi Kanamori

kanamori@is.nagoya-u.ac.jp

Department of Computer Science and Mathematical Informatics,
Nagoya University, Furocho, Chikusaku, Nagoya 464-8603, Japan

Editor: Saeys et al.

Abstract
Mutual information is useful in various data processing tasks such as feature selection or independent component analysis. In this paper, we propose a new method of approximating
mutual information based on maximum likelihood estimation of a density ratio function.
Our method, called Maximum Likelihood Mutual Information (MLMI), has several attractive properties, e.g., density estimation is not involved, it is a single-shot procedure,
the global optimal solution can be efficiently computed, and cross-validation is available formodel selection. Numerical experiments show that MLMI compares favorably with existing
methods.

1. Introduction
Detection of dependencies between random variables is highly useful in various machinelearning problems such as feature selection (Guyon and Elisseeff, 2003; Torkkola, 2003) andindependent component analysis (Comon, 1994). Although classical correlation analysis
would be still useful in these problems, it cannot be used for discovering non-linear dependencies with no correlation. On the other hand, mutual information (MI), which plays
an important role in information theory (Cover and Thomas, 1991), allows us to identify
general nonlinear dependencies. MI is defined by
I(X, Y ) :=ZZ

pxy (x, y) log



pxy (x, y)px (x)py (y)



dxdy,

(1)

and it vanishes if and only if x and y are independent. For this reason, estimating MI fromsamples has gathered a lot of attention for many years.
c 2008 Suzuki, Sugiyama, Sese and Kanamori.

Suzuki, Sugiyama, Sese and Kanamori

A naive approach to estimating MI is to use a kernel density estimator (KDE) (Silverman, 1986; Fraser and Swinney, 1986), i.e., the densities pxy (x, y), px (x), and py (y)are separately estimated from samples and the estimated densities are used for computing MI. The bandwidth of the kernel functions could be optimized based on likelihoodcross-validation (Härdle et al., 2004), so there remains no open tuning parameter in thisapproach. However, density estimation is known to be a hard problem and division by estimated densities is involved when approximating MI, which tend to magnify the estimation
error. Therefore, the KDE-based method may not be reliable in practice.
Alternative methods involve estimation of the entropies using k-nearest neighbor (KNN)
samples (Kraskov et al., 2004) or using the Edgeworth (EDGE) expansion (Hulle, 2005).
The KNN-based approach was shown to perform better than KDE (Khan et al., 2007),
given that the number k is chosen appropriately—a small (large) k yields an estimator withsmall (large) bias and large (small) variance. However, appropriately determining the valueof k so that the bias-variance trade-off is optimally controlled is not straightforward in thecontext of MI estimation. The EDGE method works well when the target density is close
to the normal distribution; otherwise it is biased and therefore not reliable.
In this paper, we propose a new MI estimator that can overcome the limitations of theexisting approaches. Our method, which we call Maximum Likelihood Mutual Information(MLMI), does not involve density estimation and directly models the density ratio:
w(x, y) :=pxy (x, y).
px (x)py (y)(2)

Thus it is a single-shot procedure without division by estimated quantities and therefore theestimation error is not further expanded. The density ratio is estimated by the maximum
likelihood method and it is cast as a convex optimization problem. Therefore, the uniqueglobal optimal solution can be obtained efficiently. Furthermore, cross-validation (CV) is
available for model selection, so the values of tuning parameters such as the kernel width canbe adaptively determined in an objective manner. Our method does not assume normalityof the target distribution and therefore is flexible. Numerical experiments show that MLMI
compares favorably with existing methods.
A sibling of MLMI, called LSMI, is presented in Suzuki et al. (2008). In that paper,
we used a least-squares method for density ratio based MI estimation and emphasized its
practical usefulness in variable selection. On the other hand, the current paper employsa maximum likelihood method and focuses more on mathematical aspects of density ratio
based MI estimation.

2. Approximating MI by Maximum Likelihood Density Ratio Estimation
In this section, we formulate the MI approximation problem as a density ratio estimation
problem and propose a new MI estimation method.
2.1 Formulation
Let DX (⊂ Rdx ) and DY (⊂ Rdy ) be data domains and suppose we are given n independentand identically distributed (i.i.d.) paired samples{(xi , yi ) | xi ∈ DX , yi ∈ DY }ni=16Approximating Mutual Information by Maximum Likelihood Density Ratio Estimation

drawn from a joint distribution with density pxy (x, y). Let us denote the marginal densities
of xi and yi by px (x) and py (y), respectively. The goal is to estimate MI defined by Eq.(1).
Our approach here is to estimate the density ratio w(x, y) defined by Eq.(2); then MI
can be approximated using a density ratio estimator w(x,
b
y) by
1
b
I(X,
Y ) :=nnX
i=1

log w(xb i , yi ).

We model the density ratio function w(x, y) by the following linear model:w(x,by) := α⊤ ϕ(x, y),(3)where α := (α1 , α2 , . . . , αb )⊤ are parameters to be learned from samples, ⊤ denotes thetranspose of a matrix or a vector, andϕ(x, y) := (ϕ1 (x, y), ϕ2 (x, y), . . . , ϕb (x, y))⊤are basis functions such thatϕ(x, y) ≥ 0b for all (x, y) ∈ DX × DY .
0b denotes the b-dimensional vector with all zeros. Note that ϕ(x, y) could be dependent
on the samples {xi , yi }ni=1 , i.e., kernel models are also allowed. We explain how the basis
functions ϕ(x, y) are designed in Section 2.4.
2.2 Maximum Likelihood Estimation of Density Ratio Function
Using an estimated density ratio w(x,by), we may estimate the joint density pxy (x, y) by
pbxy (x, y) := w(x,by)px (x)py (y).

Based on this, we determine the parameter α in the model w(x,by) so that the following
log-likelihood L is maximized:
L(α) :=nX
i=1

log pbxy (xi , yi ) =

n
X
i=1

n
n

 X
X
log py (yi ).
log px (xi ) +log α⊤ ϕ(xi , yi ) +
i=1

i=1

The second and the third terms are constants and therefore can be safely ignored. This is
our objective function to be maximized with respect to the parameters α, which is concave.
Note that this corresponds to an empirical approximation of the Kullback-Leibler divergencefrom pxy (x, y) to pbxy (x, y) up to some irrelevant constant.
w(x,by) is an estimator of the density ratio w(x, y) which is non-negative by definition.
Therefore, it is natural to impose w(x,by) ≥ 0 for all (x, y) ∈ DX × DY , which can beachieved by restricting α ≥ 0b . In addition to non-negativity, w(x,by) should be properly
normalized since pbxy (x, y) is a probability density function:
Z
Z
X
1
1 = pbxy (x, y)dxdy = α⊤ ϕ(x, y)px (x)py (y)dxdy ≈α⊤ ϕ(xi , yj ), (4)n(n − 1)1≤i6=j≤n7

Suzuki, Sugiyama, Sese and Kanamori

where we used the U -statistic (Serfling, 1980, p.171) for obtaining the empirical estimator.
Now our optimization criterion is summarized as follows.
#
" n


X
⊤
maximizelog α ϕ(xi , yi )α∈Rbsubject toi=11n(n − 1)X

1≤i6=j≤nα⊤ ϕ(xi , yj ) = 1 and α ≥ 0b .

(5)

We call the above method Maximum Likelihood Mutual Information (MLMI).
MLMI can be characterized by Legendre-Fenchel duality of the convex function ‘− log’(Rockafellar, 1970; Boyd and Vandenberghe, 2004). For f (u) = − log(u), MI is expressedasZpx (x)py (y)I(X, Y ) = pxy (x, y)fdxdy.
pxy (x, y)
Let f ∗ be the Legendre-Fenchel dual function of f , which is defined and given byf ∗ (v) := sup {uv − f (u)} = −1 − log(−v)(for v < 0).

u∈RThen Legendre-Fenchel duality implies that I(X, Y ) is obtained by solving the following
concave maximization problem (Nguyen et al., 2008):Z
px (x)py (y)∗− f (−w(x,by)) dxdy
I(X, Y ) = suppxy (x, y) −w(x,by)
pxy (x, y)
w≥0
b

Z 

= sup
− w(x,by)px (x)py (y) + log w(x,by)pxy (x, y) dxdy + 1 ,w≥0b

where the supremum is taken over all non-negative measurable functions. If the linear model
assumption (3) and the normalization constraint (4) are imposed and the expectation isapproximated by the sample average, the above formulation is reduced to MLMI.
2.3 Convergence Bound

Here, we show a non-parametric convergence rate of the solution of the optimization problem(5). The set of basis functions is denoted by
F := {ϕθ | θ ∈ Θ},
where Θ is some parameter or index set. The set of basis functions used for estimation atn samples is characterized by a subset of the parameter set Θn ⊆ Θ and denoted by
Fn := {ϕθ | θ ∈ Θn } ⊂ F,
which can behave stochastically. The set of finite linear combinations of F with positivecoefficients and its bounded subset are denoted by()
X
αl ϕθl αl ≥ 0, ϕθl ∈ F , G M := {g ∈ G | kgk∞ ≤ M } ,
G :=
l

8

Approximating Mutual Information by Maximum Likelihood Density Ratio Estimationand their subsets used for estimation at n samples are denoted by()
X
αl ϕθl αl ≥ 0, ϕθl ∈ Fn ⊂ G, GnM := {g ∈ Gn | kgk∞ ≤ M } ⊂ G M .
Gn :=l

Let Gbn be the feasible set of MLMI:(
Gbn := g ∈ Gn

1n(n − 1)X)

g(xi , yj ) = 1 .

1≤i6=j≤nThen, the solution gbn of (generalized) MLMI is given as follows:" n#
1X
gbn := arg max
log g(xi , yi ) .
n
g∈Gbn
i=1

For simplicity, we assume that the optimal solution can be uniquely determined. In order toderive the convergence rates of MLMI, we make the following assumptions. Let g0 denote
g0 (x, y) :=pxy (x, y).
px (x)py (y)We define the (generalized) Hellinger distance with respect to px py as′hQ (g, g ) :=Zp√
( g − g ′ )2 px (x)py (y)dxdy

1/2

,

where g and g ′ are non-negative measurable functions (not necessarily probability densities).
Assumption 1
1. On the support of pxy , there exists a constant η1 < ∞ such thatg0 ≤ η1 .
2. All basis functions are non-negative, and there exist constants ǫ0 , ξ0 > 0 such thatZ
ϕ(x, y)px (x)py (y)dxdy ≥ ǫ0 , kϕk∞ ≤ ξ0 , (∀ϕ ∈ F).
3. There exist constants 0 < γ < 2 and K such thatlog N[] (ǫ, G M , hQ ) ≤ K

√

M
ǫ

!γ

,

(6)

where N[] (ǫ, F, d) is the ǫ-bracketing number of F with norm d (van de Geer, 2000)1 .
Then we obtain the following theorem and corollary (see Appendix for a sketch of proof).
1. As shown in Section 2.4, we will use Gaussian kernel models in practice. The bracketing number for
Gaussian mixture models is extensively discussed, e.g., in Ghosal and van der Vaart (2001).

9

Suzuki, Sugiyama, Sese and Kanamori

Theorem 1 Let
gn∗ := max
g∈Gbn

Z

log(g(x, y))p(x, y)dxdy.

In addition to Assumption 1, if there exist constants e
c0 , c0 such that gn∗ satisfies


g0
P e
c0 ≤ ∗ ≤ c0 → 1,
gn

then

hQ (b
gn , g0 ) = Op (n

1
− 2+γ

(7)

+ hQ (gn∗ , g0 )),where Op denotes the asymptotic order in probability.
Corollary 2 If there exists N such that g0 ∈ Gn (∀n ≥ N ), thenhQ (bgn , g0 ) = Op (n

1
− 2+γ

).

In Sugiyama et al. (2008), a similar convergence result to the above theorem has been
provided, where g0 was assumed to be bounded from both below and above. On the other
hand, the current proof only requires an upper bound on g0 (see Assumption 1.1). Thus,
the above theorem is more general than the result in Sugiyama et al. (2008) in terms of g0 .
More technically, we use the bracketing number with respect to the Hellinger distancefor describing the complexity of a function class, while Sugiyama et al. (2008) used thebracketing number with respect to the ℓ2 distance or the covering number. The main mathematical device used in Sugiyama et al. (2008) was to bound hQ (bgn , gn∗ ) by the difference∗
gn + gn∗ )). On the other hand,
between the empirical mean and the expectation of log(2gn /(b∗∗the current paper employs 2b
gn /(bgn + gn ) instead of log(2gn /(bgn + gn∗ )), which enables us
to replace the lower bound of g0 with the bounds of the ratio g0 /gn∗ . Then we utilize theconvexity of Gbn and the bracketing number condition with respect to the Hellinger distance(see Section 7 of van de Geer (2000) for more details).
2.4 Model Selection by Likelihood Cross Validation
The performance of MLMI depends on the choice of the basis functions ϕ(x, y). Here weshow that model selection can be carried out based on a variant of CV.
First, the samples {zi | zi = (xi , yi )}ni=1 are divided into K disjoint subsets {Zk }K
k=1 of(approximately) the same size. Then a density ratio estimator w
bZk (x, y) is obtained using{Zj }j6=k (i.e., without Zk ) and the log-likelihood for the hold-out samples Zk is computedasX
1(K-CV)LZk=
log w
bZk (x′ , y ′ ),|Zk | ′ ′(x ,y )∈Zk

where |Zk | denotes the number of sample pairs in the set Zk . This procedure is repeatedfor k = 1, 2, . . . , K and its average L(K-CV) is outputted:
K

L(K-CV) =1 X (K-CV)
LZk
.
K
k=1

10

Approximating Mutual Information by Maximum Likelihood Density Ratio Estimation

For model selection, we compute L(K-CV) for all model candidates (the basis functions
ϕ(x, y) in the current setting) and choose the one that maximizes the hold-out loglikelihood. Note that L(K-CV) is an almost unbiased estimate of the Kullback-Leibler divergence from pxy (x, y) to pbxy (x, y) up to some irrelevant constant, where the ‘almost’-nesscomes from the fact that the number of samples is reduced in the CV procedure due to data
splitting (Schölkopf and Smola, 2002).
A good model may be chosen by the above CV procedure, given that a set of promisingmodel candidates is prepared. As model candidates, we propose using a Gaussian kernel
model: for z = (x⊤ , y ⊤ )⊤ ,kz − cℓ k2
ϕℓ (x, y) = exp −
2σ 2







kx − uℓ k2
ky − vℓ k2= exp −
exp −
,
2σ 2
2σ 2

⊤ ⊤ b
where {cℓ | cℓ = (u⊤ℓ , vℓ ) }ℓ=1 are Gaussian centers; we choose the centers randomly from⊤ ⊤ n{zi | zi = (x⊤i , yi ) }i=1 . Below, we fix the number of basis functions at b = min(200, n)
and choose the Gaussian width σ by CV.

3. Relation to Existing Methods
In this section, we discuss the characteristics of existing and proposed approaches.
3.1 Kernel Density Estimator (KDE)
KDE is a non-parametric technique to estimate a probability density function p(x) defined
on Rd from its i.i.d. samples {xi }ni=1 . For the Gaussian kernel, KDE is expressed as
pb(x) =

nXkx − xi k2
1
exp
−
.
2σ 2
n(2πσ 2 )d/2 i=1

The performance of KDE depends on the choice of the kernel width σ and it can be
optimized by likelihood CV as follows (Härdle et al., 2004): First, divide the samples {xi }ni=1into K disjoint subsets {Xk }K
k=1 of (approximately) the same size. Then obtain a density
estimate pbXk (x) from {Xj }j6=k (i.e., without Xk ) and compute its hold-out log-likelihood forXk :1 X
log pbXk (x).
|Xk |
x∈Xk

This procedure is repeated for k = 1, 2, . . . , K and choose the value of σ such that the
average of the hold-out log-likelihood over all k is maximized. Note that the average holdout log-likelihood is an almost unbiased estimate of the Kullback-Leibler divergence fromp(x) to pb(x) up to some irrelevant constant.
Based on KDE, MI can be approximated using density estimates pbxy (x, y), pbx (x), andpby (y) (obtained from {(xi , yi )}ni=1 , {xi }ni=1 , and {yi }ni=1 , respectively) as
1b
I(X,Y ) :=
n

n
X

log

i=1

11

pbxy (xi , yi ).
pbx (xi )bpy (yi )Suzuki, Sugiyama, Sese and Kanamori

However, density estimation is known to be a hard problem and division by estimateddensities may expand the estimation error. For this reason, the KDE-based approach may
not be reliable in practice.
3.2 K-nearest Neighbor Method (KNN)
MI can be expressed in terms of the entropies as
I(X, Y ) = H(X) + H(Y ) − H(X, Y ),where H(X) denotes the entropy of X:
H(X) := −Z

px (x) log px (x)dx.

Thus MI can be approximated if the entropies H(X), H(Y ), and H(X, Y ) are estimated.
Kraskov et al. (2004) developed an entropy estimator that utilizes the k-nearest neighbordistance (KNN). Let us define the norm of z = (x, y) by kzkz := max{kxk, kyk}, wherek · k denotes the Euclidean norm. Let Nk (i) be the set of k-nearest neighbor samples of(xi , yi ) with respect to the norm k · kz , and let
ǫx (i) := max{kxi − xi′ k | (xi′ , yi′ ) ∈ Nk (i)},nx (i) := |{zi′ | kxi − xi′ k ≤ ǫx (i)}|,ǫy (i) := max{kyi − yi′ k | (xi′ , yi′ ) ∈ Nk (i)},ny (i) := |{zi′ | kyi − yi′ k ≤ ǫy (i)}|.

Then the KNN-based MI estimator is given by
1
1
b
I(X,
Y ) :=ψ(k) + ψ(n) − −
k n

n h
X
i=1

i
ψ(nx (i)) + ψ(ny (i)) ,where ψ is the digamma function.
An advantage of the above KNN-based method is that it does not simply replace entropies with their estimates, but it is designed to cancel the error of individual entropy
estimation. A practical drawback of the KNN-based approach is that the estimation accuracy depends on the value of k and there seems no systematic strategy to choose the valueof k appropriately.
3.3 Edgeworth Expansion Method (EDGE)
Hulle (2005) proposed an entropy approximation method based on the Edgeworth expansion,
where the entropy of a distribution is approximated by that of the normal distribution
and some additional higher-order correction terms. More specifically, for a d-dimensionalb of the entropy H is given by
distribution, an estimator H
b =Hnormal − 1
H
12

d
X
i=1

κ2i,i,i −

1
4

d
X

i,j=1,i6=jκ2i,i,j −

1
72

d
X

κ2i,j,k ,i,j,k=1,i<j<kwhere Hnormal is the entropy of the normal distribution with covariance matrix equal to
the target distribution and κi,j,k (1 ≤ i, j, k ≤ d) is the standardized third cumulant of thetarget distribution. In practice, all the cumulants are estimated from samples.
12

Approximating Mutual Information by Maximum Likelihood Density Ratio Estimation

Table 1: Relation among existing and proposed MI estimators. If the order of the Edgeworth expansion is regarded a tuning parameter, model selection of EDGE shouldbe ‘Not available’.

KDE
KNN
EDGE
MLMI

Division by estimated quantities
Involved
Not involved
Not involved
Not involved

Model selection
Available
Not available
Not necessary
Available

Distribution
Free
Free
Nearly normal
Freebb ), andBased on EDGE, MI can be approximated using entropy estimates H(X),
H(Yb
H(X, Y ) asbbb ) − H(X,b
I(X,Y ) := H(X)
+ H(Y
Y ).

If the underlying distribution is close to the normal distribution, the above approximationis quite accurate and the EDGE method works very well. However, if the distributions
are far from the normal distribution, the approximation error gets large and therefore the
EDGE method is unreliable.
In principle, it is possible to include the fourth and even higher cumulants for further
reducing the estimation bias. However, this in turn increases the estimation variance; theexpansion up to the third cumulants seems to be reasonable.
3.4 Discussions
The characteristics of the proposed and existing MI estimators are summarized in Table 1.
KDE is distribution-free and model selection is possible by CV. However, division by estimated densities is involved, which can result in magnifying the estimation error. KNN isdistribution-free and does not involve division by estimated quantities. However, there is
no model selection method for determining the number of nearest neighbors and therefore
its practical performance is not reliable. EDGE does not involve division by estimatedquantities and any tuning parameters. However, it is based on the assumption that thetarget distributions are close to the normal distribution and the result is not reliable if thisassumption is violated. MLMI is distribution-free, it does not involve division by estimatedquantities, and model selection is possible by CV. Thus MLMI overcomes the limitationsof the existing approaches.

4. Numerical Experiments
In this section, we experimentally investigate the performance of the proposed and existingMI estimators using artificial datasets. We use the following four datasets (see Figure 1):(a) Linear dependence: y has a linear dependence on x as
x ∼ N (x; 0, 0.5) and y|x ∼ N (y; 3x, 1),
where N (x; µ, σ 2 ) denotes the normal density with mean µ and variance σ 2 .
13

Suzuki, Sugiyama, Sese and Kanamori(b) Non-linear dependence with correlation: y has a quadratic dependence on x as
x ∼ N (x; 0, 1) and y|x ∼ N (y; x2 , 1).
(c) Non-linear dependence without correlation: y has a lattice-structured dependence on x as(
N (x; 0, 1/3)
if x ≤ |1/6|,
x ∼ U(x; −0.5, 0.5) and y|x ∼0.5N (x; 1, 1/3) + 0.5N (x; −1, 1/3) otherwise,
where U(x; a, b) denotes the uniform density on (a, b).
(d) Independence: x and y are independent to each other asx ∼ U(x; 0, 0.5) and y|x ∼ N (y; 0, 1).
The task is to estimate MI between x and y. We compare the performance ofMLMI(CV), KDE(CV), KNN(k) for k = 1, 5, 15, and EDGE; the approximation errorof an MI estimate Ib is measured by |Ib − I|. Figure 2 depicts the average approximation
error—MLMI, KDE, KNN(5), and EDGE perform well for the dataset (a), MLMI tends tooutperform the other estimators for the dataset (b), MLMI and KNN(5) show the best performance against the other methods for the dataset (c), and MLMI, EDGE, and KNN(15)perform well for the dataset (d). In the above simulation, KDE works moderately well forthe datasets (a)–(c), while it performs poorly for the dataset (d). This instability would
be due to devision by estimated densities, which tends to magnify the estimation error.
KNN seems work well for all four datasets if the value of k is chosen optimally; the bestvalue of k varies depending on the datasets and thus using a prefixed value of k is not
appropriate—k needs to be chosen adaptively using the data samples. However, there isno systematic model selection strategy for KNN and therefore KNN would be unreliablein practice. EDGE works well for the datasets (a) and (b), which posses high normality2 .
However, for the datasets (c) and (d) where normality of the target distributions is low, theEDGE method performs poorly. In contrast, MLMI with CV performs reasonably well for
all four datasets in a stable manner.
These experimental results show that MLMI nicely compensates for the weaknesses ofthe existing methods and therefore we conclude that MLMI should be regarded as a usefulalternative to the existing methods of MI estimation.

5. Conclusions
In this paper, we proposed a new method of estimating mutual information. The proposed
method, called MLMI, has several useful properties, e.g., it is a single-shot procedure,
density estimation is not involved, it is equipped with a cross-validation procedure for model
selection, and the unique global solution can be computed efficiently. We have provided arigorous convergence analysis of the proposed algorithm as well as numerical experiments
illustrating the usefulness of the proposed method.
2. Note that although the Edgeworth approximation is exact when the target distributions are completely
normal, the EDGE method still suffers from some estimation error since the cumulants are estimatedfrom samples.

14

Approximating Mutual Information by Maximum Likelihood Density Ratio Estimation

6

10

4
5
y

y

2
0

0

−2
−4

−1

0
x

−5

1

0
x

2

(b) Non-linear dependence with correlation

2

4

1

2

0

0

y

y

(a) Linear dependence

−2

−1

−2

−2
−0.5

0
x

−4

0.5

0

(c) Non-linear dependence without correlation

0.1

0.2

x

0.3

0.4

0.5

(d) Independence

Figure 1: Datasets used in experiments.

0.25
0.2
0.15
0.1
0.05
0
50

100
Number of samples

150

0.5
0.4
0.3
0.2
0.1

100
Number of samples

0.4

0.3

0.2

0.1

0
50

100
Number of samples

150(b) Non-linear dependence with correlation
Average MI approximation error

Average MI approximation error(a) Linear dependence

050

KNN (15)EDGE

KNN (1)
KNN (5)
Average MI approximation error

Average MI approximation error

MLMI
KDE

150

(c) Non-linear dependence without correlation

0.25
0.2
0.15
0.1
0.05
0
50

100
Number of samples

150(d) Independence

Figure 2: MI approximation error measured by |Ib− I| averaged over 100 trials as a functionof the sample size n. The symbol ‘◦’ on a line means that the corresponding
method is the best in terms of the average error or is judged to be comparable tothe best method by the t-test at the significance level 1%.
15

Suzuki, Sugiyama, Sese and Kanamori

Appendix: Sketch of Proof of Theorem 1 and Corollary 2
Proof of Theorem 1
For notational simplicity, we define the linear operators Qn , Q, Pn , P asPPnf (xi , yi )1≤i6=j≤n f (xi , yj )Qn f :=, Qf := Epx py f, Pn f := i=1, P f := Epxy f.
n(n − 1)n√
Step1: Proof of supϕ∈Fn |(Qn − Q)ϕ| = Op (1/ n). Hoeffding (1963) derived Bernstein’s
inequality for U -statistics which, in our context, is written as




x
x
√
P |(Qn − Q)f | > √≤ 2 exp −,(8)n(2xkf k∞ / n + 6Q(f − Q(f ))2 )for all x > 0 and n ≥ 2 (see also Proposition 2.3 of Arcones and Giné, 1993). By applying theabove inequality to the proof of Theorem 5.11 of van de Geer (2000) instead of Bernsten’s
inequality for the i.i.d. sum, we obtain a “double sum” version of Theorem 5.11 of van de√Geer (2000), i.e., exponential decay of the tail probability of n supf ∈F |(Qn −Q)f |, where F
is a class of uniformly bounded measurable functions and satisfies the polynomial bracketing
condition (6) as G M . This will be used √for obtaining√(9) later. By Assumption 1.2, we have
√
√
2
′
2
ξ
′
0
Fn ⊂ G . Now Q(ϕ − ϕ ) = Q( ϕ − ϕ ) ( ϕ + ϕ′ )2 ≤ 4ξ0 hQ (ϕ, ϕ′ )2 yieldspN[] (ǫ, G ξ0 , L2 (Q)) ≤ N[] (ǫ/2 ξ0 , G ξ0 , hQ ).
Thus the uniform convergence theorem (van de Geer, 2000, Theorem 5.11) with Assumption
1.3 implies
1
sup |(Qn − Q)ϕ| = Op √
.
(9)
n
ϕ∈Fn

f = 2ξ0 /ǫ0 and G (c) be
Let M
G

(c)

:=



2g
g + ḡ

g0≤ c0 .
g ∈ G , ḡ ∈ G , e
c0 ≤
ḡ
f
M

f
M

√
f
Step2: Proof of N[] (ǫ, G (c) , L2 (P )) ≤ N[] (ǫ/8 2c0 , G M , hQ )2 .
(9) yields√
inf Qn ϕ ≥ ǫ0 − Op (1/ n),
ϕ∈Fn

and
Q(S̄n ) → 1 for S̄n :=
inf Qn ϕ ≥ ǫ0 /2 .

ϕ∈Fn

On the event S̄n , all the elements of Gbn is uniformly bounded from above:
X
X
X
X
αl ≤ 2/ǫ0 .
αl ǫ0 /2 ⇒αl Qn (ϕl ) ≥αl ϕl ) =1 = Qn (lllf
Thus on the event S̄n , Gbn ⊂ GnM is always satisfied. We defineg0
Sn := S̄n ∩ e
c0 < ∗ < c0 .
gn

16

l

(10)

Approximating Mutual Information by Maximum Likelihood Density Ratio Estimation

The rest of the proof goes through an analogous line to Theorem 7.6 of van de Geer
(2000). Since Gbn is convex, from the definition of gbn , we obtain0 ≤ Pn log

2b
gn
2b
gn
g ∗ − gbn
2b
gn
≤ Pn
− 1 = (Pn − P )
−P n
.
∗
∗
∗
gbn + gn
gbn + gn
gbn + gn
gbn + gn∗

Since gn∗ maximizes P log g in the convex set Gbn , we haved
P log(αg + (1 − α)gn∗ )
≤ 0,
dα
α=0

This yields the inequalityPg − gn∗
≤ 0,
gn∗

(11)

∀g ∈ Gbn .

∀g ∈ Gbn .

Thus the second term of the far RHS of (11) is bounded from below as
P

1 (g ∗ − gbn )2gn∗ − gbn
1 g ∗ − gbn 1 (gn∗ − gbn )2≥ P ∗n∗= P n ∗
+ P ∗ ∗∗
gbn + gn
2gn
2 gn (gn + gbn )2 gn (gn + gbn )√ ∗ p 2 √ ∗ p 2
p
p
1 ( gn − gbn ) ( gn + gbn )
1
2 g0
∗ −
g
g
b
)
= P
≥
Q
(
.
n
n
2
gn∗ (bgn + gn∗ )2
gn∗

Combining (11) and (12), on the event g0 /gn∗ ≥ e
c0 , we have

p
2b
gn
1 p
.
e
c0 hQ (gn∗ , gbn )2 = e
c0 Q( gn∗ − gbn )2 ≤ (Pn − P )2
gbn + gn∗

(12)

(13)

gn
e ∈ G be g0 /ḡ ≤ c0 andThis indicates that it suffices to bound (Pn − P ) gbn2b+g∗ . Let ḡ, gn
g0 /e
g ≤ c0 , then
p
p


√
√
4g0 g 2 ( ḡ − ge)2 ( ḡ + ge)2
2g 2
2g
4g 2 (ḡ − ge)2
−
=Q
=P
P
g + ḡ g + ge
(g + ḡ)2 (g + ge)2
(g + ḡ)2 (g + ge)2p√
≤ 16c0 Q( ḡ − ge)2 .

Similarly, for g, g ′ ∈ G, on the event g0 /gn∗ < c0 , we have

2
p
2g ′
2g
√
− ′
≤ 16c0 Q( g − g ′ )2 .
P∗∗
g + gn g + gn

Therefore, for g, g ′ ∈ G and ḡ, ge ∈ G such that g0 /ḡ ≤ c0 and g0 /e
g ≤ c0 , the following issatisfied:
s 
2
√
2g ′
2g
≤ 4 2c0 (hQ (g, g ′ ) + hQ (ḡ, ge)).
− ′
P
(14)
g + ḡ g + gegn(c) ) → 1 holds by the assumptions (7) and (10). Since 2g is
Note that P ( gbn2b+g
∗ ∈ G
g+ḡ
n
increasing with respect to g and decreasing with respect to ḡ, (14) yields that the bracketing
number of G (c) is bounded by√f(15)
N[] (ǫ, G (c) , L2 (P )) ≤ N[] (ǫ/8 2c0 , G M , hQ )2 .

17

Suzuki, Sugiyama, Sese and Kanamori−1

Step3: Proof of P (hQ (bgn , gn∗ ) > δn ) → 0 where δn = O(n 2+γ ).
Now (15) gives an upper bound of the entropy integral of G (c) as follows:
Z δq
Z δq
√
1
1
f
(c)
√
, hQ ))dǫ.
log(N[] (ǫ, G , L2 (P )))dǫ ≤ √2 log(N[] (ǫ/8 2c0 , G M
n 0
n 0Note that on the event g0 /gn∗ < c0 , we have
√
√


 √
2


√
( g − gn∗ )2 ( g + gn∗ )22g
g − gn∗ 2
P
=P
−1 =P
g + gn∗g + gn∗(g + gn∗ )2
p 2√≤ 2c0 Q( g − gn∗ ) ≤ 4c0 hQ (g, gn∗ )2 .

(16)

Theorem 5.11 of van de Geer (2000) yields that there exist constants K ′ , C, C ′ , c such thatfor all δ ≥ δn where δn satisfies
Z q
K ′ δn2log(N[] (ǫ, G (c) , L2 (P )))dǫ,δn ≥ √(17)n 0
the following inequalities are satisfied:
P (hQ (bgn , gn∗ ) > δ)

!

2g
∗ 2
−1 ≥e
c0 hQ (g, gn )≤P
sup(Pn − P )g + gn∗∗ )>δg∈Gbn :hQ (g,gn
!


∞
X
√ 2s 2√2gn(Pn − P )−1 ≥e
c0 n2 δ ∩ Sn + P (Snc )
≤
P
sup
∗
g
+
g
∗ )≤2s+1 δng∈Gbn :hQ (g,gn
s=0
!
∞
X
√ 2s 2√n(Pn − P ) (f − 1) ≥ e
c0 n2 δ + P (Snc )
P
sup
≤
≤

s=0
∞
X
s=0



√
f ∈G (c) :kf −1kP ≤ c0 2s+2 δn22s δ 2nδ 2c
C exp −+ P (Sn ) ≤ c exp − 2 + P (Snc ).
C′
c

(18)

In the above equation, Snc denotes the complement of Sn and the second inequality follows
the application of the “peeling device” (see p.70 and Theorem 7.6 of van de Geer, 2000, for− 1
more details). Now we can take δn = O(n 2+γ ) because to ensure (17) it suffices to let δnsatisfyZ q√√
K ′ δnf
δn2 ≥ √
2 log(N[] (ǫ/8 2c0 , G M, hQ ))dǫ = O(δn1−γ/2 / n).
n 0
−

1

This yields that the far RHS of (18) converges to 0 and therefore hQ (bgn , gn∗ ) = Op (n 2+γ )holds.
Proof of Corollary 2
By the definition of gn∗ , P log gn∗ ≥ P log Qgn0g0 . On the other hand, since g0 maximizes ∗ gn
P log g in {g ∈ G | Qg = 1}, we have P log Qg≤ P log g0 . Combining these inequalities
∗
n
gives




g0 /Qn g0
g0∗= P log
Qn g0 Qgn ≤ log(Qn g0 Qgn∗ ).
(19)
0 ≤ P loggn∗ /Qgn∗
gn∗
18Approximating Mutual Information by Maximum Likelihood Density Ratio Estimation

A similar argument to Lemma 1.3 of van de Geer (2000) yields


g0
2
∗
∗
2hQ (gn /Qgn , g0 ) ≤ P log.
gn∗ /Qgn∗
Moreover we havelog(Qn g0 Qgn∗ ) ≤ Qn g0 Qgn∗ − 1 = (Qgn∗ )(Qn − Q)



gn∗
g0 −
Q(gn∗ )

.

Then, by (19), we have2h2Q (gn∗ /Qgn∗ , g0 ) ≤ (Qgn∗ )(Qn − Q) g0 −gn∗
Q(gn∗ )

.

f

We can easily check that gn∗ /Q(gn∗ ) is contained in G M . Thus
s
s
!2
!2

2
2

gn∗
gn∗
gn∗
gn∗
√
√
f
−
+
=
Q
.
Q g0 −
≤
4(
M
+η
)h
g
,
g
g
0
0
0 Q
0
Q(gn∗ )
Q(gn∗ )
Q(gn∗ )
Q(gn∗ )fSimilarly, for Ge := {g − g0 | g ∈ G M }, we havepe L2 (Q)) ≤ N[] (ǫ/2 2M
f, G,e hQ ).
N[] (ǫ, G,

Applying a similar argument to (18) with the double sum version of Theorem 5.11 of van deGeer (2000), we have


gn∗
− 1
hQ g0 ,
= Op (n 2+γ ).
∗Q(gn )
∗11gn= Op (n− 2 ), we have hQ (g0 , gn∗ ) =Since Q(gn∗ ) = 1 + Op (n− 2 ) yields hQ gn∗ , Q(g
∗)
n

Op (n

1
− 2+γ

) by the triangle inequality.

References
M. A. Arcones and E. Giné. Limit theorems for U -processes. The Annals of Probability, 21(3):1494–1542, 1993.
S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, Cambridge, 2004.
P. Comon. Independent component analysis, a new concept?
287–314, 1994.

Signal Processing, 36(3):T. M. Cover and J. A. Thomas. Elements of Information Theory. John Wiley & Sons, Inc.,
N. Y., 1991.
A. M. Fraser and H. L. Swinney.
Independent coordinates for strange attractors from mutual information. Physical Review A, 33(2):1134–1140, 1986. doi:
http://dx.doi.org/10.1103/PhysRevA.33.1134.
19

Suzuki, Sugiyama, Sese and Kanamori

S. Ghosal and A. W. van der Vaart. Entropies and rates of convergence for maximum
likelihood and Bayes estimation for mixtures of normal densities. Annals of Statistics,
29:1233–1263, 2001.
I. Guyon and A. Elisseeff. An introduction to variable feature selection. Journal of Machine
Learning Research, 3:1157–1182, 2003.
W. Härdle, M. Müller, S. Sperlich, and A. Werwatz. Nonparametric and Semiparametric
Models. Springer Series in Statistics. Springer, Berlin, 2004.
W. Hoeffding. Probability inequalities for sums of bounded random variables. Journal of
the American Statistical Association, 58(301):13–30, 1963.
M. M. Van Hulle. Edgeworth approximation of multivariate differential entropy. NeuralComputation, 17(9):1903–1910, 2005.
S. Khan, S. Bandyopadhyay, A. Ganguly, and S. Saigal. Relative performance of mutual
information estimation methods for quantifying the dependence among short and noisy
data. Physical Review E, 76:026209, 2007.
A. Kraskov, H. Stögbauer, and P. Grassberger. Estimating mutual information. PhysicalReview E, 69:066138, 2004.
X. Nguyen, M. J. Wainwright, and M. I. Jordan. Estimating divergence functions and thelikelihood ratio by penalized convex risk minimization. In Advances in Neural Information
Processing Systems 20, 2008.
G. Rockafellar. Convex Analysis. Princeton University Press, Princeton, 1970.
B. Schölkopf and A. J. Smola. Learning with Kernels. MIT Press, Cambridge, MA, 2002.
R. J. Serfling. Approximation Theorems of Mathematical Statistics. Wiley, New York, 1980.
B. W. Silverman. Density Estimation for Statistics and Data Analysis.
Hall/CRC, April 1986.

Chapman &M. Sugiyama, T. Suzuki, S. Nakajima, H. Kashima, P. von Bünau, and M. Kawanabe.
Direct importance estimation for covariate shift adaptation. Annals of the Institute of
Statistical Mathematics, 60(4), 2008. to appear.
T. Suzuki, M. Sugiyama, J. Sese, and T. Kanamori. A least-squares approach to mutual
information estimation with application in variable selection. In Proceedings of the 3rdWorkshop on New Challenges for Feature Selection in Data Mining and Knowledge Discovery (FSDM2008), 2008. to appear.
K. Torkkola. Feature extraction by non-parametric mutual information maximization. Journal of Machine Learning Research, 3:1415–1438, 2003.
S. van de Geer. Empirical Processes in M-Estimation. Cambridge University Press, 2000.

20