Neural Computation. vol.25, no.10, pp.2734–2775, 2013.

1

Density-Diﬀerence Estimation∗
Masashi Sugiyama
Tokyo Institute of Technology, Japan.
sugi@cs.titech.ac.jp
http://sugiyama-www.cs.titech.ac.jp/˜sugi
Takafumi Kanamori
Nagoya University, Japan.
kanamori@is.nagoya-u.ac.jp
Taiji Suzuki
The University of Tokyo, Japan.
s-taiji@stat.t.u-tokyo.ac.jp
Marthinus Christoﬀel du Plessis
Tokyo Institute of Technology, Japan.
christo@sg.cs.titech.ac.jp
Song LiuTokyo Institute of Technology, Japan.
song@sg.cs.titech.ac.jp
Ichiro Takeuchi
Nagoya Institute of Technology, Japan.
takeuchi.ichiro@nitech.ac.jp

Abstract
We address the problem of estimating the diﬀerence between two probability densities. A naive approach is a two-step procedure of ﬁrst estimating two densities
separately and then computing their diﬀerence. However, such a two-step proceduredoes not necessarily work well because the ﬁrst step is performed without regard
to the second step and thus a small estimation error incurred in the ﬁrst stage cancause a big error in the second stage. In this paper, we propose a single-shot procedure for directly estimating the density diﬀerence without separately estimating
two densities. We derive a non-parametric ﬁnite-sample error bound for the proposed single-shot density-diﬀerence estimator and show that it achieves the optimal
convergence rate. We then show how the proposed density-diﬀerence estimator canbe utilized in L2 -distance approximation. Finally, we experimentally demonstrate
the usefulness of the proposed method in robust distribution comparison such asclass-prior estimation and change-point detection.
Keywords
density diﬀerence, L2 -distance, robustness, Kullback-Leibler divergence, kernel density estimation.
∗

An earlier version of this paper is presented at the Neural Information Processing Systems (NIPS)
conference, 2012 (Sugiyama et al., 2012c).

Density-Diﬀerence Estimation

1

2

Introduction

When estimating a quantity consisting of two elements, a two-stage approach of ﬁrst
estimating the two elements separately and then approximating the target quantity basedon the estimates of the two elements often performs poorly, because the ﬁrst stage iscarried out without regard to the second stage and thus a small estimation error incurredin the ﬁrst stage can cause a big error in the second stage. To cope with this problem,
it would be more appropriate to directly estimate the target quantity in a single-shotprocess without separately estimating the two elements.
A seminal example that follows this general idea is pattern recognition by the support
vector machine 1 (Boser et al., 1992; Cortes & Vapnik, 1995; Vapnik, 1998): Instead of
separately estimating two probability distributions of positive and negative patterns, thesupport vector machine directly learns the boundary between the positive and negativeclasses that is suﬃcient for pattern recognition. More recently, a problem of estimating
the ratio of two probability densities was tackled in a similar fashion (Qin, 1998; Sugiyama
et al., 2008; Gretton et al., 2009; Kanamori et al., 2009; Nguyen et al., 2010; Kanamori
et al., 2012; Sugiyama et al., 2012b; Sugiyama et al., 2012a): The ratio of two probability densities is directly estimated without going through separate estimation of the two
probability densities.
In this paper, we further explore this line of research, and propose a method for directly estimating the diﬀerence between two probability densities in a single-shot process.
Density ratios and density diﬀerences can both be used for comparing probability densities
via approximation of divergences such as the Kullback-Leibler (KL) divergence (Kullback& Leibler, 1951) and the L2 -distance. A divergence estimator can be used for solving various machine learning tasks including class-balance estimation under class-prior change(Saerens et al., 2002; du Plessis & Sugiyama, 2012), image segmentation and registration (Liu et al., 2010; Atif et al., 2003), target object detection and recognition (Gray& Principe, 2010; Yamanaka et al., 2013b), feature selection and extraction (Torkkola,
2003; Suzuki & Sugiyama, 2013), and change-point detection in time series (Kawahara &Sugiyama, 2012; Liu et al., 2013; Yamanaka et al., 2013a). In this divergence approximation scenario, density diﬀerences are more advantageous than density ratios in several
aspects. For example, density ratios can be unbounded even for simple cases (Cortes
et al., 2010; Yamada et al., 2013), whereas density diﬀerences are always bounded aslong as both densities are bounded. Thus, density diﬀerences are expected to be learned
more easily than density ratios. Also, density ratios are asymmetric and thus the “direction” needs to be determined by a user, whereas density diﬀerences are symmetric andthus there is no need to think about the direction. These are our primal motivations todevelop a density-diﬀerence estimator.
1

More precisely, Vapnik (1998) said as follows: If you possess a restricted amount of information
for solving some problem, try to solve the problem directly and never solve a more general problem as
an intermediate step. It is possible that the available information is suﬃcient for a direct solution but
is insuﬃcient for solving a more general intermediate problem. Two-stage density-diﬀerence estimation
corresponds to solving a more general problem of separate density estimation in the ﬁrst stage.

Density-Diﬀerence Estimation

3

Note that density ratios have their own applications beyond divergence approximation to which density diﬀerences may not be applied, such as importance sampling andconditional probability estimation (Sugiyama et al., 2012a). On the other hand, density
diﬀerences also have their own unique applications to which density ratios may not beapplied, such as the estimation of highest density-diﬀerence regions in ﬂow cytometricdata analysis (Duong et al., 2009) and unsupervised labeling (du Plessis, 2013). Thisimplies that, for density ratios and density diﬀerences, neither of them includes the otherin terms of ranges of applications. This is our additional motivation to pursue a practicalalgorithm for density-diﬀerence estimation.
For this density-diﬀerence estimation problem, we propose a single-shot method, called
the least-squares density-diﬀerence (LSDD) estimator, that directly estimates the density diﬀerence without separately estimating two densities. LSDD is derived within theframework of kernel regularized least-squares estimation, and its solution can be computed analytically in a computationally eﬃcient and stable manner. Furthermore, LSDD
is equipped with cross-validation, and thus all tuning parameters such as the kernel width
and the regularization parameter can be systematically and objectively optimized. We
derive a ﬁnite-sample error bound for the LSDD estimator in a non-parametric setup andshow that it achieves the optimal convergence rate.
We also apply LSDD to L2 -distance estimation and show that it is more accurate thanthe diﬀerence of KDEs, which tends to severely underestimate the L2 -distance (Anderson
et al., 1994). Compared with the KL divergence, the L2 -distance is more robust against
outliers (Basu et al., 1998; Scott, 2001; Besbeas & Morgan, 2004). We experimentallydemonstrate the usefulness of LSDD in robust distribution comparison such as semisupervised class-prior estimation and unsupervised change detection.
The rest of this paper is structured as follows. In Section 2, we derive the LSDD
method and investigate its theoretical properties. In Section 3, we show how LSDD
can be utilized for L2 -distance approximation. In Section 4, we illustrate the numericalbehavior of LSDD. Finally, we conclude in Section 5.

2

Density-Diﬀerence Estimation

In this section, we propose a single-shot method for estimating the diﬀerence between two
probability densities from samples, and analyze its theoretical properties.

2.1

Problem Formulation and Naive Approach

First, we formulate the problem of density-diﬀerence estimation.
Suppose that we are given two sets of independent and identically distributed samples′
X := {xi }ni=1 and X ′ := {x′i′ }ni′ =1 from probability distributions on Rd with densities p(x)Density-Diﬀerence Estimation4

and p′ (x), respectively:
i.i.d.

X := {xi }ni=1 ∼ p(x),′X ′ := {x′i′ }ni′ =1 ∼ p′ (x).
i.i.d.

Our goal is to estimate the diﬀerence f (x) between p(x) and p′ (x) from the samples X
and X ′ :f (x) := p(x) − p′ (x).
A naive approach to density-diﬀerence estimation is to use kernel density estimators(KDEs) (Silverman, 1986). For Gaussian kernels, the KDE-based density-diﬀerence estimator is given by
fe(x) := pb(x) − pb′ (x),
where
(
)
n
∑
1
∥x − xi ∥2pb(x) :=exp −,n(2πσ 2 )d/2 i=1
2σ 2
(
)
n′
∑
∥x − x′i′ ∥2
1
′
exp −
.
pb (x) := ′n (2πσ ′2 )d/2 i′ =12σ ′2The Gaussian widths σ and σ ′ may be determined based on cross-validation (Härdle et al.,2004).
However, we argue that the KDE-based density-diﬀerence estimator is not the best
approach because of its two-step nature: Small estimation error incurred in each density
estimate can cause a big error in the ﬁnal density-diﬀerence estimate. More intuitively,
good density estimators tend to be smooth and thus a density-diﬀerence estimator obtained from such smooth density estimators tends to be over-smoothed (Hall & Wand,
1988; Anderson et al., 1994, see also numerical experiments in Section 4.1.1).
To overcome this weakness, we give a single-shot procedure of directly estimating thedensity diﬀerence f (x) without separately estimating the densities p(x) and p′ (x).

2.2

Least-Squares Density-Diﬀerence Estimation

In our proposed approach, we ﬁt a density-diﬀerence model g(x) to the true densitydiﬀerence function f (x) under the squared loss2 :∫ ()2
argming(x) − f (x) dx.
(1)
g

2

Hall and Wand (1988) used a leave-one-out variant of this criterion for jointly determining the
bandwidths of two KDEs. See Section 4 for its numerical behavior.

Density-Diﬀerence Estimation

5

We use the following linear-in-parameter model as g(x):g(x) =b∑θℓ ψℓ (x) = θ ⊤ ψ(x),

(2)

ℓ=1

where b denotes the number of basis functions, ψ(x) = (ψ1 (x), . . . , ψb (x))⊤ is a bdimensional basis function vector, θ = (θ1 , . . . , θb )⊤ is a b-dimensional parameter vector,and ⊤ denotes the transpose. In practice, we use the following non-parametric Gaussian
kernel model as g(x):()∥x − cℓ ∥2g(x) =
θℓ exp −
,
2
2σ
ℓ=1
′
n+n
∑

(3)

