Direct Density Ratio Estimation with Dimensionality Reduction
Masashi Sugiyama∗, Satoshi Hara†, Paul von Bünau‡,
Taiji Suzuki§, Takafumi Kanamori¶, and Motoaki Kawanabe∥
Abstract
Methods for directly estimating the ratio of two probability
density functions without going through density estimation
have been actively explored recently since they can be usedfor various data processing tasks such as non-stationarity
adaptation, outlier detection, conditional density estimation, feature selection, and independent component analysis. However, even the state-of-the-art density ratio estimation methods still perform rather poorly in high-dimensional
problems. In this paper, we propose a new density ratio estimation method which incorporates dimensionality reduction
into a density ratio estimation procedure. Our key idea is to
identify a low-dimensional subspace in which the two densities corresponding to the denominator and the numerator in
the density ratio are signiﬁcantly diﬀerent. Then the density
ratio is estimated only within this low-dimensional subspace.
Through numerical examples, we illustrate the eﬀectivenessof the proposed method.

1 Introduction
Recently, it has been shown [28] that various data
mining and machine learning tasks can be formulated interms of the ratio of two probability density functions
pde (x) and pnu (x):r(x) =

pnu (x),pde (x)weighting the loss function according to the density ratio:∫
E [loss(x)] = loss(x)pnu (x)dx
pnu (x)∫= loss(x)r(x)pde (x)dx=

E [loss(x)r(x)].

pde (x)Thus, the inter-domain bias can be canceled bydensity-ratio weighted learning.
Applications of importance sampling include nonstationarity adaptation [23, 43, 31, 30, 27, 29],transfer learning [26, 41], and multi-task learning[2].
• Outlier detection: Let us consider an outlier detection problem of ﬁnding outliers in an evaluation
dataset based on another “model” dataset that only
contains inliers [11, 10, 40, 24]. Deﬁning the density
ratio over the two sets of samples, one can see thatthe density-ratio values for inliers are close to one,while those for outliers tend to be signiﬁcantly deviated from one. Thus the density-ratio value could
be used as an index of the degree of outlyingness.
The same technique can also be applied to changepoint detection in time-series [17].

where the subscripts ‘nu’ and ‘de’ denote ‘numerator’and ‘denominator’, respectively. Possible usage of thedensity ratio includes the following tasks.
• Importance sampling in supervised learning:
Samples in one domain (drawn from pde (x)) is utilized for learning in other domains (characterized
by pnu (x)). Such data transfer is carried out by∗ Tokyo Institute of Technology and PRESTO, Japan Scienceand Technology Agency (JST)† Tokyo Institute of Technology‡ Technical University of Berlin§ The University of Tokyo¶ Nagoya University∥ Fraunhofer FIRST.IDA

595• Conditional probability estimation: Suppose
we are given n i.i.d. paired samples {(xk , y k )}nk=1
drawn from a joint distribution with density
q(x, y).
The goal is to estimate the conditional probability q(y|x). When the domain of
x is continuous, conditional density estimation is
not straightforward since a naive empirical approximation cannot be used [4]. Let us regard{(xk , y k )}nk=1 as samples corresponding to the numerator of the density ratio and {xk }nk=1 as samples corresponding to the denominator of the density ratio, i.e., we consider the density ratio deﬁned
asq(x, y)= q(y|x),r(x, y) =q(x)Copyright © by SIAM.
Unauthorized reproduction of this article is prohibited.

