Mutual Information Approximation via
Maximum Likelihood Estimation of Density Ratio
Taiji Suzuki

Masashi Sugiyama

Toshiyuki Tanaka

Dept. of Mathematical Informatics,
The University of Tokyo
Email: s-taiji@stat.t.u-tokyo.ac.jp

Dept. of Computer Science,
Tokyo Institute of Technology.
Email: sugi@cs.titech.ac.jp

Dept. of Systems Science,
Kyoto University
Email: tt@i.kyoto-u.ac.jp

Abstract—We propose a new method of approximating mutual
information based on maximum likelihood estimation of a density
ratio function. The proposed method, Maximum Likelihood Mutual Information (MLMI), possesses useful properties, e.g., it does
not involve density estimation, the global optimal solution canbe efficiently computed, it has suitable convergence properties,
and model selection criteria are available. Numerical experiments
show that MLMI compares favorably with existing methods.

I. I NTRODUCTION
Mutual information (MI) between two random variables X
and Y∫∫
pXY (x,y)I(X, Y ) :=pXY (x, y) log pX(1)(x)pY (y) dxdy
plays a central role in information theory and statistics since itvanishes if and only if X and Y are independent. A great dealof effort has been made to estimate MI from samples. Such
MI estimators are useful in testing independence, and havebeen applied to various machine learning and signal processing
tasks such as feature selection and independent component
analysis.
In this paper, we propose a new MI estimator, Maximum
Likelihood Mutual Information (MLMI), that can overcome
the limitations of existing approaches. MLMI does not involve
density estimation, but directly learns the density ratiopXY (x,y)w(x, y) := pX(x)pY (y)(2)

via maximum likelihood (ML) estimation. Given a densityratio model g(x, y), MI can simply be approximated by
∑n
b
I(X,
Y ) := n1 i=1 log g(xi , yi ).

(3)

We prove consistency of the MLMI procedure and elucidate
its rate of convergence. Estimation of the density ratio isformulated in the ML framework as a convex optimization
problem. Therefore, the unique global optimal solution canbe obtained efficiently. Furthermore, cross validation and information criteria are available for model selection, so thatvalues of tuning parameters, such as the kernel width, can beadaptively determined in an objective manner. Numerical experiments show that MLMI compares favorably with existing
methods.

II. N EW MI E STIMATOR
In this section, we formulate the MI estimation problem
as a density ratio estimation problem and propose a new MI
estimation method.
Formulation: Let DX ⊂ RdX and DY ⊂ RdY be data
domains. Suppose we are given n independent and identically
distributed (i.i.d.) paired samples {(xi , yi ) ∈ DX × DY }ni=1drawn from a joint distribution with density pXY (x, y). Let
pX (x) and pY (y) denote the marginal densities of X and Y ,respectively. The goal is to estimate MI defined by (1).
A key factor of our approach is to avoid density estimation, since density estimation is harder than estimationof MI; instead, we estimate the density ratio w(x, y) defined by (2). An obtained density-ratio model g(x, y) isthen plugged into (3) to estimate MI. We model the density ratio by a linear model g(x, y) := α⊤ φ(x, y), where
α := (α1 , . . . , αb )⊤ are parameters to be learned fromsamples, where ⊤ denotes the transpose of a vector, and whereφ(x, y) := (φ1 (x, y), . . . , φb (x, y))⊤ are basis functions
such that φ(x, y) ≥ 0b for all (x, y) ∈ DX × DY . 0b
denotes the b-dimensional all-zero vector. The basis functions
φ(x, y) and their number b could be dependent on the samples{(xi , yi )}ni=1 , i.e., kernel models are also allowed.
Maximum Likelihood Estimation of Density Ratio: Given a density-ratio model g(x, y), one may estimate the joint density pXY (x, y) by pbXY (x, y) :=g(x, y)pX (x)pY (y), based on which we propose to estimate
the parameter α of g(x, y) from the samples in the ML
framework: Ignoring irrelevant constants, it is equivalent tomaximizing
(
)
∑n
b
J(α)
:= n1 i=1 log α⊤ φ(xi , yi ) .
This is our objective function, which is concave.
The density-ratio model g(x, y) should be non-negativeby definition. It is therefore natural to impose the constraintg(x, y) ≥ 0 for all (x, y) ∈ DX × DY . This can be achievedby imposing α ≥ 0b . In addition to non-negativity, g(x, y)should be properly normalized since pbXY (x, y) is a density:
∫∫
∑
1
⊤
1=
pbXY (x, y)dxdy ≈ n(n−1)
1≤i̸=j≤n α φ(xi , yj ),where we used the U -statistic for∑ obtaining ⊤the empirical1estimator. Let Sb := {α | n(n−1)
1≤i̸=j≤n α φ(xi , yj ) =1, α ≥ 0b }. Then our optimization criterion is summarized as