where (c1 , . . . , cn , cn+1 , . . . , cn+n′ ) := (x1 , . . . , xn , x′1 , . . . , x′n′ ) are Gaussian kernel centers.
If n + n′ is large, we may use only a subset of {x1 , . . . , xn , x′1 , . . . , x′n′ } as Gaussian kernel
centers.
For the model (2), the optimal parameter θ ∗ is given by∫ ()2∗
θ := argmin
g(x) − f (x) dx
θ
]
[∫
∫
2
= argming(x) dx − 2 g(x)f (x)dx
θ[]= argmin θ ⊤ Hθ − 2h⊤ θ
=H

θ
−1

h,

where H is the b × b matrix and h is the b-dimensional vector deﬁned as∫
H := ψ(x)ψ(x)⊤ dx,∫∫
h := ψ(x)p(x)dx − ψ(x′ )p′ (x′ )dx′ .
Note that, for the Gaussian kernel model (3), the integral in H can be computed analytically as
)
(
)
(
∫
∥x − cℓ′ ∥2∥x − cℓ ∥2
exp −dx
Hℓ,ℓ′ = exp −
2σ 2
2σ 2
(
)
∥cℓ − cℓ′ ∥22 d/2= (πσ ) exp −,4σ 2
where d denotes the dimensionality of x. This is a part of the reason why we chose
the Gaussian kernel model in practice. Another reason for this choice is its theoreticalsuperiority, as discussed in Section 2.3.2.

Density-Diﬀerence Estimation

6

Replacing the expectations in h by empirical estimators and adding an ℓ2 -regularizer
to the objective function, we arrive at the following optimization problem:[]
b := argmin θ ⊤ Hθ − 2hb ⊤ θ + λθ ⊤ θ ,θ(4)θ

b is the b-dimensional vector deﬁned
where λ (≥ 0) is the regularization parameter and h
as
′

n
n
1∑
1 ∑
b
h=
ψ(xi ) − ′ψ(x′i′ ).
n i=1n i′ =1

Taking the derivative of the objective function in Eq.(4) and equating it to zero, we canb analytically as
obtain the solution θb = (H + λI b )−1 h,bθwhere I b denotes the b-dimensional identity matrix.
Finally, a density-diﬀerence estimator fb(x) is given asb⊤ ψ(x).
fb(x) = θ(5)

We call this the least-squares density-diﬀerence (LSDD) estimator.

2.3

Theoretical Analysis

Here, we theoretically investigate the behavior of the LSDD estimator.
2.3.1

Parametric ConvergenceFirst, we consider a linear parametric setup where basis functions in our density-diﬀerence
model (2) are ﬁxed.
√√
Suppose that n/(n + n′ ) converges to η ∈ [0, 1],
and
let
λ
=
o(
1/n,
1/n′ ). Then
√
∗
nn′ b
the central limit theorem (Rao, 1965) asserts that
′ (θ − θ ) converges in law to then+nnormal distribution with mean 0 and covariance matrix
H −1 ((1 − η)V p + ηV p′ )H −1 ,where V p denotes the covariance matrix of ψ(x) under the probability density p(x):
∫
(
)(
)⊤
V p :=
ψ(x) − ψ p ψ(x) − ψ p p(x)dx,(6)and ψ p denotes the expectation of ψ(x) under the probability density p(x):∫ψ p := ψ(x)p(x)dx.
This√ result implies that the LSDD estimator has asymptotic normality with asymptotic
order 1/n + 1/n′ , which is the optimal convergence rate in the parametric setup.

Density-Diﬀerence Estimation

2.3.2

7

Non-Parametric Error Bound

Next, we consider a non-parametric setup where a density-diﬀerence function is learned
in a Gaussian reproducing kernel Hilbert space (RKHS) (Aronszajn, 1950).
Let Hγ be the Gaussian RKHS with width γ:)(∥x − x′ ∥2′.
kγ (x, x ) = exp −γ2
Let us consider a slightly modiﬁed LSDD estimator that is more suitable for nonparametric error analysis3 : For n′ = n,
[
( n
)
]
n
∑
∑
1
1
fb := argmin ∥g∥2L2 − 2g(xi ) −g(x′i′ ) + λ∥g∥2Hγ ,n i=1n i′ =1g∈Hγ
where ∥ · ∥L2 denotes the L2 -norm and ∥ · ∥Hγ denotes the norm in RKHS Hγ .
Then we can prove that, for all ρ, ρ′ > 0, there exists a constant K such that, for allτ ≥ 1 and n ≥ 1, the non-parametric LSDD estimator with appropriate choice of λ andγ satisﬁes4()
2α′∥fb − f ∥2L2 + λ∥fb∥2Hγ ≤ K n− 2α+d +ρ + τ n−1+ρ
(7)with probability not less than 1−4e−τ . Here, d denotes the dimensionality of input vector
x, and α ≥ 0 denotes the regularity of Besov space to which the true density-diﬀerencefunction f belongs (smaller/larger α means f is “less/more complex”; see Appendix A for
2α
its precise deﬁnition). Because n− 2α+d is the optimal learning rate in this setup (Eberts& Steinwart, 2011), the above result shows that the non-parametric LSDD estimatorachieves the optimal convergence rate.
It is known that, if the naive KDE with a Gaussian kernel is used for estimating aprobability density with regularity α > 2, the optimal learning rate cannot be achieved(Farrell, 1972; Silverman, 1986). To achieve the optimal rate by KDE, we should choosea kernel function speciﬁcally tailored to each regularity α (Parzen, 1962). However, sucha kernel function is not non-negative and it is diﬃcult to implement it in practice. On
the other hand, our LSDD estimator can always achieve the optimal learning rate for a
Gaussian kernel without regard to regularity α.

2.4Model Selection by Cross-Validation

The above theoretical analyses showed the superiority of LSDD in terms of the convergence
rates. However, the practical performance of LSDD depends on the choice of models (i.e.,3More speciﬁcally, the regularizer is replaced from the squared ℓ2 -norm of parameters to the squared
RKHS-norm of a learned function, which is necessary to establish consistency. Nevertheless, we use thesquared ℓ2 -norm of parameters in experiments because it is simpler and performs well in practice.
4
Because our theoretical result is highly technical, we only describe a rough idea here. More precise statement of the result and its complete proof are provided in Appendix A, where we utilize themathematical technique developed in Eberts and Steinwart (2011) for a regression problem.

Density-Diﬀerence Estimation8

the kernel width σ and the regularization parameter λ). Here, we show that the model
can be optimized by cross-validation (CV).
′More speciﬁcally, we ﬁrst divide the samples X = {xi }ni=1 and X ′ = {x′i′ }ni′ =1 into T
disjoint subsets {Xt }Tt=1 and {Xt′ }Tt=1 , respectively. Then we obtain a density-diﬀerence
estimate fbt (x) from X \Xt and X ′ \Xt′ (i.e., all samples without Xt and Xt′ ), and computeits hold-out error for Xt and Xt′ as∫
2 ∑b ′2 ∑b(t)
ft (x) + ′CV := fbt (x)2 dx −ft (x ),
|Xt | x∈X
|Xt | ′ ′x ∈Xttwhere |X | denotes the number of elements in the set X . We repeat this hold-out validation
procedure for t = 1, . . . , T , and compute the average hold-out error as1∑
CV(t) .
T t=1
T

CV :=

Finally, we choose the model that minimizes CV.
R
A MATLAB⃝
implementation of LSDD is available from‘http://sugiyama-www.cs.titech.ac.jp/~sugi/software/LSDD/’.

3

L2-Distance Estimation by LSDD

In this section, we consider the problem of approximating the L2 -distance between p(x)and p′ (x),
∫ (
)2
2
′
L (p, p ) :=p(x) − p′ (x) dx,(8)′from samples X := {xi }ni=1 and X ′ := {x′i′ }ni′ =1 drawn independently from the probability
distributions with densities p(x) and p′ (x), respectively.

3.1

Basic Forms

For an equivalent expression2′

L (p, p ) =∫∫f (x)p(x)dx −f (x′ )p′ (x′ )dx′ ,if we replace f (x) with an LSDD estimator fb(x) and approximate the expectations by
empirical averages, the following L2 -distance estimator can be obtained:b ⊤ θ.
b
L2 (p, p′ ) ≈ h(9)Density-Diﬀerence Estimation

9

Similarly, for another expression2′

L (p, p ) =∫f (x)2 dx,

replacing f (x) with an LSDD estimator fb(x) gives another L2 -distance estimator:b⊤ H θ.
b
L2 (p, p′ ) ≈ θ

3.2

(10)

Reduction of Bias Caused by Regularization

Eq.(9) and Eq.(10) themselves give approximations to L2 (p, p′ ). Nevertheless, we argue
that the use of their combination, deﬁned byb ⊤θb−θb⊤ H θ,bb2 (X , X ′ ) := 2hL(11)

is more sensible.
To explain the reason, let us consider a generalized L2 -distance estimator of the following form:b ⊤θ
b + (1 − β)θ
b⊤ H θ,bβh(12)where β is a real scalar. If the regularization parameter λ (≥ 0) is small, then Eq.(12)can be expressed asb + op (λ),
b − λ(2 − β)hb ⊤ H −2 hb ⊤θb + (1 − β)θ
b⊤ H θb=hb ⊤ H −1 h
βh

(13)

where op denotes the probabilistic order (the derivation of Eq.(13) is given in Appendix B).
Thus, up to Op (λ), the bias introduced by regularization (i.e., the second term in the righthand side of Eq.(13) that depends on λ) can be eliminated if β = 2, which yields Eq.(11).
Note that, if no regularization is imposed (i.e., λ = 0), both Eq.(9) and Eq.(10) yieldb the ﬁrst term in the right-hand side of Eq.(13).
b ⊤ H −1 h,hEq.(11) is actually equivalent to the negative of the optimal objective value of the
LSDD optimization problem without regularization (i.e., Eq.(4) with λ = 0). This canbe naturally interpreted through a lower bound of L2 (p, p′ ) obtained by Legendre-Fenchel
convex duality (Rockafellar, 1970):
[ (∫
) ∫
]
∫
2
′
′
2
L (p, p ) = sup 2g(x)p(x)dx − g(x)p (x)dx − g(x) dx ,gwhere the supremum is attained at g = f . If the expectations are replaced by empirical estimators and the linear-in-parameter model (2) is used as g, the above optimization problem is reduced to the LSDD objective function without regularization (see Eq.(4)). Thus,
LSDD corresponds to approximately maximizing the above lower bound and Eq.(11) is
its maximum value.
Through eigenvalue decomposition of H, we can show thatb ⊤θb−θb⊤ H θb≥hb ⊤θb≥θb⊤ H θ.
b2hThus, our approximator (11) is not less than the plain approximators (9) and (10).

Density-Diﬀerence Estimation

3.3

10

Further Bias Correction⊤b H −1 h,b the ﬁrst term in Eq.(13), is an essential part of the L2 -distance estimatorh(11). However, it is actually a slightly biased estimator of the target quantity h⊤ H −1 h(= θ ∗ ⊤ Hθ ∗ = h⊤ θ ∗ ):
(
(
))
⊤
1
1
⊤
−1
−1
−1
b = h H h + tr Hb H h],(14)E[hV p + ′ V p′nn′where E denotes the expectation over all samples X = {xi }ni=1 and X ′ = {x′i′ }ni′ =1 , andV p and V p′ are deﬁned by Eq.(6) (its derivation is given in Appendix C).
The second term in the right-hand side of Eq.(14) is an estimation bias that is generally
non-zero. Thus, based on Eq.(14), we can construct a bias-corrected L2 -distance estimator
as
(
(
))
⊤
⊤
1b
1 b
−1
2
′
b
b
b
b
e
L (X , X ) := 2h θ − θ H θ − tr H
V p + ′ V p′
,
(15)
n
n
where Vb p is an empirical estimator of covariance matrix V p ,
n
)(
)⊤
1 ∑(
b
b
b
ψ(xi ) − ψ p ψ(xi ) − ψ p ,V p :=
n i=1

b is an empirical estimator of the expectation ψ :
and ψ
p
p
n
1∑
b
ψ(xi ).
ψ p :=
n i=1

The true L2 -distance is non-negative by deﬁnition (see Eq.(8)), but the above biascorrected estimate can take a negative value. Following the same line as Baranchik (1964),
the positive-part estimator may be more accurate:{}2e2 (X , X ′ ) .
L (X , X ′ ) := max 0, L
2

However, in our preliminary experiments, L (X , X ′ ) does not always perform well particb2 (X , X ′ )ularly when H is ill-conditioned. For this reason, we practically propose to use Ldeﬁned by Eq.(11).

4Experiments

In this section, we experimentally evaluate the performance of LSDD.

4.1

Numerical Examples

First, we show numerical examples using artiﬁcial datasets.

Density-Diﬀerence Estimation

4.1.1

11

LSDD vs. KDE

Letp(x) = N (x; (µ, 0, . . . , 0)⊤ , (4π)−1 I d ),p′ (x) = N (x; (0, 0, . . . , 0)⊤ , (4π)−1 I d ),where N (x; µ, Σ) denotes the multi-dimensional normal density with mean vector µ andvariance-covariance matrix Σ with respect to x, and I d denotes the d-dimensional identity
matrix. Before feeding data samples to algorithms, we pre-normalize them to have unit
variance in the element-wise manner.
In LSDD, the Gaussian width σ and the regularization parameter λ are chosen by
5-fold cross-validation in terms of the LSDD criterion (see Section 2.4) from the following
grid values:σ ∈ {10−2 , 10−1.5 , 10−1 , 10−0.5 , 100 },λ ∈ {10−1 , 10−0.5 , 100 , 100.5 , 101 }.
We experimentally compare the behavior of LSDD with two methods based on the diﬀerence of kernel density estimators (KDEs):
KDEi: Two Gaussian widths are independently chosen from the above candidate values
based on 5-fold least-squares cross-validation. That is, for each density, we perform
cross-validation in terms of the L2 -distance between estimated and true densities sothat the density is optimally approximated (Härdle et al., 2004).
KDEj Two Gaussian widths are jointly chosen from the above candidate values based on
5-fold cross-validation in terms of the LSDD criterion (Hall & Wand, 1988). Thatis, we compute the cross-validated LSDD criterion as a function of two Gaussian
widths and ﬁnd the best pair that minimizes the criterion.
We ﬁrst illustrate the behavior of the LSDD and KDE-based methods under d = 1 andn = n′ = 200. Figure 1 depicts the data samples and density-diﬀerence estimation results
obtained by LSDD, KDEi, and KDEj for µ = 0 (i.e., f (x) = p(x) − p′ (x) = 0). Crossvalidation scores of LSDD are also included at the bottom of the ﬁgure. The ﬁgure showsthat LSDD and KDEj give accurate estimates of the true density diﬀerence f (x) = 0. On
the other hand, the density-diﬀerence estimate obtained by KDEi is rather ﬂuctuated,
although both densities are reasonably well approximated by KDEs. This illustrates anadvantage of directly estimating the density diﬀerence without going through separate
estimation of each density. Figure 2 depicts the results for µ = 0.5 (i.e., f (x) ̸= 0),
showing again that LSDD performs well. KDEi and KDEj give the same estimation
result for this dataset, which slightly underestimates the peaks.
Next, we compare the performance of L2 -distance approximation based on LSDD,
KDEi, and KDEj. For µ = 0, 0.2, 0.4, 0.6, 0.8 and d = 1, 5, we draw n = n′ = 200
samples from the above p(x) and p′ (x). Figure 3 depicts the mean and standard error of

Density-Diﬀerence Estimation

12

40

1.6

xi

f(x)
f(x)

1.4

20

1.2

0

0.8

^

1

−1

−0.5

0
x

0.5

1
0.6

40

0.4

x’

i’

20

0.2
0
−0.2

0
−1

−0.5

0
x

0.5

1

−0.4 b

−0.5

0
x

0.5

1

(b) LSDD (σ = 100 and λ = 101 )(a) Samples
1.6

1.6
p(x)−p’(x)^p(x)−p’(x)
p(x)
^
p(x)
p’(x)
^
p’(x)

1.4

^

1.2
1
0.8

1.2
1
0.8
0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2
−1

−0.5

0
x

0.5

1

(c) KDEi (σ = 10−1 and σ ′ = 10−1 )

f(x)
f(x)

^

1.4

0.6

−0.4 b

−1−0.4 b

−1

−0.5

0.5

1

(d) KDEj (σ = 100 and σ ′ = 100 )σ=10 -2σ=10 -1.5σ=10 -1
σ=10 -1.5σ=10 0

0.16
0.14
0.12
CV score

0
x

0.1
0.08
0.06
0.04
0.02
b

−1

−0.5

0
log10 λ0.51(e) CV scores for LSDD

Figure 1: Estimation of density diﬀerence when µ = 0 (i.e., f (x) = p(x) − p′ (x) = 0).
In the captions of (b), (c), and (d), model parameters chosen by cross-validation are also
described. σ and σ ′ in (c) and (d) denote Gaussian widths for p(x) and p′ (x), respectively.

Density-Diﬀerence Estimation

13

40

1.5

xi

f(x)
f(x)

^

1

20

0.5

0
−1

−0.5

0

0.5

1

1.5

x

0

60
x’

−0.5

i’

40
−1

20
0
−1

−0.5

0

0.5

1

−1.5 b

1.5

−1

−0.5

0

x

0.5

1

1.5

x

(b) LSDD (σ = 10−0.5 and λ = 100 )(a) Samples
b

1.5

1.5
f(x)
f(x)

^

1

1

0.5

0.5

0

0
p(x)−p’(x)^p(x)−p’(x)
p(x)
^
p(x)
p’(x)
^
p’(x)

−0.5

−1

−1.5 b

−0.5

^

−1

−0.5

−1

0

0.5

1

1.5

−1.5 b

−1

−0.5

0

(c) KDEi (σ = 10−1 and σ ′ = 10−1 )11.5(d) KDEj (σ = 10−1 and σ ′ = 10−1 )bσ=10 -2σ=10 -1.5σ=10 -1
σ=10 -1.5σ=10 0

−0.1
−0.2
−0.3
−0.4
CV score

0.5
x

x

−0.5
−0.6
−0.7
−0.8
−0.9
−1
b

−1

−0.5

0
log λ

0.5

1

10

(e) CV scores for LSDD

Figure 2: Estimation of density diﬀerence when µ = 0.5 (i.e., f (x) = p(x) − p′ (x) ̸= 0).
In the captions of (b), (c), and (d), model parameters chosen by cross-validation are also
described. σ and σ ′ in (c) and (d) denote Gaussian widths for p(x) and p′ (x), respectively.

Density-Diﬀerence Estimation

14

True
LSDD
KDE i
KDE j

1.6
1.4

1.4
1.2
L2 distanceL2 distance

1.2
1
0.8

1
0.8

0.6

0.6

0.4

0.4

0.2

0.2

0b

True
LSDD
KDE i
KDE j

1.6

0

0.1

0.2

0.3

0.4
µ

(a) d = 1

0.5

0.6

0.7

0.8

0b

0

0.1

0.2

0.3

0.4
µ

0.5

0.6

0.7

0.8

(b) d = 5

Figure 3: L2 -distance estimation by LSDD, KDEi, and KDEj for n = n′ = 200 as functionsof the Gaussian mean µ. Means and standard errors over 1000 runs are plotted.
estimated L2 -distances over 1000 runs as functions of mean µ. When d = 1 (Figure 3(a)),
the LSDD-based L2 -distance estimator gives the most accurate estimates of the trueL2 -distance, whereas the KDEi-based L2 -distance estimator slightly underestimates thetrue L2 -distance for large µ. This is caused by the fact that KDE tends to provide
smooth density estimates (see Figure 2(c) again): Such smooth density estimates areaccurate as density estimates, but the diﬀerence of smooth density estimates yields asmall L2 -distance estimate (Anderson et al., 1994). More speciﬁcally, the density p′ (x) is
estimated accurately at around x = 0.5, but negative values of the density diﬀerence f (x)are underestimated there because pb(x) is smoother than the true density p(x) and thus
its tail values at around x = 0.5 are larger. The KDEj-based L2 -distance estimator tends
to improve this drawback of KDEi to some extent, but it still slightly underestimates thetrue L2 -distance when µ is large.
When d = 5 (Figure 3(b)), the KDE-based L2 -distance estimators even severely underestimate the true L2 -distance for large µ. On the other hand, the LSDD-based L2 -distanceestimator still gives reasonably accurate estimates of the true L2 -distance even when d = 5.
However, we note that LSDD also slightly underestimates the true L2 -distance when µ
is large, because slight underestimation tends to yield smaller variance and thus suchstabilized solutions are more accurate in terms of the bias-variance trade-oﬀ.
In Figure 1, we have illustrated that LSDD and KDEj work better than KDEi when
the Gaussian mean is µ = 0. However, in Figure 3(a), KDEi is actually shown to be thebest performing method for µ = 0 in terms of the average over 1000 runs. To ﬁll this gap,
let us depict in Figure 4 histograms of L2 -distance estimates obtained by LSDD, KDEi,
and KDEj over 1000 runs for the Gaussian mean µ = 0. This graph shows that LSDD
and KDEj give exactly correct solutions (i.e., zero) about 300 times, whereas KDEi gives
estimates about 0.01 more than 300 times. The graphs plotted in Figure 1 correspond tosuch typical results where LSDD and KDEj outperform KDEi. On the other hand, KDEi

Density-Diﬀerence Estimation

15

L2 distance
LSDD

300
200
100
0

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

KDEi

300
200
100
0

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

KDEj

300
200
100
0

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

Figure 4: Histograms of L2 -distance estimation by LSDD, KDEi, and KDEj over 1000
runs for n = n′ = 200 and the Gaussian mean µ = 0.
stably gives estimates less than 0.1 almost always, whereas LSDD and KDEj occasionallygive large estimates. This rather unstable behavior of LSDD and KDEj, which wascaused by inappropriate choice of the Gaussian width (and the regularization parameter
for LSDD) by cross-validation, led to larger mean values of LSDD and KDEj than KDEi
in Figure 3(a).
Finally, we investigate the behavior of LSDD, KDEi, and KDEj in L2 -distance estimation when the numbers of samples from two distributions are imbalanced. Figure 5
plots the means and standard errors of the L2 -distance estimated by LSDD, KDEi, andKDEj for d = 1, 5 and µ = 0, 0.8 over 1000 runs as functions of n′ when n is ﬁxed to 200.
When d = 1 and µ = 0 (Figure 5(a)), all three methods behave similarly and theaccuracy tends to be improved as n′ increases. However, improvement when n′ > n =200 is moderate. When the input dimensionality is increased to d = 5 (Figure 5(b)),
LSDD and KDEj still have the same tendency. However, KDEi behaves diﬀerently andthe approximation error tends to grow as n′ increases. This implies that improving theaccuracy of one of the density estimates does not necessarily improve the overall estimation
accuracy of the density diﬀerence.

Density-Diﬀerence Estimation

16

True
LSDD
KDE i
KDE j

0.16
0.14

0.4
0.35
L2 distance

0.12
L2 distance

True
LSDD
KDE i
KDE j

0.45

0.1
0.08
0.06

0.3
0.25
0.2
0.15

0.04

0.1

0.02

0.05

0b
50

100

150

200
n’

250

300

0b
50

350

(a) d = 1 and µ = 0

100

150

200
n’

250

300

350

(b) d = 5 and µ = 0
b

True
LSDD
KDE i
KDE j

1.8
1.78
1.76

1.6

1.74

1.2
L2 distanceL2 distance

True
LSDD
KDE i
KDE j

1.4

1.72
1.7
1.68

1
0.8

1.66
1.64

0.6

1.62

0.4

1.6
b

50

b

100

150

200
n’

250

(c) d = 1 and µ = 0.8

300

350

50

100

150

200
n’

250

300

350

(d) d = 5 and µ = 0.8

Figure 5: L2 -distance estimation by LSDD, KDEi, and KDEj for n = 200 as functions ofn′ . Means and standard errors over 1000 runs are plotted.
When d = 1 and µ = 0.8 (Figure 5(c)), LSDD tends to provide better estimates as
n′ increases, whereas KDEi and KDEj keep underestimating the true L2 -distance evenwhen n′ is increased. Finally, when d = 5 and µ = 0.8 (Figure 5(d)), LSDD stably
provides reasonably good results and its performance does not change signiﬁcantly whenn′ is increased. On the other hand, KDEi and KDEj tend to give better results as n′increases.
Overall, LSDD sometimes gives slightly better results for n′ > n, but its performance
is not signiﬁcantly diﬀerent from those for n′ = n. On the other hand, the accuracy of
KDEi and KDEj when n′ is increased gets better or worse depending on the situation.
Thus, having more data samples from one of the distributions does not seem to always
improve the estimation accuracy in density-diﬀerence estimation.

Density-Diﬀerence Estimation

4.1.2

17

L2 -Distance vs. KL-Divergence

The Kullback-Leibler (KL) divergence (Kullback & Leibler, 1951) is a popular divergencemeasure for comparing probability distributions. The KL-divergence from p(x) to p′ (x)is deﬁned as∫p(x)′KL(p∥p ) := p(x) log ′
dx.
p (x)
First, we illustrate the diﬀerence between the L2 -distance and the KL-divergence. For
d = 1, letp(x) = (1 − η)N (x; 0, 12 ) + ηN (x; µ, 1/42 ),p′ (x) = N (x; 0, 12 ).
Implications of the above densities are that samples drawn from N (x; 0, 12 ) are inliers,
whereas samples drawn from N (x; µ, 1/42 ) are outliers. We set the outlier rate at η = 0.1
and the outlier mean at µ = 0, 2, 4, . . . , 10 (see Figure 6).
Figure 7(a) depicts the L2 -distance and the KL-divergence for outlier mean µ =0, 2, 4, . . . , 10. This shows that both the L2 -distance and the KL-divergence increase as µincreases. However, the L2 -distance is bounded from above, whereas the KL-divergence
diverges to inﬁnity as µ tends to inﬁnity. This result implies that the L2 -distance is lesssensitive to outliers than the KL-divergence, which well agrees with the observation given
in Basu et al. (1998).
Next, we draw n = n′ = 100 samples from p(x) and p′ (x), and estimate the L2 distance by LSDD and the KL-divergence by the Kullback-Leibler importance estimation
procedure 5 (KLIEP) (Sugiyama et al., 2008). Figure 7(b) depicts estimated L2 -distanceand KL-divergence for outlier mean µ = 0, 2, 4, . . . , 10 over 100 runs. This shows thatboth LSDD and KLIEP reasonably capture the proﬁles of the true L2 -distance and the
KL-divergence, although the scale of KLIEP values is much diﬀerent from the true values(see Figure 7(a)) because the estimated normalization factor was unreliable.
Finally, based on the permutation test procedure (Efron & Tibshirani, 1993), we conduct hypothesis testing of the null hypothesis whether densities p and p′ are the same.
More speciﬁcally, we ﬁrst compute a distance estimate for the original datasets X andb , X ′ ). Next, we randomly permute theX ′ and obtain a distance/divergence estimate D(X|X ∪ X ′ | samples, and assign the ﬁrst |X | samples to a set Xe and the remaining |X ′ |samples to another set Xe′ . Then we compute a distance/divergence estimate again usinge Xe, Xe′ ). Because Xe and Xe′ can
the randomly permuted datasets Xe and Xe′ and obtain D(e Xe, Xe′ ) would take a value close
be regarded as being drawn from the same distribution, D(
5

Estimation of the KL-divergence from data has been extensively studied recently (Wang et al., 2005;
Sugiyama et al., 2008; Pérez-Cruz, 2008; Silva & Narayanan, 2010; Nguyen et al., 2010). Among them,
KLIEP was shown to possess a superior convergence property and demonstrated to work well in practice(Sugiyama et al., 2008). KLIEP is based on direct estimation of density ratio p(x)/p′ (x) without density
estimation of p(x) and p′ (x). See also Nguyen et al. (2010), which proposes essentially the sameprocedure.

Density-Diﬀerence Estimation

18

b

µ=0

p(x)

0.4
µ=2
0.2
0b

−2

0

2

µ=4

µ=6

µ=8

µ=10

4
x

6

8

10
b

0.4
p’(x)
0.3
0.2
0.1
0b

−2

0

2

4
x

6

8

10

Figure 6: Comparing two densities in the presence of outliers. p(x) includes outliers atµ = 0, 2, 4, . . . , 10.
b

0.02
0.01

0.01
0.005

L2−LSDD

L2−distance
0b
0

2

4

6

8

0
0

10

2

4

6

8

10

6

8

10

µ

µ
b

6

KL−KLIEP

KL−divergence
4

0.2

2

0.1
0

0b
0

2

4

6

8

0

10

2

4
µ

µ

(a) True L2 -distance and true KL-divergence (b) Means and standard errors of L2 -distanceestimation by LSDD and KL-divergence estimation by KLIEP over 100 runs.

Figure 7: Comparison of L2 -distance and KL-divergence for outlier rate η = 0.1 as functions of outlier mean µ.
1

0.9

0.8

0.8

0.7

0.7
Rejection rate

Rejection rate

1

L2−LSDD
KL−KLIEP

0.9

0.6
0.5
0.4

0.6
0.5
0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0

2

4

µ

6

(a) η = 0.1

810

L2−LSDD
KL−KLIEP
0

0.1

0.2
η

0.3

0.4

(b) µ = 10

Figure 8: Results of two-sample test over 1000 runs for (a) outlier rate η = 0.1 as functionsof outlier mean µ and for (b) outlier mean µ = 10 as functions of outlier rate η.

Density-Diﬀerence Estimation

19

to zero. This random permutation procedure is repeated many times (100 times in thee Xe, Xe′ ) under the null hypothesis (i.e.,
following experiments), and the distribution of D(
the two distributions are the same) is constructed. Finally, the p-value is approximatedb , X ′ ) in the histogram of D(e Xe, Xe′ ). We set theby evaluating the relative ranking of D(Xsigniﬁcance level at 5%.
Figure 8(a) depicts the rejection rate of the null hypothesis (i.e., p = p′ ) over 1000
runs for outlier rate η = 0.1 and outlier mean µ = 0, 2, 4, . . . , 10, based on the L2 -distanceestimated by LSDD and the KL-divergence estimated by KLIEP. This shows that the
KLIEP-based test rejects the null hypothesis more frequently for large µ, whereas therejection rate of the LSDD-based test is almost unchanged with respect to µ.
Figure 8(b) depicts the rejection rate of the null hypothesis for outlier mean µ = 10 andoutlier rate η = 0, 0.1, 0.2, . . . , 0.4. When η = 0 (i.e., no outliers), both the LSDD-based
test and the KLIEP-based test accept the null hypothesis with the designated signiﬁcancelevel approximately. When η = 0.1, the LSDD-based test still keeps a low rejection rate,
whereas the KLIEP-based test tends to reject the null hypothesis more frequently. Whenη ≥ 0.3, both the LSDD-based test and the KLIEP-based test always reject the null
hypothesis.
Overall, the above results imply that the two-sample homogeneity test by LSDD ismore robust against outliers (i.e., two distributions tend to be regarded as the same evenin the presence of outliers) than the KLIEP-based test.

4.2

Applications

Here, we apply LSDD to semi-supervised class-balance estimation under class-prior change
and change-point detection in time series.
4.2.1

Semi-Supervised Class-Balance Estimation

In real-world pattern recognition tasks, changes in class balance between the training andtest phases are often observed. In such cases, naive classiﬁer training produces signiﬁcantestimation bias because the class balance in the training dataset does not properly reﬂectthat of the test dataset. Here, we consider the problem of learning the class balance of atest dataset in a semi-supervised learning setup where unlabeled test samples are provided
in addition to labeled training samples (Chapelle et al., 2006).
More formally, we consider the binary classiﬁcation problem of classifying patternx ∈ Rd to class y ∈ {+1, −1} under class-prior change, where the class-prior probability
for training data ptrain (y) and that for test data ptest (y) are diﬀerent:
ptrain (y) ̸= ptest (y).
However, we assume that the class-conditional density for training data ptrain (x|y) andthat for test data ptest (x|y) are unchanged:
ptrain (x|y) = ptest (x|y).

Density-Diﬀerence Estimation

20|
ptrain(x|y = −1)

ptrain(x|y = +1)|
ptest (x)train

xFigure 9: Schematic illustration of semi-supervised class-balance estimation.
Note that training and test joint densities ptrain (x, y) and ptest (x, y) as well as training
and test input densities ptrain (x) and ptest (x) are generally diﬀerent under this setup.
Here, our objective is to estimate ptest (y) from labeled training samples {(xi , yi )}ni=1′drawn independently from ptrain (x, y) and unlabeled test samples {x′i′ }ni′ =1 drawn inde′
pendently from ptest (x). Given test labels {yi′′ }ni′ =1 , ptest (y) can be naively estimated byn′y /n′ , where n′y is the number of test samples in class y. Here, however, we want to′estimate ptest (y) without {yi′′ }ni′ =1 .
The class balance in the test set can be estimated by matching a mixture of class-wisetraining input densities,
qtest (x; π) := πptrain (x|y = +1) + (1 − π)ptrain (x|y = −1),to the test input density ptest (x) (Saerens et al., 2002), where π ∈ [0, 1] is a mixing coeﬃcient to learn. See Figure 9 for schematic illustration. Here, we use the L2 -distanceestimated by LSDD, KDEi, and KDEj (see Section 4.1.1) for this distribution matching. Note that, when LSDD is used to estimate the L2 -distance, separate estimationof ptrain (x|y = ±1) is not involved, but the diﬀerence between ptest (x) and qtest (x; π) is
directly estimated.
As an additional baseline, we include the EM-based class-prior estimation method(Saerens et al., 2002), which actually corresponds to distribution matching under the KL
divergence. More speciﬁcally, in the EM-based algorithm, test class-prior estimate pbtest (y)and test class-posterior estimate pbtest (y|x) are iteratively estimated as follows:
1. Obtain an estimate of the training class-posterior probability, pbtrain (y|x), from labeled training samples {(xi , yi )}ni=1 , for example, by kernel logistic regression (Hastie
et al., 2001) or its squared-loss variant (Sugiyama, 2010).
2. Obtain an estimate of the training class-prior probability from training data{(xi , yi )}ni=1 as pbtrain (y) = ny /n, where ny is the number of training samples inclass y.
3. Set the initial estimate of the test class-posterior probability aspbtest (y) ←− pbtrain (y).

Density-Diﬀerence Estimation

21

4. Compute a new test class-posterior estimate pbtest (y|x) based on the current testclass-prior estimate pbtest (y) as
pbtest (y)b
ptrain (y|x)/b
ptrain (y).
′′btest (y )bptrain (y |x)/b
ptrain (y ′ )y ′ =1 p

pbtest (y|x) ←− ∑c

5. Compute a new test class-prior estimate pbtest (y) based on the current test class-priorestimate pbtest (y|x) asn′

1 ∑pbtest (y) ←− ′pbtest (y|x′i′ ).
n i′ =1
6. Iterate 4. and 5. until convergence.
We use four UCI benchmark datasets6 for experiments, where we randomly choose
10 labeled training samples from each class and 50 unlabeled test samples following trueclass-prior
π ∗ = 0.1, 0.2, . . . , 0.9.
The left graphs in Figure 10 plot the mean and standard error of the squared diﬀerence
between true and estimated class-balances π. These graphs show that LSDD tends toprovide better class-balance estimates than alternative approaches.
Next, we use the estimated class balance to train a classiﬁer. We use a weighted ℓ2 regularized least-squares classiﬁer (Rifkin et al., 2003). That is, a class label yb for a test
input x is estimated by)( n∑
yb = signαbℓ K(x, xℓ ) ,ℓ=1

where K(x, x′ ) is the Gaussian kernel function with kernel width κ. {αbℓ }nℓ=1 are learned
parameters given by


( n
)2
n
n
∑
∑
∑
πyi
(b
α1 , . . . , αbn ) := argmin 
αℓ2  ,αℓ K(xi , xℓ ) − yi + δ
nyi /n ℓ=1
α1 ,...,αni=1ℓ=1where π+1 = π, π−1 = 1 − π, and δ (≥ 0) is the regularization parameter. The Gaussian
width κ and the regularization parameter δ are chosen by 5-fold weighted cross-validation(Sugiyama et al., 2007) in terms of the misclassiﬁcation error.
The right graphs in Figure 10 plot the test misclassiﬁcation error over 1000 runs. The
results show the LSDD-based method provides lower classiﬁcation errors, which would be
brought by good estimates of test class-balances.
6

http://archive.ics.uci.edu/ml/

Density-Diﬀerence Estimation

LSDD
KDEi

0.5

KDEj
EM

0.2

Misclassification rate

Class balance squared error

0.25

22

0.15
0.1
0.05
0

0.2

0.4

π*

0.6

KDEj
EM

0.4
0.3
0.2
0.1

0.8

LSDD
KDEi

0.2

0.4

π*

0.6

0.8

LSDD
KDEi

0.15

KDEj
EM

0.1

0.05
0.2

0.4

π*

0.6

0.45
Misclassification rate

Class balance squared error(a) Australian dataset
0.4
0.35
0.3
0.25

LSDD
KDEi

0.2

0.8

0.2

0.4

KDEj
EM
π*

0.6

0.8

LSDD
KDEi

0.3

KDEj
EM

0.5
Misclassification rate

Class balance squared error(b) Diabetes dataset

0.25
0.2
0.15
0.1
0.05
0

0.2

0.4

π*

0.6

0.4

0.3

0.2

0.8

LSDD
KDEi
0.2

0.4

π*

KDEj
EM
0.6

0.8

(c) German dataset
LSDD
KDEi

0.2

KDEj
EM

0.15
0.1
0.05
0

0.2

0.4

π*

0.6

0.8

0.5
Misclassification rate

Class balance squared error

0.25

LSDD
KDEi

KDEj
EM

0.4
0.3
0.2
0.1

0.2

0.4

π*

0.6

0.8

(d) Statlogheart dataset

Figure 10: Results of semi-supervised class-balance estimation. Left: Squared error ofclass balance estimation. Right: Misclassiﬁcation error by a weighted ℓ2 -regularized leastsquares classiﬁer with weighted cross-validation.

Density-Diﬀerence Estimation

23

Y (t + 2ry(t + r)

Y
y(t)

Time
a

r

b

c

Y (t) a b
Y (t + 1) b

c
c

d

d

e

f

g

h

f

g

h

g

h

Y(
r
k
f
g
Y (t + r − 1) e Y(
+
Y(t)

i

j

k

l

( + −
Y (t + r)i
Y (t + r + 1)jl Y (t + 2r − 1)(k +)Y(t + r)Figure 11: Schematic illustration of unsupervised change detection.
4.2.2

Unsupervised Change Detection

The objective of change detection is to discover abrupt property changes behind timeseries data.
Let y(t) ∈ Rm be an m-dimensional time-series sample at time t, and let
Y (t) := [y(t)⊤ , y(t + 1)⊤ , . . . , y(t + k − 1)⊤ ]⊤ ∈ Rkm
be a subsequence of time series at time t with length k. We treat the subsequence Y (t)as a sample, instead of a single point y(t), by which time-dependent information can beincorporated naturally (Kawahara & Sugiyama, 2012). Let Y(t) be a set of r retrospective
subsequence samples starting at time t:
Y(t) := {Y (t), Y (t + 1), . . . , Y (t + r − 1)}.
Our strategy is to compute a certain dissimilarity measure between two consecutive
segments Y(t) and Y(t + r), and use it as the plausibility of change points (see Figure 11).
As a dissimilarity measure, we use the L2 -distance estimated by LSDD and the KLdivergence estimated by the KL importance estimation procedure (KLIEP) (Sugiyama
et al., 2008). We set k = 10 and r = 50.
We use two datasets. One is the IPSJ SIG-SLP Corpora and Environments for Noisy
Speech Recognition (CENSREC) dataset7 . This dataset was provided by the National
Institute of Informatics, Japan that records human voice in a noisy environment suchas a restaurant. Another dataset is the Human Activity Sensing Consortium (HASC)
challenge 2011 8 , which provides human activity information collected by portable threeaxis accelerometers. Because the orientation of the accelerometers is not necessarily ﬁxed,
we take the ℓ2 -norm of the 3-dimensional data. The HASC dataset is relatively simple,
so we artiﬁcially add zero-mean Gaussian noise with standard deviation 5 at each time
point with probability 0.005.
78

http://research.nii.ac.jp/src/en/CENSREC-1-C.html
urlhttp://hasc.jp/hc2011/

Density-Diﬀerence Estimation

24

0.1

Original data

0
−0.1
−0.2

0

500

1000

1500

1000

1500

1000

1500

KLIEP score

40
20
0

0

500
LSDD score

1
0

0

500
Time

(a) CENCREC dataset

5

Original
data

0
−5
0

500

1000

1500

2000

2500

3000

2000

2500

3000

2000

2500

3000

KLIEP score

40
20
0

0

500

1000

1500
LSDD score

2
1
0

0

500

1000

1500
Time

(b) HASC dataset

Figure 12: Illustrative results of unsupervised change detection for (a) CENSREC speech
data and (b) HASC dataset. Original time-series data is plotted in the top graphs, andchange scores obtained by KLIEP and LSDD are plotted in the bottom graphs.

Density-Diﬀerence Estimation

25

The top graphs in Figure 12 display the original time-series, where true change points
were manually annotated. The time-series data in Figure 12(a) corresponds to a sequenceof “noise”, “speech”, “noise”, “speech”, and “noise”, whereas that in Figure 12(b) corresponds to a sequence of actions “jog”, “stay”, “stair down”, “stay”, and “stair up”.
The bottom graphs in Figure 12 plot change scores obtained by each method. The
results show that the LSDD-based change score indicates the existence of change points
more clearly than the KLIEP-based approach. The superior performance of LSDD over
the KLIEP-based change score would be brought by its robustness against outliers (seeSection 4.1.2).
Finally, we compare the change-detection performance more systematically using thereceiver operating characteristic (ROC) curves (i.e., the false positive rate vs. the true
positive rate) and the area under the ROC curve (AUC) values. In addition to LSDD and
KLIEP, we also test the L2 -distance estimated by KDEi and KDEj (see Section 4.1.1).
Moreover, in our comparison, we also include native change detection methods based on
autoregressive models (AR) (Takeuchi & Yamanishi, 2006), subspace identiﬁcation (SI)
(Kawahara et al., 2007), singular spectrum transformation (SST) (Moskvina & Zhigljavsky, 2003), one-class support vector machine (SVM) (Desobry et al., 2005), kernel
Fisher discriminant analysis (KFD) (Harchaoui et al., 2009), and kernel change-point detection (KCP) (Arlot et al., 2012). Tuning parameters included in these methods weremanually optimized.
We use 10 datasets taken from each of the CENSREC and HASC data collections.
Mean ROC curves are plotted in Figure 13 and AUC values are described in Table 1.
The experimental results show that LSDD tends to outperform other methods and iscomparable to state-of-the-art native change-detection methods.

5Conclusions

In this paper, we proposed a method for directly estimating the diﬀerence between two
probability density functions without density estimation. The proposed method, called
the least-squares density-diﬀerence (LSDD), was derived within the framework of kernel
regularized least-squares estimation, and its solution can be computed analytically in acomputationally eﬃcient and stable manner. Furthermore, LSDD is equipped with crossvalidation, and thus all tuning parameters such as the kernel width and the regularization
parameter can be systematically and objectively optimized. We showed the asymptotic
normality of LSDD in a parametric setup and derived a ﬁnite-sample error bound for
LSDD in a non-parametric setup. In both cases, LSDD was shown to achieve the optimal
convergence rates.
We also proposed an L2 -distance estimator based on LSDD, which nicely cancels thebias caused by regularization. The LSDD-based L2 -distance estimator was experimentally shown to be more accurate than diﬀerences of kernel density estimators and morerobust against outliers than a Kullback-Leibler divergence estimator. However, we also experimentally observed that cross-validation of LSDD is sometimes rather unstable whenDensity-Diﬀerence Estimation

26

1
0.9

True positive rate

0.8
0.7
0.6

LSDD
KDEi
KDEj
KLIEP
AR
SI
SST
SVM
KFD
KCP

0.5
0.4
0.3
0.2
0.1
0

0

0.2

0.4
0.6
False positive rate

0.81(a) CENCREC dataset
1
0.9

True positive rate

0.8
0.7
0.6

LSDD
KDEi
KDEj
KLIEP
AR
SI
SST
SVM
KFD
KCP

0.5
0.4
0.3
0.2
0.1
0

0

0.2

0.4
0.6
False positive rate

0.81(b) HASC dataset

Figure 13: Mean ROC curves of unsupervised change detection.

Density-Diﬀerence Estimation

27

Table 1: AUC values of unsupervised change detection. The best method and comparable
ones in terms of mean AUC values by the t-test (Henkel, 1976) at the signiﬁcance level
5% are indicated with boldface. “SE” stands for “Standard error”.
(a) CENSREC dataset

Data ID
1
2
3
4
5
6
7
8
9
10
Mean
SE

LSDD KDEi
.888
.737
.871
.803
.910
.753
.936
.823
.878
.712
.830
.732
.813
.727
.889
.841
.828
.739
.943
.687
.879
.755
.014
.016

KDEj
.731
.706
.690
.578
.799
.711
.737
.734
.586
.773
.705
.023

KLIEP
.437
.618
.744
.683
.667
.696
.513
.691
.609
.692
.635
.030

AR
.769
.777
.762
.776
.768
.679
.727
.783
.776
.670
.749
.013

SI SST.739 .507.736 .541.821 .616.816 .723.701 .625.727 .484.733 .612.775 .526.770 .609.747 .551.756 .580.012 .023

SVM
.604
.612
.886
.871
.843
.781
.779
.698
.819
.835
.773
.032

KFD
.881
.912
.876
.981
.880
.841
.938
.934
.922
.889
.905
.013

KCP
.917
.879
.743
.826
.945
.947
.968
.935
.980
.984
.913
.024

SI SST.690 .806.800 .745.899 .807.735 .685.823 .809.740 .736.759 .797.704 .682.744 .781.725 .790.762 .764.020 .016

SVM
.800
.725
.932
.751
.840
.838
.829
.800
.790
.850
.815
.018

KFD
.885
.904
.707
.903
.961
.871
.770
.852
.842
.866
.856
.023

KCP
.874
.826
.641
.759
.725
.800
.532
.661
.697
.787
.730
.032

(b) HASC dataset

Data ID
1
2
3
4
5
6
7
8
9
10
Mean
SE

LSDD KDEi
.792
.823
.842
.665
.773
.605
.921
.839
.838
.849
.834
.755
.841
.763
.878
.833
.864
.850
.847
.663
.843
.764
.013
.029

KDEj
.753
.741
.536
.837
.859
.781
.598
.857
.866
.680
.751
.036

KLIEP AR
.650
.860
.712
.733
.708
.910
.587
.816
.565
.831
.676
.868
.657
.807
.581
.629
.693
.738
.554
.796
.638
.799
.020
.026

Density-Diﬀerence Estimation28

the target density diﬀerence is zero (i.e., two distributions are equivalent). This canpotentially cause performance degradation in two-sample homogeneity testing because
estimation of zero density-diﬀerence is repeatedly executed when approximating the null
distribution in the permutation-test framework. Stabilizing cross-validation and improving the accuracy of density-diﬀerence estimation when the target density-diﬀerence is zerois a remaining future work.
It is straightforward to extend the proposed LSDD method to the diﬀerence of weighted
densities
νp(x) − ν ′ p′ (x),where ν and ν ′ are scalars. Also, LSDD can be easily extended to estimate the weighted
L2 -distance:∫ ()2p(x) − p′ (x) w(x)dx,
where w(x) > 0 is a weight function.
A related line of research to density-diﬀerence estimation is density-ratio estimation(Sugiyama et al., 2012a), which directly estimates the ratio of probability densities without
separate density estimation (Qin, 1998; Huang et al., 2007; Bickel et al., 2007; Sugiyama
et al., 2008; Kanamori et al., 2009; Sugiyama et al., 2012b). Potential weaknesses ofdensity-ratio estimation are that density ratios can be unbounded even for simple cases(Cortes et al., 2010) and their estimation may suﬀer from outliers (Basu et al., 1998;
Scott, 2001; Besbeas & Morgan, 2004).
To mitigate these weaknesses, the concept of relative density ratios was introducedp(x)recently, which “ﬂatten” the density ratio pp(x)′ (x) as βp(x)+(1−β)p′ (x) for 0 ≤ β < 1 (Yamada
et al., 2013). Even when the plain density ratio is unbounded, the relative density ratio
is always bounded by β1 for β > 0. Although estimation of relative density ratios as well
as approximation of relative divergences was demonstrated to be more reliable (Yamada
et al., 2013), there is no systematic method to choose the relativity parameter β, which
is a critical limitation in practice.
On the other hand, density-diﬀerence estimation is more advantageous than densityratio estimation in the senses that density diﬀerences are always bounded as long as
each density is bounded, their estimation is robust against outliers (Basu et al., 1998;
Scott, 2001; Besbeas & Morgan, 2004), and there exist no tuning parameters such as therelativity parameter β. However, a potential weakness of density diﬀerences is that theycannot be used for importance sampling (Sugiyama & Kawanabe, 2012) and conditional
probability estimation (Sugiyama et al., 2010b; Sugiyama, 2010), which are promising
usages of density-ratio estimation. Thus, further exploring usages of density-diﬀerence
estimation, particularly in the tasks that density-ratio estimation cannot be used for, isa promising future work.
A simple application of density-diﬀerence estimation would be probabilistic patternrecognition, because the sign of the density diﬀerence gives the Bayes-optimal decision(Duda et al., 2001). Furthermore, in the context of pattern recognition with a reject
option, the density diﬀerence can be used for ﬁnding the optimal rejection thresholdDensity-Diﬀerence Estimation

29(Chow, 1970). In the future work, we will investigate the behavior of LSDD in probabilistic
pattern recognition theoretically and experimentally.
Density-diﬀerence estimation is a novel research paradigm in machine learning, andwe have proposed a simple but useful method for this emerging topic. Our future work
will develop more powerful algorithms for density-diﬀerence estimation. For example,
considering more general loss functions than the squared loss (Sugiyama et al., 2012b)
and incorporating dimension reduction (von Bünau et al., 2009; Sugiyama et al., 2010a;
Sugiyama et al., 2011; Yamada & Sugiyama, 2011) would be interesting directions topursue. Exploring a wide variety of real-world applications is also an important futurework.

Acknowledgments
The authors would like to thank Wittawat Jitkrittum and John Quinn for their comments
and Zaı̈d Harchaoui for providing us a program code of kernel change-point detection.
Masashi Sugiyama was supported by MEXT KAKENHI 23300069 and AOARD, Takafumi
Kanamori was supported by MEXT KAKENHI 24500340, Taiji Suzuki was supported
by MEXT KAKENHI 22700289, the Aihara Project, the FIRST program from JSPS
initiated by CSTP, and the Global COE Program “The research and training center for
new development in mathematics”, MEXT, Japan, Marthinus Christoﬀel du Plessis was
supported by MEXT Scholarship, Song Liu was supported by the JST PRESTO program,and Ichiro Takeuchi was supported by MEXT KAKENHI 23700165.

A

Technical Details of Non-Parametric Convergence
Analysis in Section 2.3.2

First, we deﬁne the linear operators Pn , P, Pn′ , P ′ , Qn , Q as∫n
1∑
Pn f :=f (xi ), P f :=f (x)p(x)dx,
n i=1
Rd
∫
n
1∑
′
′
′
Pn f :=
f (x)p′ (x)dx,
f (xi ), P f :=n i=1
RdQn f := Pn f − Pn′ f, Qf := P f − P ′ f.
Let Hγ be an RKHS endowed with the Gaussian kernel with width γ:()∥x − x′ ∥2′kγ (x, x ) = exp −.
γ2
A density-diﬀerence estimator fb is obtained as
]
[
2
2
b
f := argmin ∥f ∥L2 (Rd ) − 2Qn f + λ∥f ∥Hγ .
f ∈Hγ

Density-Diﬀerence Estimation30

We assume the following conditions:
Assumption 1. The densities are bounded: There exists a constant M such that∥p∥∞ ≤ M and ∥p′ ∥∞ ≤ M.
The density diﬀerence f = p − p′ is a member of Besov space with regularity α. That is,ααf ∈ B2,∞where B2,∞is the Besov space with regularity α, andα∥f ∥B2,∞:= ∥f ∥L2 (Rd ) + sup(t−α ωr,L2 (Rd ) (f, t)) < ct>0

for r = ⌊α⌋ + 1,

where ⌊α⌋ denotes the largest integer less than or equal to α and ωr,L2 (Rd ) is the r-th
modulus of smoothness (see Eberts and Steinwart (2011) for the deﬁnitions).
Then we have the following theorem.
Theorem 2. Suppose Assumption 1 is satisﬁed. Then, for all ϵ > 0 and p ∈ (0, 1), there
exists a constant K > 0 depending on M, c, ϵ, p such that for all n ≥ 1, τ ≥ 1, and λ > 0,
the LSDD estimator fb in Hγ satisﬁes()
2(1−p)d− 1+p (1+ϵ+ 1−p)−(1−p)(1+ϵ)d
4
γ
γ
τ
τ
∥fb − f ∥2L2 (Rd ) +λ∥fb∥2Hγ ≤ K λγ −d +γ 2α +++ 2 + ,
3p−p2
2
λp n
nλ n
λ 1+p n 1+pwith probability not less than 1 − 4e−τ .
To prove this, we utilize the technique developed in Eberts and Steinwart (2011) fora regression problem.
Proof. First, note that∥fb∥2L2 (Rd ) − 2Qn fb + ∥f ∥2L2 (Rd ) + λ∥fb∥2Hγ ≤ ∥f0 ∥2L2 (Rd ) − 2Qn f0 + ∥f ∥2L2 (Rd ) + λ∥f0 ∥2Hγ .
Therefore, we have∥fb − f ∥2L2 (Rd ) + λ∥fb∥2Hγ= ∥fb∥2L2 (Rd ) − 2Qn fb + ∥f ∥2L2 (Rd ) + 2(Qn − Q)fb + λ∥fb∥2Hγ≤ ∥f0 ∥2L2 (Rd ) − 2Qn f0 + ∥f ∥2L2 (Rd ) + 2(Qn − Q)fb + λ∥fb∥2Hγ= ∥f0 ∥2L2 (Rd ) − 2Qf0 + ∥f ∥2L2 (Rd ) + 2(Qn − Q)(fb − f0 ) + λ∥fb∥2Hγ= ∥f0 − f ∥2L2 (Rd ) + 2(Qn − Q)(fb − f ) + 2(Qn − Q)(f − f0 ) + λ∥fb∥2Hγ .
Let

(
)2
(
)
r ( )
∑
2
2∥x∥2
r
1−j 1
√
exp − 2 2 ,K(x) :=(−1)
jd γ π
j γ
j
j=1
d

(16)

Density-Diﬀerence Estimation31√dand fe(x) := (γ π)− 2 f . Using K and fe, we deﬁne∫
f0 := K ∗ fe :=fe(y)K(x − y)dy,Rdi.e., f0 is the convolution of K and fe. Because of Lemma 2 in Eberts and Steinwart(2011), we have f0 ∈ Hγ and∥f0 ∥Hγ ≤ (2r − 1)∥fe∥L2 (Rd ) (∵ Lemma 2 of Eberts and Steinwart (2011))√d≤ (2r − 1)(γ π)− 2 ∥f ∥L2 (Rd )√d≤ (2r − 1)(γ π)− 2 (∥p∥L2 (Rd ) + ∥p′ ∥L2 (Rd ) )√d √≤ (2r − 1)(γ π)− 2 2 M .

(17)

Moreover, Lemma 3 in Eberts and Steinwart (2011) gives∥f0 ∥∞ ≤ (2r − 1)∥f ∥∞ ≤ (2r − 1)M,

(18)

and Lemma 1 in Eberts and Steinwart (2011) yields that there exists a constant Cr,2 suchthat
γ2∥f0 − f ∥2L2 (Rd ) ≤ Cr,2 ωr,L
) ≤ Cr,2 c2 γ 2α .
(19)
2 (Rd ) (f,2
Now, following a similar line to Theorem 3 in Eberts and Steinwart (2011), we canshow that, for all ϵ > 0 and p ∈ (0, 1), there exists a constant Cϵ,p such that|(Pn − P )(fb − f )| ≤ fb − f.
To bound this, we derive the tail probability of()
fb − f(Pn − P ),∥fb − f ∥2L2 (Rd ) + λ∥fb∥2Hγ + r
where r > 0 is a positive real number such that r > r∗ forr∗ = min ∥f − f ∥2L2 (Rd ) + λ∥f ∥2Hγ .
f ∈Hγ

Let
gf,r =f −f∥f − f ∥2L2 (Rd ) + λ∥f ∥2Hγ + r

for f ∈ Hγ and r > r∗ . Then we have∥gf,r ∥∞ ≤≤∥f ∥Hγ + ∥f ∥∞∥f ∥∞ + ∥f ∥∞≤2
2∥f − f ∥L2 (Rd ) + λ∥f ∥Hγ + r∥f − f ∥2L2 (Rd ) + λ∥f ∥2Hγ + r
M
1
1
M
+
≤ √ + ,
λ∥f ∥Hγ + r/∥f ∥Hγ
r
r
2 rλ

Density-Diﬀerence Estimation

32

and
2
P gf,r=M ∥f − f ∥2L2 (Rd )M
P (f − f )2
≤
≤
.
2
2
2
2
2
2
(∥f − f ∥L2 (Rd ) + λ∥f ∥Hγ + r)(∥f − f ∥L2 (Rd ) + λ∥f ∥Hγ + r)r

Here, let
Fr := {f ∈ Hγ | ∥f − f ∥2L2 (Rd ) + λ∥f ∥2Hγ ≤ r},
and we assume that there exists a function such that[]E sup |(Pn − P )(f − f )| ≤ φn (r),f ∈Frwhere E denotes the expectation over all samples. Then, by the peeling device (seeTheorem 7.7 in Steinwart & Christmann, 2008), we have
E sup |(Pn − P )gf,r | ≤
f ∈Hγ

8φ(r)
.
r

Therefore, by Talagrand’s concentration inequality, we have
[
√
(
)]
1
10φn (r)
2M τ
14τ
M
√ +
Pr sup |(Pn − P )gf,r | <++≥ 1 − e−τ ,rnr
3n 2 rλ
r
f ∈Hγ

(20)

where Pr[·] denotes the probability of an event.
From here on, we give an upper bound of φn . The RKHS Hγ can be embedded inarbitrary Sobolev space W m (Rd ). Indeed, by the proof of Theorem 3.1 in Steinwart andScovel (2007), we havemd∥f ∥W m (Rd ) ≤ Cm γ − 2 + 4 ∥f ∥Hγ
for all f ∈ Hγ . Moreover, the theories of interpolation spaces give that, for all f ∈W m (Rd ), the supremum norm of f can be bounded as
1− d

d

′
2m
2m
∥f ∥∞ ≤ Cm
∥f ∥L2 (Rd ) ∥f ∥W m (Rd ) ,d
if d < 2m. Here we set m = 2p
. Then we have∥f ∥∞ ≤ Cp′′ ∥f ∥1−p
∥f ∥pHγ γ −L2 (Rd )d(1−p)
4.

Now, because Fr ⊂ (r/λ)1/2 BHγ and
P (f − f )2 ≤ M ∥f − f ∥2L2 (Rd ) ≤ M r

for f ∈ Fr

hold from Theorem 7.16 and Theorem 7.34 in Steinwart and Christmann (2008), we cantake{( r ) p2(1−p)(1+ϵ)d1−p2(M r) 2 n−1/2 ,φn (r) = max C1,p,ϵ γ −
λ
}
] 1−p
p [(
) p2
( r ) 1+p1+p(1−p)(1+ϵ)dd(1−p) 1−pr
C2,p,ϵ γ − 1+pn−1/(1+p) ,γ− 4 r 2λλDensity-Diﬀerence Estimation

33

where ϵ > 0 and p ∈ (0, 1) are arbitrary and C1,p,ϵ , C2,p,ϵ are constants depending on p, ϵ.
In the same way, we can also obtain a bound of supf ∈Hγ |(Pn′ − P ′ )gf,r |.
If we set r to satisfy√()
1
10φn (r)
14τ
1
2M τ
M
√ +
≥
+
+
,
(21)
8
r
nr
3n 2 rλrthen we have|(Qn − Q)(fb − f )| ≤)1(
r + ∥fb − f ∥2L2 (Rd ) + λ∥fb∥Hγ
4with probability 1 − 2e−τ . To satisfy Eq.(21), it suﬃces to set)(2(1−p)d)− 1+p (1+ϵ+ 1−p−(1−p)(1+ϵ)d
4
γ
τ
τ
γ
r=C
+
+ 2 +
,
3p−p2
2
p
λn
nλ n
λ 1+p n 1+p

(22)

(23)

where C is a suﬃciently large constant depending on M, ϵ, p.
Finally, we bound the term (Qn − Q)(f0 − f ). By Bernstein’s inequality, we have√()τ
2r M τ
|(Pn − P )(f0 − f )| ≤ C ∥f − f0 ∥L2 (P )
+
n
n
√
(
)
√
τ
2r M τ≤C
2M ∥f − f0 ∥L2 (Rd )
+
n
n
)
(
r
2 Mτ
2M τ
2
+
,
(24)
≤ C ∥f − f0 ∥L2 (Rd ) +nnwith probability 1 − e−τ , where C is a universal constant. In a similar way, we can also
obtain()
2M τ
2r M τ
′′2|(Pn − P )(f0 − f )| ≤ C ∥f − f0 ∥L2 (Rd ) ++.
nnCombining these inequalities, we have()
2r M τ
2|(Qn − Q)(f0 − f )| ≤ C ∥f − f0 ∥L2 (Rd ) +,(25)nwith probability 1 − 2e−τ , where C is a universal constant.
Substituting Eqs.(22) and (25) into Eq.(16), we have∥fb − f ∥2L2 (Rd ) + λ∥fb∥2Hγ)}{(
2r M τ
22+ r + λ∥f0 ∥Hγ ,≤ 2 ∥f0 − f ∥L2 (Rd ) + C ∥f − f0 ∥L2 (Rd ) +nwith probability 1 − 4e−τ . Moreover, by Eqs.(19) and (17), the right-hand side is further
bounded by
{
}
1+τ
2α
−d
2
2
b
b
,
∥f − f ∥L2 (Rd ) + λ∥f ∥Hγ ≤ C γ + r + λγ +nDensity-Diﬀerence Estimation

34

Finally, substituting Eq.(23) into the right-hand side, we have∥fb − f ∥2L2 (Rd ) + λ∥fb∥2Hγ}{2(1−p)d)− 1+p (1+ϵ+ 1−p−(1−p)(1+ϵ)d
4
γ
γ
τ
τ
+
+ λγ −d + 2 +≤ C γ 2α +
3p−p2
2
λp n
λn
n
λ 1+p n 1+pwith probability 1 − 4e−τ for τ ≥ 1. This gives the assertion.
If we set

λ = n− (2α+d)(1+p)+(ϵ−p+ϵp) , γ = n− (2α+d)(1+p)+(ϵ−p+ϵp) ,2α+d1

and take ϵ and p suﬃciently small, then we immediately have the following corollary.
Corollary 1. Suppose Assumption 1 is satisﬁed. Then, for all ρ, ρ′ > 0, there exists a
constant K > 0 depending on M, c, ρ, and ρ′ such that, for all n ≥ 1 and τ ≥ 1, thedensity-diﬀerence estimator fb with appropriate choice of γ and λ satisﬁes(
2ατ )∥fb − f ∥2L2 (Rd ) + λ∥fb∥2Hγ ≤ K n− 2α+d +ρ + 1−ρ′nwith probability not less than 1 − 4e−τ .
αNote that n− 2α+d is the optimal learning rate to estimate a function in B2,∞
(Eberts& Steinwart, 2011). Therefore, the density-diﬀerence estimator with a Gaussian kernel
achieves the optimal learning rate by appropriately choosing the regularization parameterand the Gaussian width. Because the learning rate depends on α, the LSDD estimator
has adaptivity to the smoothness of the true function.
Our analysis heavily relies on the techniques developed in Eberts and Steinwart (2011)
for a regression problem. The main diﬀerence is that the analysis in their paper involvesa clipping procedure, which stems from the fact that the analyzed estimator requires an
empirical approximation of the expectation of the square term. The Lipschitz continuity ofthe square function f 7→ f 2 is utilized to investigate this term, and the clipping procedure
is used to ensure the Lipschitz continuity. On the other hand, in the current paper, we
can exactly compute ∥f ∥2L2 (Rd ) so that we do not need the Lipschitz continuity.
2αB

Derivation of Eq.(13)When λ (≥ 0) is small, (H + λI b )−1 can be expanded as(H + λI b )−1 = H −1 − λH −2 + op (λ),

Density-Diﬀerence Estimation

35

where op denotes the probabilistic order. Then Eq.(12) can be expressed asb ⊤θ
b + (1 − β)θ
b⊤ H θbβh⊤⊤b (H + λI b )−1 hb + (1 − β)hb (H + λI b )−1 H (H + λI b )−1 hb= βhbb − λβ hb ⊤ H −2 hb ⊤ H −1 h= βhb + op (λ)
b − 2λ(1 − β)hb ⊤ H −2 hb ⊤ H −1 h
+ (1 − β)hb + op (λ),
b − λ(2 − β)hb ⊤ H −2 hb ⊤ H −1 h=hwhich concludes the proof.

C

Derivation of Eq.(14)b = h, we have
Because E[h]b − h⊤ H −1 h] = E[(hb − h)⊤ H −1 (hb − h)]b ⊤ H −1 hE[h()
b − h)(hb − h)⊤ ]= tr H −1 E[(h
(
(
))
1
1
−1
= tr HV p [ψ] + ′ V p′ [ψ],nnwhich concludes the proof.

References
Anderson, N., Hall, P., & Titterington, D. (1994). Two-sample test statistics for measuring
discrepancies between two multivariate probability density functions using kernel-based
density estimates. Journal of Multivariate Analysis, 50, 41–54.
Arlot, S., Celisse, A., & Harchaoui, Z. (2012). Kernel change-point detection (Technical
Report 1202.3878). arXiv.
Aronszajn, N. (1950). Theory of reproducing kernels. Transactions of the American
Mathematical Society, 68, 337–404.
Atif, J., Ripoche, X., & Osorio, A. (2003). Non-rigid medical image registration by maximisation of quadratic mutual information. IEEE 29th Annual Northeast Bioengineering
Conference (pp. 32–40).
Baranchik, A. J. (1964). Multiple regression and estimation of the mean of a multivariate
normal distribution (Technical Report 51). Department of Statistics, Stanford University, Stanford, CA, USA.

Density-Diﬀerence Estimation

36

Basu, A., Harris, I. R., Hjort, N. L., & Jones, M. C. (1998). Robust and eﬃcient estimation
by minimising a density power divergence. Biometrika, 85, 549–559.
Besbeas, P., & Morgan, B. J. T. (2004). Integrated squared error estimation of normal
mixtures. Computational Statistics & Data Analysis, 44, 517–526.
Bickel, S., Brückner, M., & Scheﬀer, T. (2007). Discriminative learning for diﬀering
training and test distributions. Proceedings of the 24th International Conference on
Machine Learning (ICML2007) (pp. 81–88).
Boser, B. E., Guyon, I. M., & Vapnik, V. N. (1992). A training algorithm for optimal
margin classiﬁers. Proceedings of the Fifth Annual ACM Workshop on Computational
Learning Theory (pp. 144–152). ACM Press.
Chapelle, O., Schölkopf, B., & Zien, A. (Eds.). (2006). Semi-supervised learning. Cambridge, MA, USA: MIT Press.
Chow, C. K. (1970). On optimum recognition error and reject tradeoﬀ. IEEE Transactions
on Information Theory, 16, 41–46.
Cortes, C., Mansour, Y., & Mohri, M. (2010). Learning bounds for importance weighting.
Advances in Neural Information Processing Systems 23 (pp. 442–450).
Cortes, C., & Vapnik, V. (1995). Support-vector networks. Machine Learning, 20, 273–
297.
Desobry, F., Davy, M., & Doncarli, C. (2005). An online kernel change detection algorithm. IEEE Transactions on Signal Processing, 53, 2961–2974.
du Plessis, M. C. (2013). Class-prior change labeling. Unpublished research memo.
du Plessis, M. C., & Sugiyama, M. (2012). Semi-supervised learning of class balance
under class-prior change by distribution matching. Proceedings of 29th International
Conference on Machine Learning (ICML2012) (pp. 823–830). Edinburgh, Scotland.
Duda, R. O., Hart, P. E., & Stork, D. G. (2001). Pattern classiﬁcation. New York, NY,
USA: Wiley. Second edition.
Duong, T., Koch, I., & Wand, M. P. (2009). Highest density diﬀerence region estimation
with application to ﬂow cytometric data. Biometrical Journal, 51, 504–521.
Eberts, M., & Steinwart, I. (2011). Optimal learning rates for least squares SVMs using
Gaussian kernels. Advances in Neural Information Processing Systems 24 (pp. 1539–
1547).
Efron, B., & Tibshirani, R. J. (1993). An introduction to the bootstrap. New York, NY,
USA: Chapman & Hall/CRC.

Density-Diﬀerence Estimation

37

Farrell, R. H. (1972). On the best obtainable asymptotic rates of convergence in estimationof a density function at a point. The Annals of Mathematical Statistics, 43, 170–180.
Gray, D. M., & Principe, J. C. (2010). Quadratic mutual information for dimensionality
reduction and classiﬁcation. Proceedings of SPIE (p. 76960D).
Gretton, A., Smola, A., Huang, J., Schmittfull, M., Borgwardt, K., & Schölkopf, B.
(2009). Covariate shift by kernel mean matching. Dataset Shift in Machine Learning(pp. 131–160). Cambridge, MA, USA: MIT Press.
Hall, P., & Wand, M. P. (1988). On nonparametric discrimination using density diﬀerences. Biometrika, 75, 541–547.
Harchaoui, Z., Bach, F., & Moulines, E. (2009). Kernel change-point analysis. Advancesin Neural Information Processing Systems 21 (pp. 609–616).
Härdle, W., Müller, M., Sperlich, S., & Werwatz, A. (2004). Nonparametric and semiparametric models. Berlin, Germany: Springer.
Hastie, T., Tibshirani, R., & Friedman, J. (2001). The elements of statistical learning:
Data mining, inference, and prediction. New York, NY, USA: Springer.
Henkel, R. E. (1976). Tests of signiﬁcance. Beverly Hills, CA, USA.: SAGE Publication.
Huang, J., Smola, A., Gretton, A., Borgwardt, K. M., & Schölkopf, B. (2007). Correcting
sample selection bias by unlabeled data. Advances in Neural Information Processing
Systems 19 (pp. 601–608). Cambridge, MA, USA: MIT Press.
Kanamori, T., Hido, S., & Sugiyama, M. (2009). A least-squares approach to direct
importance estimation. Journal of Machine Learning Research, 10, 1391–1445.
Kanamori, T., Suzuki, T., & Sugiyama, M. (2012). Statistical analysis of kernel-based
least-squares density-ratio estimation. Machine Learning, 86, 335–367.
Kawahara, Y., & Sugiyama, M. (2012). Sequential change-point detection based on direct
density-ratio estimation. Statistical Analysis and Data Mining, 5, 114–127.
Kawahara, Y., Yairi, T., & Machida, K. (2007). Change-point detection in time-series data
based on subspace identiﬁcation. Proceedings of the 7th IEEE International Conferenceon Data Mining (pp. 559–564).
Kullback, S., & Leibler, R. A. (1951). On information and suﬃciency. The Annals of
Mathematical Statistics, 22, 79–86.
Liu, B., Cheng, H. D., Huang, J., Tian, J., Tang, X., & Liu, J. (2010). Probability density
diﬀerence-based active contour for ultrasound image segmentation. Pattern Recognition,
43, 2028–2042.

Density-Diﬀerence Estimation

38

Liu, S., Yamada, M., Collier, N., & Sugiyama, M. (2013). Change-point detection intime-series data by relative density-ratio estimation. Neural Networks, 43, 72–83.
Moskvina, V., & Zhigljavsky, A. A. (2003). An algorithm based on singular spectrum
analysis for change-point detection. Communication in Statistics: Simulation & Computation, 32, 319–352.
Nguyen, X., Wainwright, M. J., & Jordan, M. I. (2010). Estimating divergence functionals
and the likelihood ratio by convex risk minimization. IEEE Transactions on Information
Theory, 56, 5847–5861.
Parzen, E. (1962). On the estimation of a probability density function and mode. The
Annals of Mathematical Statistics, 33, 1065–1076.
Pérez-Cruz, F. (2008). Kullback-Leibler divergence estimation of continuous distributions.
Proceedings of IEEE International Symposium on Information Theory (pp. 1666–1670).
Toronto, Canada.
Qin, J. (1998). Inferences for case-control and semiparametric two-sample density ratio
models. Biometrika, 85, 619–630.
Rao, C. R. (1965). Linear statistical inference and its applications. New York, NY, USA:
Wiley.
Rifkin, R., Yeo, G., & Poggio, T. (2003). Regularized least-squares classiﬁcation. In
J. A. K. Suykens, G. Horvath, S. Basu, C. Micchelli and J. Vandewalle (Eds.), Advancesin learning theory: Methods, models and applications, vol. 190 of NATO Science Series
III: Computer & Systems Sciences, 131–154. Amsterdam, the Netherlands: IOS Press.
Rockafellar, R. T. (1970). Convex analysis. Princeton, NJ, USA: Princeton University
Press.
Saerens, M., Latinne, P., & Decaestecker, C. (2002). Adjusting the outputs of a classiﬁer
to new a priori probabilities: A simple procedure. Neural Computation, 14, 21–41.
Scott, D. W. (2001). Parametric statistical modeling by minimum integrated square error.
Technometrics, 43, 274–285.
Silva, J., & Narayanan, S. S. (2010). Information divergence estimation based on datadependent partitions. Journal of Statistical Planning and Inference, 140, 3180–3198.
Silverman, B. W. (1986). Density estimation for statistics and data analysis. London,
UK: Chapman and Hall.
Steinwart, I., & Christmann, A. (2008). Support vector machines. New York, NY, USA:
Springer.

Density-Diﬀerence Estimation

39

Steinwart, I., & Scovel, C. (2007). Fast rates for support vector machines using Gaussian
kernels. The Annals of Statistics, 35, 575–607.
Sugiyama, M. (2010). Superfast-trainable multi-class probabilistic classiﬁer by leastsquares posterior ﬁtting. IEICE Transactions on Information and Systems, E93-D,
2690–2701.
Sugiyama, M., & Kawanabe, M. (2012). Machine learning in non-stationary environments: Introduction to covariate shift adaptation. Cambridge, Massachusetts, USA:
MIT Press.
Sugiyama, M., Kawanabe, M., & Chui, P. L. (2010a). Dimensionality reduction for density
ratio estimation in high-dimensional spaces. Neural Networks, 23, 44–59.
Sugiyama, M., Krauledat, M., & Müller, K.-R. (2007). Covariate shift adaptation byimportance weighted cross validation. Journal of Machine Learning Research, 8, 985–
1005.
Sugiyama, M., Suzuki, T., & Kanamori, T. (2012a). Density ratio estimation in machinelearning. Cambridge, UK: Cambridge University Press.
Sugiyama, M., Suzuki, T., & Kanamori, T. (2012b). Density ratio matching under the
Bregman divergence: A uniﬁed framework of density ratio estimation. Annals of theInstitute of Statistical Mathematics, 64, 1009–1044.
Sugiyama, M., Suzuki, T., Kanamori, T., du Plessis, M. C., Liu, S., & Takeuchi, I. (2012c).
Density-diﬀerence estimation. Advances in Neural Information Processing Systems 25.
Sugiyama, M., Suzuki, T., Nakajima, S., Kashima, H., von Bünau, P., & Kawanabe,
M. (2008). Direct importance estimation for covariate shift adaptation. Annals of theInstitute of Statistical Mathematics, 60, 699–746.
Sugiyama, M., Takeuchi, I., Suzuki, T., Kanamori, T., Hachiya, H., & Okanohara, D.
(2010b). Least-squares conditional density estimation. IEICE Transactions on Information and Systems, E93-D, 583–594.
Sugiyama, M., Yamada, M., von Bünau, P., Suzuki, T., Kanamori, T., & Kawanabe, M.
(2011). Direct density-ratio estimation with dimensionality reduction via least-squares
hetero-distributional subspace search. Neural Networks, 24, 183–198.
Suzuki, T., & Sugiyama, M. (2013). Suﬃcient dimension reduction via squared-loss mutual
information estimation. Neural Computation, 25, 725–758.
Takeuchi, Y., & Yamanishi, K. (2006). A unifying framework for detecting outliers andchange points from non-stationary time series data. IEEE Transactions on Knowledge
and Data Engineering, 18, 482–489.

Density-Diﬀerence Estimation

40

Torkkola, K. (2003). Feature extraction by non-parametric mutual information maximization. Journal of Machine Learning Research, 3, 1415–1438.
Vapnik, V. N. (1998). Statistical learning theory. New York, NY, USA: Wiley.
von Bünau, P., Meinecke, F. C., Király, F. J., & Müller, K.-R. (2009). Finding stationary
subspaces in multivariate time series. Physical Review Letters, 103, 214101.
Wang, Q., Kulkarmi, S. R., & Verdú, S. (2005). Divergence estimation of continuous
distributions based on data-dependent partitions. IEEE Transactions on Information
Theory, 51, 3064–3074.
Yamada, M., & Sugiyama, M. (2011). Direct density-ratio estimation with dimensionality
reduction via hetero-distributional subspace analysis. Proceedings of the Twenty-Fifth
AAAI Conference on Artiﬁcial Intelligence (AAAI2011) (pp. 549–554). San Francisco,California, USA: The AAAI Press.
Yamada, M., Suzuki, T., Kanamori, T., Hachiya, H., & Sugiyama, M. (2013). Relative
density-ratio estimation for robust distribution comparison. Neural Computation, 25,1324–1370.
Yamanaka, M., Matsugu, M., & Sugiyama, M. (2013a). Detection of activities and events
without explicit categorization. IPSJ Transactions on Mathematical Modeling and Its
Applications. to appear.
Yamanaka, M., Matsugu, M., & Sugiyama, M. (2013b). Salient object detection based on
direct density-ratio estimation. IPSJ Transactions on Mathematical Modeling and Its
Applications. to appear.