where q(x) is the marginal density of x. Then a However, model selection methods are not available fordensity-ratio estimation method directly gives an kernel mean matching. Thus, several tuning parameestimator of the conditional density.
ters such as the kernel width and the regularization paThe problem is conditional density estimation when rameter need to be hand-tuned using some heuristics,
y is continuous [34], while it is probabilistic classi- which is highly unreliable in practice. Furthermore, thecomputation of kernel mean matching is rather expenﬁcation when y is categorical.
sive since a quadratic programming problem has to be• Estimation of divergence functionals/mutual solved.
information: Suppose we are given n i.i.d. paired
An alternative method based on logistic regression
samples {(xk , y k )}nk=1 drawn from a joint distri- [22, 6, 3] formulates the density ratio estimation probbution with density q(x, y). Let us denote the lem as the problem of separating samples from the two
marginal densities of x and y by q(x) and q(y), sets by logistic regression. An advantage of the logistic
respectively. Then mutual information I(X, Y ) be- regression formulation is that standard cross-validation
tween random variables X and Y is deﬁned by(CV) is available for model selection since the problem
one needs to solve is a standard supervised classiﬁca∫∫q(x, y)tion problem. Thus, all the tuning parameters can beI(X, Y ) =q(x, y) log
dxdy,
q(x)q(y)objectively determined by CV. However, it is still computationally rather demanding due to non-linearity ofwhich can be used for measuring independence bethe optimization problem.
tween X and Y . Let us regard {(xk , y k )}nk=1 as
Maximum likelihood estimation of density ratio
samples corresponding to the numerator of the denfunctions
is another line of methods that allows us tosity ratio and {(xk , y k′ )}nk,k′ =1 as samples correavoid density estimation [20, 19, 32, 33]. An advantage
sponding to the denominator of the density raof the maximum likelihood approach is that it is alsotio. Then mutual information can be directly esequipped with CV and thus model selection is possible
timated using a density-ratio estimation method[32, 33]. However, this approach is also computationally[20, 19, 32, 33, 38, 13, 39, 14].
rather expensive due to non-linearity of the optimization
Such an independence measure can be used for problem.
various purposes such as variable selection (inputRecently, a least-squares method of density ratio esoutput dependency) [38, 37], supervised dimen- timation was proposed [13, 14]. This is also equippedsionality reduction (input-output dependency) [36], with a build-in CV method, and hence all tuning paand independent component analysis (input-input rameters can be objectively determined. Furthermore,
dependency) [35].
its solution can be computed analytically just by solving
a system of linear equations. Thus it is highly advantaBecause of the wide applicability of density ratios,
geous in terms of computation time. The least-squares
the problem of estimating the density ratios is attracting
method was also shown to be numerically stable [15] una great deal of attention recently and various methods
der condition number analysis. Thus the least-squares
have been explored [22, 6, 12, 3, 19, 32, 33, 13, 14]. A
method is a reliable density ratio estimator.
naive approach is to estimate the two densities in the
As described above, various methods have been deratio (corresponding to the denominator and the nuveloped for directly estimating the density ratios. The
merator) separately using a ﬂexible technique such as
success of these direct density-ratio estimation methods
non-parametric kernel density estimation [8] and thencould be intuitively understood through Vapnik’s printake the ratio of the estimated densities. However, thisciple [42]:naive two-step approach is not reliable in practical situations since kernel density estimation performs poorly“When solving a problem of interest, do notin high-dimensional problems; furthermore, division bysolve a more general problem as an intermedian estimated density tends to magnify the estimation
ate step”.
error.
To improve the estimation accuracy, various meth- The support vector machine would be a successful
ods have been developed for directly estimating the example of this principle—instead of estimating thedensity ratio without going through density estimation. data generation model, it directly models the decision
The moment matching method based on reproducing boundary which is simpler and suﬃcient for patternkernels [1, 25] called kernel mean matching [12] uses recognition. In the current context, estimating the two
the kernel trick to eﬃciently match the mean of two densities is more general than estimating the density
sets of samples in a reproducing kernel Hilbert space. ratio since knowing the two densities implies knowing

596

Copyright © by SIAM.
Unauthorized reproduction of this article is prohibited.

Knowing two densities

pnu (x), pde (x)The problem we address in this paper is to estimate thedensity ratiopnu (x)r(x) =pde (x)

Knowing ratior(x) =

pnu (x)pde (x)nde
nu nnufrom samples {xdei }i=1 and {xj }j=1 .
Our basic idea is to ﬁrst identify a lowerdimensional hetero-distributional subspace in which thetwo densities corresponding to the denominator and thenumerator are signiﬁcantly diﬀerent, and then performdensity ratio estimation only in this subspace.

Figure 1: Density ratio estimation is substantially easier 2.2 Hetero-distributional Subspace.
than density estimation. The density ratio r(x) can be Let u be an m-dimensional vector (1 ≤ m ≤ d) and v
computed if two densities pnu (x) and pde (x) are known. is a (d − m)-dimensional vector deﬁned as
However, even if the density ratio is known, the two[ ] [ ]
densities cannot be computed in general.
u
U
=
x.
v
V
the density ratio, but not vice versa. Thus directly
estimating the density ratio would be more promisingthan density ratio estimation via density estimation(Figure 1). Rigorous theoretical analysis along this line
was carried out in the paper [16].
Although the above density ratio estimators wereshown to compare favorably with naive kernel density
estimation through extensive experiments, density ratio
estimation in high-dimensional problems is still challenging. In this paper, we propose to incorporate dimensionality reduction in a density ratio estimation procedure. More speciﬁcally, our idea is to identify a subspace in which the two densities are signiﬁcantly diﬀerent (called the hetero-distributional subspace); Then we
perform density ratio estimation only in this subspace.
We derive an analytic estimator of a divergence between
the two densities and this estimator is used for searching
the hetero-distributional subspace. Through numericalexamples, we illustrate the usefulness of the proposed
method.
2

U is an m × d matrix and V is a (d − m) × d matrix;
furthermore, without loss of generality, it is possible toassume that the row vectors of U and V form an orthonormal basis, i.e., U and V correspond to “projection” matrices that are orthogonally complementary toeach other (see Figure 2). Using the decomposition ofx into u and v, we can express the two densities pde (x)and pnu (x) as
pde (x) = pde (v|u)pde (u),pnu (x) = pnu (v|u)pnu (u).
Our key theoretical assumption which forms thebasis of our proposed algorithm is that the conditional
densities pde (v|u) and pnu (v|u) agree with each other,i.e., the two densities pde (x) and pnu (x) are decomposedaspde (x) = p(v|u)pde (u),pnu (x) = p(v|u)pnu (u),Problem Formulation