Let E denote expectation over random samples{(xi , yi )}ni=1 . The behavior of E[I(bg )] as a function ofbb := arg maxα∈Sb J(α).
α
the number of samples n is called the learning curve. Theb ⊤ φ(x,Let gb(x, y) := αThen our MI estimator is given by learning curve has been a quantity of interest in statistics∑y).
n
1(MLMI)bI(X, Y ) := n i=1 log gb(xi , yi ). We call this method since it can be used for deriving information criteria. Here,in order to avoid technical difficulties, we assume in the
Maximum Likelihood Mutual Information (MLMI).
analysis of the learning curve of MLMI that C⊥ is a linear
space. Then the learning curve can be expressed as follows.
III. R ATES OF C ONVERGENCE†)c
Theorem 3: E[I(bg )] = I(g ∗ ) − tr(RGR(RQR)
+ 2n
+
2n
Parametric Cases: We start with parametric cases where
3
−2
),
where
†
denotestheMoore-Penrose
pseudo-inverse.
O(nthe number of basis functions b is fixed. We derive the asympNon-Parametric Cases: Now we go on to elucidation ofb the convergence rate
totic distribution of learned parameter α,
thenon-parametric convergence rates.
of MI estimation, and the learning curve of MLMI. We assumeFirst,
we define some notations. The set of basis functions is
the same regularity conditions as those used in [11], which isdenotedby F := {φθ | θ ∈ Θ}, where Θ is a parameter/index
a common setup in developing asymptotic parametric theory.
set.
The
set of basis functions used for estimation with nFirst, we define some notations. Let J and S be the ‘ideal’samplesischaracterized by a subset of the parameter set
objective function and feasible set of MLMI, respectively:
Θ
⊆
Θ
and
denoted by Fn := {φθ | θ ∈ Θn } ⊂ F,n∫∫J(α) :=pXY (x, y) log(α⊤ φ(x, y))dxdy,
which can be stochastic. The set of finite linear combinations∫∫
S := {α |pX (x)pY (y)α⊤ φ(x, y)dxdy = 1, α ≥ 0}. of F with non-negative coefficients and its bounded subset are
denoted respectively by∑
Let α∗ := arg maxα∈S J(α) and g ∗ (x, y) := α∗ ⊤ φ(x, y).
G := { l αl φθl | αl ≥ 0, φθl ∈ F},
S is a convex polytope and the approximating cone C at α∗
G M := {g ∈ G | ∥g∥∞ ≤ M }.
is defined as C := {λ(α − α∗ ) | α ∈ S, λ ≥ 0}. Let ∇α
be a partial derivative operator with respect to α. We denote Their subsets used for estimation with n samples are denoted∇α log g|g=α∗⊤φ by ∇α log g ∗ ; similarly we use ∇α log gb. Let by G := {∑ α φ | α ≥ 0, φ ∈ F } ⊂ G. Let Gb and
l
θl
n
n
n
l l θl
∫∫
g
b
be
the
feasible
set
and
the
solution
of
MLMI:
∗
⊤
n
pXY (x, y)∇α log g (x, y)dxdy = 0},C⊥ := C ∩ {δ | δ∑∫∫1
Gbn := {g ∈ Gn | n(n−1)
G :=pXY (x, y)∇α log g ∗ (x, y)∇α⊤ log g ∗ (x, y)dxdy,1≤i̸=j≤n g(xi , yj ) = 1},∫∫][ 1 ∑n⊤∗Q := −pXY (x, y)∇α ∇α log g (x, y)dxdy,
gbn := arg maxg∈Gbn n i=1 log g(xi , yi ) .
∫∫∫c := n1pX (x)pX (x′ )pY (y)g ∗(x, y)g ∗(x′ , y)dxdx′ dy
For simplicity, we assume that the solution gbn is unique.
∫∫∫∫
′ ′
+ n2
pXY (x, y)pX(x′)pY (y ′)g ∗(x, y ′ )g ∗(x′ , y)dxdydxdy
Here we use the (generalized) Hellinger distance with
∫∫∫
1
4
′
∗
∗
′
′
respect
to pX pY as the error metric, since this allows us to+npX (x)pY (y )pY (y)g (x, y)g (x, y )dxdydy − n
∫∫
avoid
some
technical difficulties:
∗
′
′
ϕ :=
pXY (x, y)pY (y ′ ) gφ(x,y)∗ (x,y) g (x, y )dxdydy
(∫∫
)1/2
√
√
′
hQ (g, g ) :=pX (x)pY (y)( g − g ′ )2 dxdy,∫∫
′ φ(x,y) ∗′′+pXY (x, y)pX (x ) g∗ (x,y) g (x , y)dx dxdy.
where g and g ′ are non-negative measurable functions (notLet R be the orthogonal projection matrix onto the linear hull necessarily probability densities). We further make the folof C⊥ . Let N (µ, Σ) be the normal distribution with mean µ lowing assumptions.
and covariance matrix Σ. Then we have the following theorem. 1. On the support of pXY , there exists a constant η < ∞ such√b − α∗ ) converges in law to that the true density ratio w is upper-bounded as w(x, y) ≤ η.
Theorem 1: n(αarg minδ∈C⊥ ∥δ − Z∥Q + α∗ Z ′ , where ∥α∥2Q := α⊤ Qα, 2. All basis functions are non-negative, and there exist conand [where∈ R are random variables]) such ∫∫] Z ∈(Rb and[ Z ′ −1stants ϵ, ξ > 0 such that ∀φ ∈ F, ∥φ∥∞ ≤ ξ and
RQ RGRQ−1 R RϕZ
pX (x)pY (y)φ(x, y)dxdy ≥ ϵ.
that
∼
N
0
,
.
b+1
Z′
ϕ⊤ R
c
3. There exist constants 0√ < γ < 2 and K > 0 such that
We analyze the convergence rate of our MI estimator. Let
log N[] (ϵ, G M , hQ ) ≤ K( M /ϵ)γ , where N[] (ϵ, F, h) is the∫∫
ϵ-bracketing number of F with distance h [9].
I(g) :=pXY (x, y) log g(x, y)dxdy,∑
Assumption 3 ensures that the model is not so complicated.
n
1
b :=
I(g)
∗
i=1 log g(xi , yi ).
nGaussian mixture∫∫ models satisfy this condition [2]. Let gn :=pXY (x, y) log g(x, y)dxdy. Then we have
Let Op be the asymptotic order in probability. Then, we have arg maxg∈Gbn
the following theorems.
b g ) − I(w)| = |I(g ∗ ) − I(w)| + Op (n− 21 ).
Theorem 2: |I(bTheorem 4: If there exist c0 , c1 such that P{(xi ,yi )} (c0 ≤Corollary 1: If the parametric model contains the true
1
1
w
∗
−
b g ) − I(w)| = Op (n 2 ).
gn , w) = Op (n− 2+γ + hQ (gn∗ , w)).
∗ ≤ c1 ) → 1, then hQ (bdensity ratio, i.e., g = w, then |I(bgn
Theorem 5: If there exists δ > 0 such that g(x, y) ≥ δ for
The corollary means that MLMI retains optimality in terms
1b gn ) − I(w)| =of the order of convergence in n, since Op (n− 2 ) is the optimal ∀(x, y) ∈ DX × DY and1 ∀g ∈ Gbn , then |I(bconvergence rate in the parametric setup.
|I(gn∗ ) − I(w)| + Op (n− 2+γ ).

Corollary 2: If there exists N such that w ∈ Gn for∀n ≥ N , and if there exists δ > 0 such that g(x, y) ≥ δb gn ) − I(w)| =for ∀(x, y) ∈ DX × DY and ∀g ∈ Gbn , then |I(b
1
− 2+γ
Op (n
).
This corollary shows that the convergence rate of nonparametric MLMI is slightly slower than the parametric counterpart, but the non-parametric method requires a milder model
assumption for eliminating the modeling error. According
to [5], the above convergence rate achieves the optimal minimax rate under some setup. Thus the convergence property of
non-parametric MLMI would be optimal in the same sense.
IV. M ODEL S ELECTION
We have shown that MLMI possesses preferable convergence properties. However, the practical performance of
MLMI may strongly depend on the choice of the basis
functions φ(x, y). Here we show how model selection of
MLMI can be carried out.
Parametric Cases: For parametric models, we provide an
information criterion which can be used for model selection.
b be theAssume that C⊥ is a linear space for simplicity. Let Hbb and letinclusionwise minimal face of S which contains α,b |λ∈Cb⊥ be the linear space defined by Cb⊥ := {λ(α − α)b be the orthogonal projection onto Cb⊥ andb Let RR, α ∈ H}.
let
bG
b R(
b R
bQ
b R)
b † ),b g ) − 1 tr(R
Ib(IC) (bg ) := I(bnb := 1 ∑n ∇α log gb(xi , yi )∇⊤ log gb(xi , yi ),
G
α
i=1
n
∑n
1
⊤
b
Q := − n i=1 ∇α ∇α log gb(xi , yi ).
Then we have the following theorem:
3
Theorem 6: E[Ib(IC) (bg )] = E[I(bg )] + O(n− 2 ).
The above theorem shows that Ib(IC) is an asymptotic
unbiased estimator of E[I(bg )] up to O(n−1 ). On the otherb g ) is asymptotically
hand, the naive empirical estimator I(b− 12bunbiased only up to O(n ): E[I(bg )] = E[I(bg )] + O(n−1 ).
(IC)bThus I
would be a more accurate estimator of I(bg ) thanthe naive one.
For model selection, we prepare a set of model candidates(the basis functions φ(x, y) in the current setting) and choose
the one that has the largest value of Ib(IC) as the ‘best’ model.
Non-Parametric Cases: For non-parametric models, wecannot unfortunately utilize the theoretical results given in theprevious section since the coefficient of the convergence rate
is not explicitly given. Here, we propose to use numericalestimation via cross validation for model selection.
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 gbZk (x, y) is obtained using{Zj }j̸=k (i.e., without Zk ) and the score for the hold-out
samples Zk is computed as∑(K-CV)
IbZk= |Z1k | (x′ ,y′ )∈Zk log gbZk (x′ , y ′ ),where |Zk | denotes the number of sample pairs in the set Zk .
Repeat this procedure for k = 1, . . . , K and output its average∑K b(K-CV)
1
Ib(K-CV) = K.
k=1 IZk

Then we have the following theorem.
Theorem 7: The leave-one-out cross-validation scoreIb(n-CV) is an unbiased estimate of MI learned from (n − 1)
samples: E[Ib(n-CV) ] = E[I(bgn−1 )].
Given that E[I(bgn−1 )] ≈ E[I(bgn )], cross validation gives analmost unbiased estimate of MI: Similar almost unbiasedness
can be established for general K-fold cross validation, but thebias may be larger. For model selection, we compute Ib(K-CV)
for all model candidates and choose the one that maximizes
the cross-validation score.
Basis Function Design: A good model may be chosen by
an information criterion or cross validation, given that a set ofpromising model candidates is prepared. As model candidates,
we propose to use a Gaussian kernel model:ℓ∥ℓ∥φℓ (x, y) = exp(− ∥x−u
) exp(− ∥y−v
),
2σ 2
2σ 2
2

2

where {(uℓ , vℓ )}bℓ=1 are Gaussian centers; we choose thecenters randomly from {zi | zi = (xi , yi )}ni=1 .
pXY (x,y)
By definition, the density ratio pX(x)pY (y) tends to takelarge values if pXY (x, y) is large and pX (x)pY (y) is small;conversely, the ratio tends to be small (i.e., close to zero)
if pX (x)pY (y) is large and pXY (x, y) is small. When a
non-negative function is approximated by a Gaussian mixture
model in general, many kernels may be needed in the region
where the output of the target function is large; on the other
hand, only a small number of kernels would be enough in theregion where the output of the target function is close to zero.
On the basis of this idea, we propose a heuristic to allocate
many kernels in regions where pXY (x, y) is large, which canbe achieved by setting the Gaussian centers at {(xi , yi )}ni=1 .
Alternatively, we might locate Gaussian kernels on{(xi , yj )}ni,j=1 , which, however, requires n2 basis functions
and therefore is computationally prohibitive when n is large.
Our preliminary experiments showed that using n2 kernels
did not improve the performance, but significantly increased
the computation time. For this reason, in the experiments tofollow, we decided to use the above heuristic with the numberof basis functions fixed to b = min(200, n) and to choose theGaussian width σ by cross validation.
V. R ELATION TO E XISTING M ETHODSIn this section, we discuss the characteristics of existing andthe proposed approaches.
Legendre-Fenchel Duality of KL-Divergence: MI (or,more generally, KL-divergence for arbitrary two probability
densities) can be characterized by Legendre-Fenchel dualityof the convex function ‘− log’ [7]. More precisely, I(X, Y )can be characterized as the solution of the following concave
maximization problem [5]:∫I(X, Y ) = supg≥0 [ (−pX (x)pY (y)g(x, y)+ pXY (x, y) log g(x, y))dxdy + 1],
where the supremum is taken over all non-negative measurable
functions. In [5], an empirical version of this optimization

problem is solved by restricting the search space within a(Gaussian) reproducing kernel Hilbert space with a regularizer.
Our MLMI method is closely related to the method givenby [5]: Indeed, if the linear model assumption and the normalization constraint are imposed, and if the expectation isapproximated by the sample average, the above formulationis reduced to MLMI. However, our approach would be moreadvantageous in the following respects.
• Our ML formulation is equipped with built-in regularization
effects due to non-negativity and normalization constraints.
This allows us to avoid introducing an additional regularization
parameter and contributes to reducing the computational cost.
• A non-parametric convergence rate was investigated forgeneral KL divergence estimation in [5]. On the other hand,
we elucidated the non-parametric convergence rate specificallyin the context of MI estimation; this includes proper treatmentof the ‘decoupling’ effect in pX (x)pY (y). Furthermore, thestrong positivity of the true density ratio w is not required in
our proof and our result also covers a situation where the true
density ratio is not contained in the model.
• We derived an asymptotic distribution of the estimator in theparametric setup and elucidated the parametric convergence
rate of MI estimation and the asymptotic learning curve.
• We gave a practical heuristic for designing basis functions
for approximating the density ratio. This allows us to reducethe computational cost significantly.
Kernel Density Estimation (KDE) and Adaptive Histogram Methods: KDE is a non-parametric technique to estimate a density function p(x) from its i.i.d. samples {xi }ni=1 .
For the Gaussian kernel, KDE is expressed as2∑ni∥pb(x) = n(2πσ12 )d/2 i=1 exp(− ∥x−x).
2σ 2
The performance of KDE depends on the choice of thekernel width σ, which can be optimized by cross validation.
KDE-based estimation of MI can be performed using density estimates pbXY (x, y), pbX (x), and pbY (y) obtained from{(xi , yi )}ni=1 , {xi }ni=1 , and {yi }ni=1 , respectively as∑npbXY (xi ,yi )Ib = n1 i=1 log pbX(xi )bpY (yi ) .
However, density estimation itself is known to be a hard
problem, and division by estimated densities may increase
estimation errors. For this reason, the KDE-based approach
may not be reliable in practice.
Histogram-based estimators with data-dependent partition
would be more adaptive density estimation schemes. In thecontext of KL divergence estimation, consistency propertiesof histogram-based methods, which could be regarded as
pXY (x,y)
implicitly estimating the ratio pX(x)pY (y) , have been investigated in [10], [8]. MI estimation following this line has
been explored in [1]. However, such histogram-based methods may seriously suffer from the curse of dimensionality,
and are therefore not reliable in high-dimensional problems.
Furthermore, the convergence rate seems to be unexplored yet.
K-Nearest Neighbor (KNN) Based Method: If estimatesof the entropies H are obtained, MI can be estimated via the

MI-entropy identity as
I(X, Y ) =H(X) + H(Y ) − H(X, Y ),∫H(X) := − pX (x) log pX (x)dx.
In [4], an entropy estimator that utilizes the KNN distancehas been developed. Let us define the norm of z = (x, y) by∥z∥z := max{∥x∥, ∥y∥}, where ∥ · ∥ denotes the Euclidean
norm. Let Nk (i) be the set of KNN samples of (xi , yi ) withrespect to the norm ∥ · ∥z , and letϵX (i) := max{∥xi − xi′ ∥ | (xi′ , yi′ ) ∈ Nk (i)},nX (i) := |{zi′ | ∥xi − xi′ ∥ ≤ ϵX (i)}|,ϵY (i) := max{∥yi − yi′ ∥ | (xi′ , yi′ ) ∈ Nk (i)},nY (i) := |{zi′ | ∥yi − yi′ ∥ ≤ ϵY (i)}|.
Then the KNN-based MI estimator is given by∑n
Ib = ψ(k) + ψ(n) − k1 − n1 i=1 [ψ(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 errors of individual entropy estimates.
A practical drawback of the KNN-based approach is that theestimation accuracy depends on the value of k and there seems
no systematic strategy to choose the value of k appropriately.
Recently a KL divergence estimator utilizing KNN density
estimation is proposed in [6] and its consistency has been
investigated. A notable property of this estimator is thatconsistency of density estimation is not necessary to establish
consistency of the KL divergence estimator. However, the rateof convergence seems to be still an open research issue.
Edgeworth Expansion (EDGE) Based Method: In [3], an
entropy estimator based on the Edgeworth expansion was proposed. The basic idea is to approximate the entropy by that of
the normal distribution and additional higher-order correction
terms. More specifically, for a d-dimensional distribution, anb of the entropy H is given by
estimator H2
κ2i,i,j ∑dκ2i,j,kb = H ∗ −∑d κi,i,i −∑dH
i=1 12i,j=1,i̸=j 4 − i,j,k=1,i<j<k 72 ,where H ∗ is the entropy of the normal distribution withcovariance matrix equal to the target distribution and κi,j,k isthe standardized third cumulant of the target distribution. An
estimate of MI can be obtained via the MI-entropy identity. In
practice, all the cumulants should be estimated from samples.
If the underlying distribution is close to normal, the above
approximation is accurate and the EDGE-based method works
well. However, if the distributions are far from the normal distribution, the approximation error becomes large and therefore
the EDGE-based method is no longer reliable.
VI. N UMERICAL E XPERIMENTS
In this section, we experimentally investigate the performance of the proposed and existing MI estimators. The task
is to estimate MI between X ∈ R and Y ∈ R. We used thefollowing four datasets (see Figure 1):

6

10

4
5
y

y

2
0

0

−2
−4

−1

0
x

−5

1

2

4

1

2

0

0

−1

0
x

−4

0.5

0

(c) Non-linear dependence without
correlation

Fig. 1.

0.2
0.15
0.1
0.05

100
Number of samples

150(a) Linear dependence

Average MI approximation error

0.4
0.3
0.2
0.1

100
Number of samples

x

0.3

0.4

0.5

VII. C ONCLUSIONS

0.4

0.3

0.2

0.1

0
50

100
Number of samples

150(b) Non-linear dependence with correlation

0.5

0
50

0.2

KNN (15)EDGE

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

0.25

0
50

0.1

(d) Independence

Datasets used in experiments.
MLMI
KDE

Average MI approximation error

2

−2

−2
−0.5

Average MI approximation error

0x(b) Non-linear dependence with correlationyy(a) Linear dependence

−2

150

(c) Non-linear dependence without
correlation

0.25
0.2
0.15
0.1
0.05
0
50

the dataset (d). KDE worked moderately well on the datasets(a)–(c), while it performed poorly on the dataset (d). Thisinstability would be ascribed to division by estimated densities,
which tends to magnify estimation errors. KNN seems to work
well on all four datasets if the value of k is chosen optimally.
As already mentioned, however, there is no systematic model
selection strategy for KNN, so that KNN would be unreliablein practice. EDGE worked well on the datasets (a), (b), and(d), which possess high normality. However, for the dataset(c), where normality of the target distributions is low, theEDGE method performed poorly. In contrast, MLMI withcross validation performed reasonably well for all four datasetsin a stable manner.
These experimental results have shown that MLMI nicelycompensates for the weaknesses of the existing methods, andwe therefore conclude that MLMI should be regarded as a
useful alternative to the existing methods of MI estimation.

100
Number of samples

150(d) Independence

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

(a) Linear dependence: Y has a linear dependence on X as
X ∼ N (0, 0.5) and Y |X ∼ N (3X, 1).
(b) Non-linear dependence with correlation: Y has a quadratic
dependence on X as X ∼ N (0, 1) and Y |X ∼ N (X 2 , 1).
(c) Non-linear dependence without correlation: Y has a latticestructured dependence on X as X ∼ U(− 21 , 12 ), Y |X ∼N (0, 13 ) if X ≤ | 16 |, and Y |X ∼ 12 (N (1, 13 ) + N (−1, 13 ))otherwise, where U(a, b) is the uniform density on (a, b).
(d) Independence: X and Y are independent of each other as
X ∼ U(0, 12 ) and Y |X ∼ N (0, 1).
We compared MLMI, KDE, KNN (k = 1, 5, 15), and
EDGE. Figure 2 depicts the average approximation errors—
MLMI, KDE, KNN with k = 5, and EDGE performed well on
the dataset (a), MLMI tends to outperform the other methods
on the dataset (b), MLMI and KNN with k = 5 showed thebest performance against the other methods on the dataset (c),and MLMI, EDGE, and KNN with k = 15 performed well on

We have 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 crossvalidation procedure for model selection, and the unique global
solution can be computed efficiently. We have provided arigorous convergence analysis of the proposed algorithm aswell as numerical experiments illustrating the usefulness ofthe proposed method.
R
A MATLAB⃝
implementation of MLMI is available from ‘http://sugiyama-www.cs.titech.ac.jp/˜sugi/software/
MLMI/index.html’. TS acknowledges support by MEXTGCOE Program (G05). MS acknowledges supports by MEXT
(2068000), JFE 21st Century Foundation, AOARD, and SCAT.
R EFERENCES[1] G. A. Darbellay and I. Vajda. Estimation of the information by anadaptive partitioning of the observation space. IEEE Transactions on
Information Theory, 45(4):1315–1321, 1999.
[2] S. Ghosal and A. W. van der Vaart. Entropies and rates of convergencefor maximum likelihood and Bayes estimation for mixtures of normal
densities. Annals of Statistics, 29:1233–1263, 2001.
[3] M. M. Van Hulle. Edgeworth approximation of multivariate differential
entropy. Neural Computation, 17(9):1903–1910, 2005.
[4] A. Kraskov, H. Stögbauer, and P. Grassberger. Estimating mutual
information. Physical Review E, 69:066138, 2004.
[5] X. Nguyen, M. J. Wainwright, and M. I. Jordan. Nonparametricestimation of the likelihood ratio and divergence functionals. IEEE
International Symposium on Information Theory, pages 2016–2020,
Nice, France, 2007.
[6] F. Pérez-Cruz. Kullback-Leibler divergence estimation of continuous
distributions. IEEE International Symposium on Information Theory,
pages 1666–1670, Toronto, Canada, 2008.
[7] G. Rockafellar. Convex Analysis. Princeton University Press, 1970.
[8] J. Silva and S. Narayanan. Universal consistency of data-driven
partitions for divergence estimation. IEEE International Symposium on
Information Theory, pages 2021–2025, Nice, France, 2007.
[9] A. W. van der Vaart and J. A. Wellner. Weak Convergence and Empirical
Processes. With Applications to Statistics. Springer, New York, 1996.
[10] Q. Wang, S. R. Kulkarmi, and S. Verdú. Divergence estimation ofcontiunous distributions based on data-dependent partitions. IEEE
Transactions on Information Theory, 51(9):3064–3074, 2005.
[11] H. White. Maximum likelihood estimation of misspecified models.
Econometrica, 50(1):1–25, 1982.