where p(v|u) is the common conditional density. ThisIn this section, we formulate the problem of density ratio assumption implies that the marginal densities of u are
diﬀerent, but the conditional density of v given u isestimation with dimensionality reduction.
common. Then the density ratio is simpliﬁed as
2.1 Density Ratio Estimation.
pnu (u)r(x) = r(u) =.
Let D (⊂ Rd ) be the data domain and suppose wepde (u)are given independent and identically distributed (i.i.d.)ndesamples {xdei }i=1 from a distribution with density Thus, the density ratio does not have to be estimatednnu
pde (x) and i.i.d. samples {xnuj }j=1 from another distri- in the entire d-dimensional space, but only in the mbution with density pnu (x). We assume that the density dimensional subspace.
pde (x) is strictly positive, i.e.,
Let us consider the set of all subspaces such thatthe conditional density p(v|u) is common to pde (x) andpde (x) > 0 for all x ∈ D.
pnu (x). We refer to the intersection of such subspaces

597

Copyright © by SIAM.
Unauthorized reproduction of this article is prohibited.

Hetero-distributional subspace

Figure 2: Hetero-distributional subspace.
as the hetero-distributional subspace. Thus the hetero- Lemma 1.
distributional subspace is the ‘smallest’ subspace outPD[pnu (x), pde (x)] − PD[pnu (u), pde (u)]side which the conditional density p(v|u) is common to)2∫ (pnu (x) pnu (u)pde (x) and pnu (x).
−
(3.1)
=
pde (x)dx
For the moment, we assume that the true dipde (x)pde (u)mensionality m of the hetero-distributional subspace is≥ 0.
known. Later, we explain how m can be estimated from
[Proof:]data in practice.
)2∫ (pnu (x) pnu (u)−pde (x)dx0≤
3 Direct Density Ratio Estimation withpde (x)pde (u)Dimensionality Reduction)2∫ (pnu (x)pnu (u)In this section, we propose a new density ratio estimator=−1−+ 1 pde (x)dx
pde (x)pde (u)which involves dimensionality reduction.
)2∫ (pnu (x)=− 1 pde (x)dx
3.1 Characterizing the Hetero-distributionalpde (x))2∫ (Subspace by the Pearson Divergence.
pnu (u)+−1pde (x)dx
We use the Pearson divergence (PD) as our criterion forpde (u)1
evaluating the discrepancy between two distributions .
)()∫ (pnu (x)pnu (u)PD from pnu (u) to pde (u) is deﬁned and expressed as−2−1− 1 pde (x)dx
pde (x)pde (u))2∫ (= PD[pnu (x), pde (x)] + PD[pnu (u), pde (u)]pnu (u)∫∫PD[pnu (u), pde (u)] =− 1 pde (u)dupnu (u)pde (u)
−
2
p
(x)dx
+
2
pnu (x)dxnu∫
pde (u)pnu (u)∫∫=pnu (u)du − 1.
pde (u)+ 2 pnu (u)du − 2 pde (x)dx
PD[pnu (u), pde (u)] vanishes if and only if pnu (u) =pde (u) for all u.
The following lemma characterizes the heterodistributional subspace in terms of PD.

1 It

is also possible to characterize the hetero-distributional
subspace by the Kullback-Leibler divergence [18]. However, as
shown later, PD allows us to obtain an analytic-form estimatorof the divergence which is useful in hetero-distributional subspace
search.

598= PD[pnu (x), pde (x)] + PD[pnu (u), pde (u)]∫pnu (u)−2pnu (u)du + 2
pde (u)= PD[pnu (x), pde (x)] − PD[pnu (u), pde (u)].
Eq.(3.1) is non-negative and it vanishes if andonly if pnu (v|u) = pde (v|u) for all u, v.
Since
PD[pnu (x), pde (x)] is a constant and does not dependon U , maximizing PD[pnu (u), pde (u)] with respect to U
leads to pnu (v|u) = pde (v|u) for all u, v (see Figure 3).
That is, the hetero-distributional subspace can be characterized as the maximizer of PD[pnu (u), pde (u)].

Copyright © by SIAM.
Unauthorized reproduction of this article is prohibited.

 

pnu x
p nu u
−
pde x
pde u

2

pde x xp nu u , p de u
pnu x , pde x(constant)Figure 3: Since PD[pnu (x), pde (x)] is a constant, minimizing
maximizing PD[pnu (u), pde (u)].

∫ ( pnu (x)pnu (u)pde (x) − pde (u))2pde (x)dx is equivalent to

3.2 Estimation of PD.
where λ (≥ 0) is a regularization parameter and
It is not possible to directly ﬁnd the maximizer ofnde∑dePD[pnu (u), pde (u)] since pde (u) and pnu (u) are unb ℓ,ℓ′ = 1
ψℓ (udeHi )ψℓ′ (ui ),nde i=1
known. According to the Legendre-Fenchel duality [5],we have
nnu
1 ∑
b
hℓ =
ψℓ (unuj ).
nnu j=1(3.2)
PD[pnu (u), pde (u)] = max J(g),g

By setting the derivative of the above objective function
to zero and solving it, we can obtain the maximizer
analytically as

where
(3.3)
J(g) := −∫∫2g(u) pde (u)du + 2

g(u)pnu (u)du − 1.

where ⊤ denotes the transpose of a matrix or a vector
and I b is the b-dimensional identity matrix.
Then an analytic estimator of the Pearson divergence dPD[pnu (u), pde (u)] is given as

Let us employ a parametric modelg(u) =b∑b
c + λI b )−1 h,b = (bαα1 , . . . , αb b )⊤ = (Hαℓ ψℓ (u),ℓ=1where {ψℓ (x)}bℓ=1 are basis functions such that ψℓ (x) ≥0 for all x (∈ D) and for ℓ = 1, . . . , b. In our
experiments, we use the Gaussian kernel model:(3.4)

g(u) =nnu∑

αℓ Kσ (u, unu
ℓ ),

ℓ=1

where

(
)
∥u − u′ ∥2
Kσ (u, u′ ) = exp −
.
2σ 2

d
PD[pnu (u), pde (u)] =

b
∑

α
bℓ b
hℓ − 1.

ℓ=1

We note that the tuning parameters in the above
procedure (i.e., the Gaussian width σ and the regularization parameter λ) can be determined by crossvalidation (CV) over the score function J(g) (seeEq.(3.3)). Using the Sherman-Woodbury-Morrison formula [7], we can actually compute the leave-one-out CV
score analytically, which is computationally very eﬃcient. However, we omit the details.

3.3 Hetero-distributional Subspace Search.
Let us maximize an empirical and regularized variant of
Given
the
Pearsondivergenceestimator
J(g) (see Eq.(3.3)) over the parametric model.
dPD[pnu (u), pde (u)], our next task is to ﬁnd a maximizer
d nu (u), pde (u)] with respect to U and identifyof PD[p
b
b
b
∑
∑
∑
2
b
the
hetero-distributional subspace (cf. Lemma 1).
bmax −αℓ αℓ′ Hℓ,ℓ′ + 2
αℓ hℓ − λαℓ ,{αℓ }bℓ=1
A gradient descent approach would be a standard′ℓ,ℓ =1
ℓ=1
ℓ=1

599

Copyright © by SIAM.
Unauthorized reproduction of this article is prohibited.

Rotation acrossthe subspace

Rotation within
the subspace
Hetero-distributional
subspace

Figure 4: In the hetero-distributional subspace search,rotation which changes the subspace only matters (thesolid arrow); rotation within the subspace (dotted arrow) can be ignored since this does not change the subspace. Similarly, rotation within the orthogonal complement of the hetero-distributional subspace can also
be ignored (not depicted in the ﬁgure).
choice for optimization.
U ←− U + t

∂d
PD
,
∂U

where O d×d′ is the d × d′ matrix with all zeros. For a
skew-symmetric matrix M (∈ Rd×d ), i.e., M ⊤ = −M ,rotation of U can be expressed as
[ ]
[
]
U
I m O m×(d−m) eM
(3.5)
,
V
where eM is the matrix exponential of M ; M = O d×d
corresponds to no rotation. Our idea is not to update
U directly, but through M . The derivative of d
PD withrespect to M at M = O d×d is given by


c
∂ PD
[
]
dd
PD
∂U
 U⊤ V ⊤
=
c
dM
∂ PD
M =O d×d∂V[ ] [(
) (
) ]
U
c ⊤
c ⊤
∂ PD
∂ PD
−
∂U
∂V
V


c
dPD
O m×m
V
dU
,
=
c
PD
V )⊤ O (d−m)×(d−m)−( ddUc

PD
where we used the fact that ddV= O (d−m)×d . Thenthe gradient ascent update rule of M is given by

where t is the step size andM ←− tbb∑∑b ℓ,ℓ′
∂d
PD
∂H
∂b
hℓ
=−
α
bℓ α
bℓ′
+2
α
bℓ
,
∂U
∂U
∂U
ℓ,ℓ′ =1ℓ=1(
nde
b ℓ,ℓ′1 ∑∂ψℓ (ude∂Hi )=ψℓ′ (udei )∂U
nde i=1∂U)dede ∂ψℓ′ (ui )+ ψℓ (ui )
,
∂U

dd
PD
dM

,
M =O d×d

where t is a step size. Then U (and also V ) are
updated by Eq.(3.5). See the paper [21] for the details
of geometric structures.

3.4 Estimating the Density Ratio in Heterodistributional Subspace.
nnu
nu
∑
b
Finally, we estimate the density ratio in the hetero∂ψℓ (uj )
∂ hℓ
1
=
.
distributional
subspace. A notable fact of our algo∂Unnu j=1∂Urithm is that the density ratio estimator in the heteroFor the Gaussian kernel model (3.4) which we use in the distributional subspace has already been obtained during the hetero-distributional subspace search; thus weℓ (u)experiments, ∂ψ∂U
is given by
do not need additional computation. More speciﬁcally,∂ψℓ (u)1
the solution of the variational problem maxg J(g) (see⊤= − 2 (u − cℓ )(x − xnuτ (ℓ) ) ψℓ (u).
(u)∂UσEq.(3.2)) is given by ppnu
[19]. Thus, our ﬁnal solution
de (u)
By the gradient ascent iteration over the m × is simply given byd matrix U , we may ﬁnd a local maximizer ofbd∑PD[pnu (u), pde (u)]. On the other hand, the number of
b x),
r
b
(x)
=
α
bℓ ψℓ (Uparameters to be optimized in the gradient algorithmℓ=1
can be actually reduced in the current setup since weare searching for a subspace—rotation within the sub- where U
b is a projection matrix obtained by the heterospace can be ignored (Figure 4). This idea is explained distributional subspace search algorithm.
below in detail.
The above result implies that if the dimensionality
The matrix U can be expressed as
is not reduced (i.e., m = d), the proposed method
[ ]
[
] U
agrees with the density ratio estimator proposed in the
U = I m O m×(d−m),
papers [13, 14]. Thus, the proposed method could be
V

600

Copyright © by SIAM.
Unauthorized reproduction of this article is prohibited.

nnu
de nded
Input: Two sets of samples {xnui }i=1 and {xj }j=1 on R
Output: Density ratio estimator rb(x)For each reduced dimensionality m = 1, . . . , dInitialize embedding matrix U m (∈ Rm×d );Repeat until U m converges
Choose Gaussian width σ and regularization parameter λ by CV;
Update U by the gradient method (see Section 3.3);endb m and corresponding density ratio estimator rbm (x);Obtain embedding matrix U
Compute its CV value as a function of m;
end
Choose the best reduced dimension mb based on the CV score;
Set rb(x) = rbm
b (x);
Figure 5: Pseudo code of the proposed density-ratio estimation algorithm.
regarded as a natural extension of the existing density Dataset (b) (Figure 7(a) and Figure 7(b)):
ratio estimator.
p(v|u) = N (v; u, 12 ),The dimensionality of the hetero-distributional subspace may be chosen by the CV score used for optimizpnu (u) = N (u; 0, 12 ),ing the Gaussian width σ and the regularization parampde (u) = 0.5N (u; −2, 1.52 ) + 0.5N (u; 2, 1.52 ).
eter λ. The entire procedure is summarized in Figure 5.
The true and estimated hetero-distributional sub4 Numerical Examples
spaces are depicted by the dashed and solid lines inIn this section, we illustrate the behavior of the proFigure 6(c) and Figure 7(c). These plots show that
posed method through numerical examples.
the proposed method gives good estimates of the true
hetero-distributional subspace. In Figure 6(e) and Fig4.1 Illustrative Examples.
ure 7(e), density-ratio functions estimated without diLet us consider two-dimensional examples (i.e., d = 2) mensionality reduction by the baseline method proposedand suppose that the two distributions pnu (x) and in the papers [13, 14] are depicted. In Figure 6(f) andpde (x) are diﬀerent only in the one-dimensional sub- Figure 7(f), density-ratio functions estimated with dispace (i.e., m = 1) spanned by (1, 0)⊤ :mensionality reduction by the proposed method are depicted. Compared with the true density ratio functions
depicted in Figure 6(d) and Figure 7(d), we can obx = (x(1) , x(2) )⊤ = (u, v)⊤ ,serve that the proposed method captures the redundantpnu (x) = p(v|u)pnu (u),structure of the true density ratio functions appropripde (x) = p(v|u)pde (u).
ately. Consequently, the propose method gives muchbetter estimates of the density ratio functions than the
Let nnu = nde = 1000. We use the following two baseline method. This illustrates the usefulness of dimensionality reduction in density ratio estimation.
datasets.
4.2 Performance Evaluation using Artificial
Datasets.
Next, we systematically investigate the behavior of thep(v|u) = p(v) = N (v; 0, 12 ),proposed method for high-dimensional data.
For the two datasets used in the previous experpnu (u) = 0.5N (u; −2, 1.52 ) + 0.5N (u; 2, 1.52 ),iments,we increase the entire dimensionality as d =pde (u) = N (u; 0, 12 ),2, 3, . . . , 10 by adding dimensions consisting of standard normal noise. The dimensionality of the heterowhere N (u; µ, σ 2 ) denotes the Gaussian density distributional subspace is estimated based on CV (seewith mean µ and variance σ 2 with respect to u.
Section 3.4).

Dataset (a) (Figure 6(a) and Figure 6(b)):

601

Copyright © by SIAM.
Unauthorized reproduction of this article is prohibited.

0.15
5

0.1

0.04

5

0.02

0.05
0

0
−10

−5

0

5
x(1)

10 −5

0
−10

x

0
−5

(2)

0
x(1)

(a) pnu (x)

(2)

5

10 −5

x

(2)

(b) pde (x)

x de
nu
x

0

x

5

3
5

2
1
0
−10

True Subspace
Estimated Subspace
−5
−10

−5

0
x (1)

5

0
−5

0

5

10

x(1)

(c) Hetero-distributional subspace

10 −5

x(2)

(d) r(x)

4
3
5

2

5

2
1

0
−10

0
−5

0

5
x(1)

10 −5

x

0
−10
(2)

0
−5

0

5
x(1)

(e) rb(x) by baseline method.

10 −5

x(2)

(f) rb(x) by proposed method.

Figure 6: Numerical results for dataset (a).

602

Copyright © by SIAM.
Unauthorized reproduction of this article is prohibited.

0.06

0.15
10

0.1

0.04

10

0.02

0.05
0

0
−10

−5

0

5
x(1)

10 −10

x

0
−10

0
−5

(2)

0
x(1)

(a) pnu (x)105

10 −10x(2)(b) pde (x)

de

x
nu
x

(2)

5

x

0

3
10

2
1

−5

0
−10

True Subspace
Estimated Subspace
−10
−10

−5

0
x (1)

5

0
−5

0

5

10

x(1)

(c) Hetero-distributional subspace

10 −10x(2)(d) r(x)

4
3
10

2

10

2
1

0
−10

0
−5

0

5
x(1)

10 −10

0
−10

x(2)

0
−5

0

5
x(1)

(e) rb(x) by baseline method.

10 −10x(2)(f) rb(x) by proposed method.

Figure 7: Numerical results for dataset (b).

603

Copyright © by SIAM.
Unauthorized reproduction of this article is prohibited.

0.45

0.35
Baseline Method
Proposed Method

0.4

Baseline Method
Proposed Method

0.3

0.35
0.25
0.2

0.25

Error

Error

0.3

0.2

0.15

0.15
0.1
0.1
0.05

0.05
0

2

3

4
5
6
7
8
Entire Data Dimensionality d

9

0

10

2

3

4
5
6
7
8
Entire Data Dimensionality d

(a)

9

10

(b)

Figure 8: Density ratio estimation error (4.6) averaged over 10 runs as a function of the entire data dimensionality
d for the artiﬁcial datasets. The best method in terms of the mean error and comparable methods according to
the t-test at the signiﬁcance level 1% are speciﬁed by ‘◦’; otherwise methods are speciﬁed by ‘×’.
We evaluate the error of a density ratio estimatorrb(x) by
∫
1
2
(4.6)
Error :=
(br(x) − r(x)) pde (x)dx.
2

such as non-stationarity adaptation, outlier detection,
feature selection, and independent component analysis. Improving the computational eﬃciency of heterodistributional subspace search is another important issue to be further investigated.

Figure 8 shows the density ratio estimation error averaged over 10 runs as functions of the entire input dimensionality d. The best method in terms of the mean error
and comparable methods according to the t-test [9] atthe signiﬁcance level 1% are speciﬁed by ‘◦’; otherwise
methods are speciﬁed by ‘×’. This shows that, while theerror of the baseline method without dimensionality reduction increases rapidly as the entire dimensionality dincreases, that of the proposed method is kept moderate. Consequently, the proposed method consistently
outperforms the baseline method.

Acknowledgments

5 Conclusions
The density ratio is becoming a quantity of interestin the machine learning and data mining communities
since it can be used for solving various data processing
tasks. In this paper, we tackled a challenging problem of
estimating density ratios in high-dimensional spaces andgave a new procedure. Our key idea was to estimate theratio only in a subspace in which two distributions (corresponding to the denominator and numerator of thedensity ratio) are signiﬁcantly diﬀerent. The proposed
method was shown to be promising in experiments.
Our future work includes the application of the
proposed method to various data processing tasks

604This work has been supported by SCAT, AOARD, andthe JST PRESTO program.
References[1] N. Aronszajn. Theory of reproducing kernels. Transactions of the American Mathematical Society, 68:337–
404, 1950.
[2] S. Bickel, J. Bogojeska, T. Lengauer, and T. Scheffer. Multi-task learning for HIV therapy screening.
In A. McCallum and S. Roweis, editors, Proceedingsof 25th Annual International Conference on Machine
Learning (ICML2008), pages 56–63, Helsinki, Finland,
Jul. 5–9 2008. Omnipress.
[3] S. Bickel, M. Brückner, and T. Scheﬀer. Discriminative
learning for diﬀering training and test distributions.
In Proceedings of the 24th International Conference on
Machine Learning, pages 81–88, 2007.
[4] C. M. Bishop. Pattern Recognition and Machine
Learning. Springer, New York, NY, USA, 2006.
[5] S. Boyd and L. Vandenberghe. Convex Optimization.
Cambridge University Press, Cambridge, 2004.
[6] K. F. Cheng and C. K. Chu. Semiparametric density
estimation under a two-sample density ratio model.
Bernoulli, 10(4):583–604, 2004.

Copyright © by SIAM.
Unauthorized reproduction of this article is prohibited.

[7] G. H. Golub and C. F. Van Loan. Matrix Computations. Johns Hopkins University Press, Baltimore, MD,
1996.
[8] W. Härdle, M. Müller, S. Sperlich, and A. Werwatz.
Nonparametric and Semiparametric Models. Springer,
Berlin, 2004.
[9] R. E. Henkel. Tests of Significance. SAGE Publication,
Beverly Hills, 1979.
[10] S. Hido, Y. Tsuboi, H. Kashima, M. Sugiyama, and
T. Kanamori. Statistical outlier detection using direct
density ratio estimation. Knowledge and Information
Systems. to appear.
[11] S. Hido, Y. Tsuboi, H. Kashima, M. Sugiyama, and
T. Kanamori. Inlier-based outlier detection via direct
density ratio estimation. In F. Giannotti, D. Gunopulos, F. Turini, C. Zaniolo, N. Ramakrishnan, andX. Wu, editors, Proceedings of IEEE International
Conference on Data Mining (ICDM2008), pages 223–
232, Pisa, Italy, Dec. 15–19 2008.
[12] J. Huang, A. Smola, A. Gretton, K. M. Borgwardt,
and B. Schölkopf. Correcting sample selection bias
by unlabeled data. In B. Schölkopf, J. Platt, andT. Hoﬀman, editors, Advances in Neural Information
Processing Systems 19, pages 601–608. MIT Press,
Cambridge, MA, 2007.
[13] T. Kanamori, S. Hido, and M Sugiyama. Eﬃcient direct density ratio estimation for non-stationarity adaptation and outlier detection. In D. Koller, D. Schuurmans, Y. Bengio, and L. Botton, editors, Advances in
Neural Information Processing Systems 21, pages 809–
816, Cambridge, MA, 2009. MIT Press.
[14] T. Kanamori, S. Hido, and M Sugiyama. A leastsquares approach to direct importance estimation.
Journal of Machine Learning Research, 10:1391–1445,
Jul. 2009.
[15] T. Kanamori, T. Suzuki, and M. Sugiyama. Condition number analysis of kernel-based density ratio estimation. Technical Report TR09-0006, Department ofComputer Science, Tokyo Institute of Technology, Feb.
2009.
[16] T. Kanamori, T. Suzuki, and M. Sugiyama. Theoretical analysis of density ratio estimation. IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences, 2010. to appear.
[17] Y. Kawahara and M. Sugiyama. Change-point detection in time-series data by direct density-ratio estimation. In H. Park, S. Parthasarathy, H. Liu, and
Z. Obradovic, editors, Proceedings of 2009 SIAM International Conference on Data Mining (SDM2009),
pages 389–400, Sparks, Nevada, USA, Apr. 30–May 2
2009.
[18] S. Kullback and R. A. Leibler. On information andsuﬃciency. Annals of Mathematical Statistics, 22:79–
86, 1951.
[19] X. Nguyen, M. Wainwright, and M. Jordan. Estimating divergence functionals and the likelihood ratio by
penalized convex risk minimization. In J. C. Platt,
D. Koller, Y. Singer, and S. Roweis, editors, Advances

605

in Neural Information Processing Systems 20, pages
1089–1096. MIT Press, Cambridge, MA, 2008.
[20] X. Nguyen, M. J. Wainwright, and M. I. Jordan. Nonparametric estimation of the likelihood ratio and divergence functionals. In Proceedings of IEEE International Symposium on Information Theory, pages 2016–2020, Nice, France, 2007.
[21] M. D. Plumbley.
Geometrical methods for nonnegative ICA: Manifolds, Lie groups and toral subalgebras. Neurocomputing, 67(Aug.):161–197, 2005.
[22] J. Qin. Inferences for case-control and semiparametric two-sample density ratio models. Biometrika,
85(3):619–639, 1998.
[23] H. Shimodaira. Improving predictive inference undercovariate shift by weighting the log-likelihood function. Journal of Statistical Planning and Inference,90(2):227–244, 2000.
[24] A. Smola, L. Song, and C. H. Teo. Relative novelty detection. In D. van Dyk and M. Welling, editors, Twelfth
International Conference on Artificial Intelligence andStatistics, volume 5 of JMLR Workshop and Conference Proceedings, pages 536–543, 2009.
[25] I. Steinwart. On the inﬂuence of the kernel on theconsistency of support vector machines. Journal of
Machine Learning Research, 2:67–93, 2001.
[26] A. Storkey and M. Sugiyama. Mixture regression forcovariate shift. In B. Schölkopf, J. C. Platt, andT. Hoﬀmann, editors, Advances in Neural Information
Processing Systems 19, pages 1337–1344, Cambridge,
MA, 2007. MIT Press.
[27] M. Sugiyama, B. Blankertz, M. Krauledat, G. Dornhege, and K.-R. Müller. Importance-weighted crossvalidation for covariate shift. In K. Franke, K.-R.
Müller, B. Nickolay, and R. Schäfer, editors, Pattern
Recognition, volume 4174 of Lecture Notes in Computer
Science, pages 354–363, Berlin, 2006. Springer.
[28] M. Sugiyama, T. Kanamori, T. Suzuki, S. Hido, J. Sese,
I. Takeuchi, and L. Wang. A density-ratio frameworkfor statistical data processing. IPSJ Transactions on
Computer Vision and Applications, 1:183–208, 2009.
[29] M. Sugiyama, M. Krauledat, and K.-R. Müller. Covariate shift adaptation by importance weighted cross validation. Journal of Machine Learning Research, 8:985–
1005, May 2007.
[30] M. Sugiyama and K.-R. Müller. Input-dependent
estimation of generalization error under covariate shift.
Statistics & Decisions, 23(4):249–279, 2005.
[31] M. Sugiyama and K.-R. Müller. Model selection undercovariate shift. In W. Duch, J. Kacprzyk, E. Oja,and S. Zadrozny, editors, Artificial Neural Networks:
Formal Models and Their Applications, volume 3697
of Lecture Notes in Computer Science, pages 235–240,
Berlin, 2005. Springer.
[32] M. Sugiyama, S. Nakajima, H. Kashima, P. von Bünau,
and M. Kawanabe. Direct importance estimation withmodel selection and its application to covariate shift
adaptation. In J. C. Platt, D. Koller, Y. Singer, and
S. Roweis, editors, Advances in Neural Information

Copyright © by SIAM.
Unauthorized reproduction of this article is prohibited.

Processing Systems 20, pages 1433–1440, Cambridge,
MA, 2008. MIT Press.
[33] M. Sugiyama, T. Suzuki, S. Nakajima, H. Kashima,P. von Bünau, and M. Kawanabe. Direct importanceestimation for covariate shift adaptation. Annals of
the Institute of Statistical Mathematics, 60(4):699–746,
2008.
[34] M. Sugiyama, I. Takeuchi, T. Suzuki, T. Kanamori,
H. Hachiya, and D. Okanohara. Least-squares conditional density estimation. IEICE Transactions on Information and Systems, E93-D(3), 2010. to appear.
[35] T. Suzuki and M. Sugiyama. Estimating squaredloss mutual information for independent component
analysis. In T. Adali, C. Jutten, J. M. T. Romano,
and A. K. Barros, editors, Independent Component
Analysis and Signal Separation, volume 5441 of Lecture
Notes in Computer Science, pages 130–137, Berlin,
2009. Springer.
[36] T. Suzuki and M. Sugiyama. Suﬃcient dimension reduction via squared-loss mutual information estimation. Technical Report TR09-0005, Department ofComputer Science, Tokyo Institute of Technology, Feb.
2009.
[37] T. Suzuki, M. Sugiyama, T. Kanamori, and J. Sese.
Mutual information estimation reveals global associations between stimuli and biological processes. BMC
Bioinformatics, 10(1):S52, 2009.
[38] T. Suzuki, M. Sugiyama, J. Sese, and T. Kanamori.
Approximating mutual information by maximum likelihood density ratio estimation. In Y. Saeys, H. Liu,
I. Inza, L. Wehenkel, and Y. Van de Peer, editors,
JMLR Workshop and Conference Proceedings, volume 4 of New Challenges for Feature Selection in Data
Mining and Knowledge Discovery, pages 5–20, 2008.
[39] T. Suzuki, M. Sugiyama, and T. Tanaka. Mutual
information approximation via maximum likelihood
estimation of density ratio. In Proceedings of 2009
IEEE International Symposium on Information Theory(ISIT2009), pages 463–467, Seoul, Korea, Jun. 28–
Jul. 3 2009.
[40] M. Takimoto, M. Matsugu, and M. Sugiyama. Visualinspection of precision instruments by least-squares
outlier detection. In Proceedings of The Fourth International Workshop on Data-Mining and Statistical Science (DMSS2009), pages 22–26, Kyoto, Japan, Jul. 7–8
2009.
[41] Y. Tsuboi, H. Kashima, S. Hido, S. Bickel, and
M. Sugiyama. Direct density ratio estimation for largescale covariate shift adaptation. Journal of Information Processing, 17:138–155, 2009.
[42] V. N. Vapnik. Statistical Learning Theory. Wiley, New
York, 1998.
[43] B. Zadrozny. Learning and evaluating classiﬁers undersample selection bias. In Proceedings of the TwentyFirst International Conference on Machine Learning,
pages 903–910, New York, NY, 2004. ACM Press.

606

Copyright © by SIAM.
Unauthorized reproduction of this article is prohibited.