Neural Networks, vol.24, no.2, pp.183–198, 2011.

Direct Density-ratio Estimation
with Dimensionality Reduction via
Least-squares Hetero-distributional Subspace Search
Masashi Sugiyama
Tokyo Institute of Technology
and PRESTO, Japan Science and Technology Agency,
2-12-1 O-okayama, Meguro-ku, Tokyo 152-8552, Japan.
sugi@cs.titech.ac.jp http://sugiyama-www.cs.titech.ac.jp/˜sugi
Makoto Yamada
Tokyo Institute of Technology
2-12-1 O-okayama, Meguro-ku, Tokyo 152-8552, Japan.
yamada@sg.cs.titech.ac.jp
Paul von Bünau
Technical University of Berlin
Franklinstr. 28/29, 10587 Berlin, Germany.
buenau@cs.tu-berlin.de
Taiji Suzuki
The University of Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan.
s-taiji@stat.t.u-tokyo.ac.jp
Takafumi Kanamori
Nagoya University
Furocho, Chikusaku, Nagoya 464-8603, Japan.
kanamori@is.nagoya-u.ac.jp
Motoaki Kawanabe
Fraunhofer FIRST.IDA
Kekuléstr. 7, 12489 Berlin, Germany.
motoaki.kawanabe@ﬁrst.fraunhofer.de1Direct Density-ratio Estimation with Dimensionality Reduction

2

Abstract
Methods for directly estimating the ratio of two probability density functions havebeen actively explored recently since they can be used for various data processing
tasks such as non-stationarity adaptation, outlier detection, and feature selection. In
this paper, we develop a new method which incorporates dimensionality reduction
into a direct density-ratio estimation procedure. Our key idea is to ﬁnd a lowdimensional subspace in which densities are signiﬁcantly diﬀerent and perform density ratio estimation only in this subspace. The proposed method, D3 -LHSS (Direct
Density-ratio estimation with Dimensionality reduction via Least-squares Heterodistributional Subspace Search), is shown to overcome the limitation of baseline
methods.
Keywords
density ratio estimation, dimensionality reduction, unconstrained least-squares importance ﬁtting

1

Introduction

Recently, it has been demonstrated that various machine learning and data mining tasks
can be formulated in terms of the ratio of two probability density functions (Sugiyama
et al., 2009; Sugiyama et al., 2011). Examples of such tasks include covariate shift adaptation (Shimodaira, 2000; Zadrozny, 2004; Sugiyama et al., 2007; Sugiyama & Kawanabe,2010), transfer learning (Storkey & Sugiyama, 2007), multi-task learning (Bickel et al.,2008), outlier detection (Hido et al., 2008; Smola et al., 2009; Hido et al., 2010), conditional density estimation (Sugiyama et al., 2010c), probabilistic classification (Sugiyama,2010), variable selection (Suzuki et al., 2009a), independent component analysis (Suzuki& Sugiyama, 2009), supervised dimensionality reduction (Suzuki & Sugiyama, 2010), andcausal inference (Yamada & Sugiyama, 2010), For this reason, estimating the density
ratio has been attracting a great deal of attention, and various approaches have been
explored (Silverman, 1978; Ćwik & Mielniczuk, 1989; Gijbels & Mielniczuk, 1995; Sun &Woodroofe, 1997; Jacob & Oliveira, 1997; Qin, 1998; Cheng & Chu, 2004; Huang et al.,2007; Bensaid & Fabre, 2007; Bickel et al., 2007; Sugiyama et al., 2008; Kanamori et al.,2009a; Chen et al., 2009; Sugiyama et al., 2010b; Nguyen et al., 2010).
A naive approach to density ratio estimation is to approximate the two densities in
the ratio (i.e., the numerator and the denominator) separately using a ﬂexible techniquesuch as non-parametric kernel density estimation (Silverman, 1986; Härdle et al., 2004),
and then take the ratio of the estimated densities. However, this naive two-step approachis not reliable in practical situations since kernel density estimation performs poorly inhigh-dimensional cases; furthermore, division by an estimated density tends to magnify
the estimation error. To improve the estimation accuracy, various methods have beendeveloped for directly estimating the density ratio without going through density estimation, e.g., the moment matching method using reproducing kernels (Aronszajn, 1950;

Direct Density-ratio Estimation with Dimensionality Reduction

Knowing two densities

pnu (x), pde (x)3Knowing ratior(x) =

pnu (x)pde (x)Figure 1: Density ratio estimation is substantially easier than density estimation. The
density ratio r(x) can be computed if two densities pnu (x) and pde (x) are known. However,even if the density ratio is known, the two densities cannot be computed in general.
Steinwart, 2001) called kernel mean matching (KMM) (Huang et al., 2007; QuiñoneroCandela et al., 2009), the method based on logistic regression (LR) (Qin, 1998; Cheng& Chu, 2004; Bickel et al., 2007), the distribution matching method under the KullbackLeibler (KL) divergence (Kullback & Leibler, 1951) called the KL importance estimation
procedure (KLIEP) (Sugiyama et al., 2008; Nguyen et al., 2010), and the density-ratio
matching methods under the squared-loss called least-squares importance fitting (LSIF)
and unconstrained LSIF (uLSIF) (Kanamori et al., 2009a). These methods have beenshown to compare favorably with naive kernel density estimation through extensive experiments.
The success of these direct density-ratio estimation methods could be intuitively understood through Vapnik’s principle (Vapnik, 1998): “When solving a problem of interest,
one should not solve a more general problem as an intermediate step”. The support vector machine would be a successful example following this principle—instead of estimating
the data generation model, it directly models the decision boundary which is simpler andsuﬃcient for pattern recognition. In the current context, estimating the densities is moregeneral than estimating the density ratio since knowing the two densities implies knowing
the ratio, but not vice versa (Figure 1). Thus directly estimating the density ratio would
be more promising than density ratio estimation via density estimation.
However, density ratio estimation in high-dimensional cases is still challenging evenwhen the ratio is estimated directly without going through density estimation. Recently,
an approach called Direct Density-ratio estimation with Dimensionality reduction (D3 )has been proposed (Sugiyama et al., 2010a). The basic idea of D3 is the following twostep procedure: First a subspace in which the numerator and denominator densities are
signiﬁcantly diﬀerent (called the hetero-distributional subspace) are identiﬁed, and thendensity ratio estimation is performed in this subspace. The rationale behind this approachis that, in practice, the distribution change does not occur in the entire space, but isoften conﬁned in a subspace. For example, in non-stationarity adaptation scenarios, thedistribution change often occurs only for some attributes and other variables are stable; inDirect Density-ratio Estimation with Dimensionality Reduction

4

outlier detection scenarios, only a small number of attributes would cause a data sample
to be an outlier.
In the D3 algorithm, the hetero-distributional subspace is identiﬁed by searching asubspace in which samples drawn from the two distributions (i.e., the numerator and thedenominator of the ratio) are separated from each other—this search is carried out ina computationally eﬃcient manner using a supervised dimensionality reduction method
called local Fisher discriminant analysis (LFDA) (Sugiyama, 2007). Then, within theidentiﬁed hetero-distributional subspace, a direct density-ratio estimation method called
unconstrained least-squares importance Fitting (uLSIF)—which was shown to be computationally eﬃcient (Kanamori et al., 2009a) and numerically stable (Kanamori et al.,2009b)—is employed for obtaining the ﬁnal density-ratio estimator. Through experiments, this D3 procedure (which we refer to as D3 -LFDA/uLSIF) was shown to improve
the performance in high-dimensional cases.
Although the framework of D3 is promising, the above D3 -LFDA/uLSIF method possesses two fundamental weaknesses: the restrictive deﬁnition of the hetero-distributional
subspace and the limiting ability of its search method. More speciﬁcally, the component
inside the hetero-distributional subspace and its complementary component are assumed
to be statistically independent in the original formulation (Sugiyama et al., 2010a). However, this assumption is rather restrictive and may not be fulﬁlled in practice. Also, in
the above D3 procedure, the hetero-distributional subspace is identiﬁed by searching asubspace in which samples drawn from the numerator and denominator distributions are
separated from each other. If samples from the two distributions are separable, the two
distributions would be signiﬁcantly diﬀerent. However, the opposite may not be always
true, i.e., non-separability does not necessarily imply that the two distributions are different (consider two similar distributions with the common support). Thus LFDA (andany other supervised dimensionality reduction methods) does not necessarily identify thecorrect hetero-distributional subspace.
The goal of this paper is to give a new procedure of D3 that can overcome the above
weaknesses. First, we adopt a more general deﬁnition of the hetero-distributional subspace. More precisely, we remove the independence assumption between the component
inside the hetero-distributional subspace and its complementary component. This allows
us to apply the concept of D3 to a wider class of problems. However, this general definition in turn makes the problem of searching the hetero-distributional subspace morechallenging—supervised dimensionality reduction methods for separating samples drawn
from the two distributions cannot be used anymore, but we need an alternative methodthat identiﬁes the largest subspace such that the two conditional distributions are equivalent in its complementary subspace.
We prove that the hetero-distributional subspace can be identiﬁed by ﬁnding a subspace in which two marginal distributions are maximally diﬀerent under the Pearson
divergence, which is a squared-loss variant of the Kullback-Leibler divergence and is an
instance of the f -divergences (Ali & Silvey, 1966; Csiszár, 1967). Then we propose a
new method, which we call Least-squares Hetero-distributional Subspace Search (LHSS),
for searching a subspace such that the Pearson divergence between two marginal distri-Direct Density-ratio Estimation with Dimensionality Reduction

Goal: Estimate density ratio r x

pnu x
pde x

∼  Naïve density estimation
nu
{{xnu
}ni=1
i
nde
{xde
j }j=1{{xnu}}nnu i.i.d.
∼
i=1 ∼ pnu xifrom samples
i.i.d.
nde{xdej }j=1 ∼ pde x}

Estimating
Estimatingtwo
two
densities
densitiesby
byKDE
KDE

TakingTakingthetheratioratio∼  Direct density-ratio estimation (the above two steps are merged into a single process)
nu
{{xnu
}ni=1
i
Directly
Directlyestimating
estimatingthe
theratio
ratiobyby
KMM,
LogReg,KLIEP,
LSIF,
KMM, LogReg, KLIEP, LSIF,ororuLSIF
uLSIF

nde
{xde
j }j=1∼  Direct density-ratio estimation with dimensionality reduction (D3-LFDA/uLSIF)
nu
{{xnu
}ni=1
i
Identifying hetero-distributional
Directly estimating the ratio by uLSIFnde{xdej }j=1

Identifying hetero-distributional
subspacesubspaceby
byLFDA
LFDA

5

Directly estimating the ratio by uLSIF
ininhetero-distributional
hetero-distributionalsubspacesubspace

r x

}
r x

}
r x

 Proposed approach: Direct density-ratio estimation with dimensionality reduction(D3-LHSS; the above two steps are merged into a single process)
∼
nu
{{xnu
}ni=1
i
nde
{xde
j }j=1

Simultaneously
Simultaneouslyidentifying
identifyinghetero-distributional
hetero-distributionalsubspacesubspaceandanddirectly
directlyestimating
estimatingthe
theratio
ratioininthe
thesubspacesubspaceby
byuLSIF
uLSIF

}
r x

Figure 2: Existing and proposed density-ratio estimation approaches.
butions are maximized. An advantage of the LHSS method is that the subspace search(divergence estimation within a subspace) is carried out also using the density-ratio estimation method uLSIF. Thus the two steps in the D3 procedure (ﬁrst identifying thehetero-distributional subspace and then estimating the density ratio within the subspace)are merged into a single step. Thanks to this, the ﬁnal density-ratio estimator can be automatically obtained without additional computation. We call the combined single-shot
density-ratio estimation procedure D3 via LHSS (D3 -LHSS). Through experiments, weshow that the weaknesses of the existing approach can be successfully overcome by theD3 -LHSS approach.
Relation among the existing and proposed density-ratio estimation methods is summarized in Figure 2.

2

Formulation of Density-ratio Estimation Problem

In this section, we formulate the problem of density ratio estimation and review a relevant
density-ratio estimation method. We brieﬂy summarize possible usage of density ratios
in various data processing tasks in Appendix A.

Direct Density-ratio Estimation with Dimensionality Reduction

2.1

6

Problem Formulation

Let D (⊂ Rd ) be the data domain and suppose we are given independent and identinnu
cally distributed (i.i.d.) samples {xnui }i=1 from a distribution with density pnu (x) andndei.i.d. samples {xdej }j=1 from another distribution with density pde (x). We assume thatthe latter density pde (x) is strictly positive, i.e.,pde (x) > 0 for all x ∈ D.
The problem we address in this paper is to estimate the density ratior(x) :=

pnu (x)pde (x)nnu
de ndefrom samples {xnui }i=1 and {xj }j=1 . The subscripts ‘nu’ and ‘de’ denote ‘numerator’ and‘denominator’, respectively.

2.2

Directly Estimating Density Ratios by Unconstrained
Least-squares Importance Fitting (uLSIF)As described in Appendix A, density ratios are useful in various data processing tasks.
Since the density ratio is usually unknown and needs to be estimated from data, methodsof estimating the density ratio have been actively explored recently (Qin, 1998; Cheng &Chu, 2004; Huang et al., 2007; Bickel et al., 2007; Sugiyama et al., 2008; Kanamori et al.,2009a). Here, we brieﬂy review a direct density-ratio estimation method called unconstrained least-squares importance fitting (uLSIF) proposed by Kanamori et al. (2009a).
For convenience in later sections, we replace the symbol x with u, i.e., let us consider the
problem of estimating the density ratior(u) :=

pnu (u)pde (u)nnu
de ndefrom the i.i.d. samples {unui }i=1 and {uj }j=1 .

2.2.1

Linear Least-squares Estimation of Density Ratios

Let us model the density ratio r(u) by the following linear model:rb(u) :=b∑αℓ ψℓ (u),

ℓ=1

where
α := (α1 , α2 , . . . , αb )⊤Direct Density-ratio Estimation with Dimensionality Reduction

7

are parameters to be learned from data samples, b denotes the number of parameters, ⊤denotes the transpose of a matrix or a vector, and {ψℓ (u)}bℓ=1 are basis functions suchthat
ψℓ (u) ≥ 0 for all u and for ℓ = 1, 2, . . . , b.
nnu
de ndeNote that b and {ψℓ (u)}bℓ=1 could be dependent on the samples {unui }i=1 and {uj }j=1 ,meaning that kernel models are also allowed. We explain how the basis functions{ψℓ (u)}bℓ=1 are designed in Section 2.2.2.
The parameters {αℓ }bℓ=1 in the model rb(u) are determined so that the following squared
error J0 is minimized:∫
1J0 (α) :=(br(u) − r(u))2 pde (u)du
2
∫
∫
∫
1
1
2
=
rb(u) pde (u)du − rb(u)pnu (u)du +r(u)pnu (u)du,2
2

where the last term is a constant and therefore can be safely ignored. Let us denote theﬁrst two terms by J:∫∫
12
rb(u) pde (u)du − rb(u)pnu (u)du.
(1)
J(α) :=
2
Note that the same objective function can be obtained via the Legendre-Fenchel dualityof a divergence (Nguyen et al., 2010).
Approximating the expectations in J by empirical averages, we obtain
b
J(α)
:=

nde
nnu
1 ∑
1 ∑
de 2
rb(unurb(uj ) −i )2nde j=1
nnu i=1

1 cb ⊤ α,= α⊤ Hα−h2c is the b × b matrix with the (ℓ, ℓ′ )-th element
where Hde∑deb ℓ,ℓ′ := 1Hψℓ (udej )ψℓ′ (uj ),nde j=1n(2)b is the b-dimensional vector with the ℓ-th element
and h
nnu
1 ∑
b
hℓ :=
ψℓ (unui ).
nnu i=1

Now the optimization problem is formulated as follows.
[
]
⊤
λ ⊤
1 ⊤c
b
b := argmin α Hα − h α + α α ,
α
2
2
α∈Rb

(3)

(4)

Direct Density-ratio Estimation with Dimensionality Reduction

8

where a penalty term λα⊤ α/2 is included for regularization purposes, and λ (≥ 0) is aregularization parameter that controls the strength of regularization. It is easy to conﬁrm
b can be analytically computed asthat the solution αbc + λI b )−1 h,b = (Hα(5)where I b is the b-dimensional identity matrix. Thanks to this analytic-form expression,
uLSIF is computationally eﬃcient compared with other density-ratio estimators which
involve non-linear optimization (Qin, 1998; Cheng & Chu, 2004; Huang et al., 2007;
Bickel et al., 2007; Sugiyama et al., 2008; Nguyen et al., 2010).
In the original uLSIF paper (Kanamori et al., 2009a), the above solution is further
modiﬁed asαbℓ ←− max(0, αbℓ ).
This modiﬁcation may improve the estimation accuracy in ﬁnite sample cases since thetrue density ratio is non-negative. Even so, we still use Eq.(5) as it is since it is diﬀerentiable with respect to U , where u = U x. This diﬀerentiability will play a crucial role in
the next section. Note that, even without the above round-up modiﬁcation, the solution
is guaranteed to converge to the optimal vector asymptotically both in parametric andnon-parametric cases (Kanamori et al., 2009a; Kanamori et al., 2009b). Thus omitting
the above modiﬁcation step may not have a strong eﬀect.
It was theoretically shown that uLSIF possesses superior theoretical properties instatistical convergence and numerical stability (Kanamori et al., 2009a; Kanamori et al.,2009b).
2.2.2

Basis Function Design

The performance of uLSIF depends on the choice of the basis functions {ψℓ (u)}bℓ=1 . As
explained below, the use of Gaussian basis functions would be reasonable:rb(u) =nnu∑

αℓ K(u, unu
ℓ ),

ℓ=1

where K(u, u′ ) is the Gaussian kernel with kernel width σ (> 0):()∥u − u′ ∥2′K(u, u ) = exp −.
2σ 2
By deﬁnition, the density ratio r(u) tends to take large values if pnu (u) is large andpde (u) is small; conversely, r(u) tends to be small (i.e., close to zero) if pnu (u) is small
and pde (u) is large. When a non-negative function is approximated by a Gaussian kernel
model, many kernels may be needed in the region where the output of the target function
is large; on the other hand, only a small number of kernels would be enough in the region
where the output of the target function is close to zero (see Figure 3). Following this9

OutputDirect Density-ratio Estimation with Dimensionality Reduction

Input

Figure 3: Heuristic of Gaussian kernel allocation.
heuristic, we allocate many kernels in the region where pnu (u) takes large values, whichnnu
may be approximately achieved by setting the Gaussian centers at {unui }i=1 .
nnu
Alternatively, we may locate (nnu + nde ) Gaussian kernels at both {unui }i=1 andnde{udej }j=1 . However, in our preliminary experiments, this did not further improve theperformance, but slightly increased the computational cost. When nnu is very large, justnnuusing all the test input points {unui }i=1 as Gaussian centers is already computationallynnurather demanding. To ease this problem, a subset of {unui }i=1 may be used as Gaussian
centers for computational eﬃciency, i.e., for a preﬁxed b (∈ {1, 2, . . . , nnu }), we use
rb(u) =b∑αℓ K(u, cℓ ),ℓ=1nnu
where {cℓ }bℓ=1 are template points randomly chosen from {unui }i=1 without replacement.
The performance of uLSIF depends on the kernel width σ and the regularization
parameter λ. Model selection of uLSIF is possible based on cross-validation (CV) withrespect to the error criterion (1) (Kanamori et al., 2009a).

3

Direct Density-ratio Estimation with Dimensionality Reduction

Although uLSIF was shown to be a useful density ratio estimation method (Kanamori
et al., 2009a), estimating the density ratio in high-dimensional spaces is still challenging.
In this section, we propose a new method of direct density-ratio estimation that involves
dimensionality reduction.

3.1

Hetero-distributional Subspace

Our basic idea is to ﬁrst ﬁnd a low-dimensional subspace in which the two densities are
signiﬁcantly diﬀerent from each other, and then perform density ratio estimation onlyin this subspace. Although a similar framework has been explored in Sugiyama et al.
(2010a), the current formulation is substantially more general than the previous approach,as explained below.

Direct Density-ratio Estimation with Dimensionality Reduction

10

Let u be an m-dimensional vector (1 ≤ m ≤ d) and v be a (d − m)-dimensional vector
deﬁned as
[ ]
[ ]
u
U
:=
x,
v
V
where U is an m × d matrix and V is a (d − m) × d matrix. In order to ensure theuniqueness of the decomposition, we assume (without loss of generality) 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 to each other (see Figure 4). Then the two
densities pnu (x) and pde (x) can be decomposed aspnu (x) = pnu (v|u)pnu (u),pde (x) = pde (v|u)pde (u).
The key theoretical assumption which forms the basis of our proposed algorithm isthat the conditional densities pnu (v|u) and pde (v|u) agree with each other, i.e., the two
densities pnu (x) and pde (x) are decomposed aspnu (x) = p(v|u)pnu (u),pde (x) = p(v|u)pde (u),where p(v|u) is the common conditional density. This assumption implies that themarginal densities of u are diﬀerent, but the conditional density of v given u is common to pnu (x) and pde (x). Then the density ratio is simpliﬁed asr(x) =

pnu (u)=: r(u).
pde (u)Thus, the density ratio does not have to be estimated in the entire d-dimensional space,but it is suﬃcient to estimate the ratio only in the m-dimensional subspace speciﬁed by
U.
Below, we will use the term, the hetero-distributional subspace, for indicating thesubspace speciﬁed by U in which pnu (u) and pde (u) are diﬀerent. More precisely, let S
be a subspace speciﬁed by U and V such thatS = {U ⊤ U x | pnu (v|u) = pde (v|u), u = U x, v = V x}.
Then the hetero-distributional subspace is deﬁned as the intersection of all subspaces
S. Intuitively, the hetero-distributional subspace is the ‘smallest’ subspace speciﬁed by
U such that pnu (v|u) and pde (v|u) agree with each other. We refer to the orthogonal
complement of the hetero-distributional subspace as the homo-distributional subspace (seeFigure 4).
This formulation is a generalization of the one proposed in Sugiyama et al. (2010a) inwhich the components in the hetero-distributional subspace and its complimentary subspace are assumed to be independent of each other. On the other hand, we do not impose

Direct Density-ratio Estimation with Dimensionality Reduction

11∈
Homo-distributional subspace

∈
m−d

v∈R

pnu v|u

pde v|u∈V ∈ R(m−d)×dp v|ux ∈ Rd∈U ∈ Rm×d∈u ∈ Rm|| subspace
Hetero-distributional
pnu u pde uFigure 4: Hetero-distributional subspace.
such an independence assumption in the current paper. As will be demonstrated in Section 4.1, this generalization has a remarkable eﬀect in extending the range of applicationsof direct density-ratio estimation with dimensionality reduction.
For the moment, we assume that the true dimensionality m of the hetero-distributional
subspace is known. Later, we explain how m is estimated from data.

3.2

Estimating Pearson Divergence Using uLSIF

Here, we introduce a criterion for hetero-distributional subspace search and how it isestimated from data.
We use the Pearson divergence (PD) as our criterion for evaluating the discrepancy
between two distributions. PD is a squared-loss variant of the Kullback-Leibler divergence(Kullback & Leibler, 1951), and is an instance of the f -divergences, which are also known
as the Csiszár f -divergences (Csiszár, 1967) or the Ali-Silvey distances (Ali & Silvey,1966). PD from pnu (x) to pde (x) is deﬁned and expressed as)2∫ (1
pnu (x)PD[pnu (x), pde (x)] :=− 1 pde (x)dx2
pde (x)∫11pnu (x)=pnu (x)dx − .
2pde (x)2
PD[pnu (x), pde (x)] vanishes if and only if pnu (x) = pde (x).
The following lemma (called the “data processing” inequality) characterizes the heterodistributional subspace in terms of PD.
Lemma 1 Let)2∫ (pnu (u)1− 1 pde (u)duPD[pnu (u), pde (u)] =2
pde (u)∫11pnu (u)=pnu (u)du − .
2pde (u)2(6)

Direct Density-ratio Estimation with Dimensionality Reduction

 

pnu x
p nu u
−
pde x
pde u

2

pde x x12

p nu u , p de u
pnu x , pde x(constant)Figure 5: Since PD[pnu (x), pde (x)] is constant, minimizing 12∫ ( pnu (x)(u)− ppnupde (x)de (u))2pde (x)dx

is equivalent to maximizing PD[pnu (u), pde (u)].
Then we have1
PD[pnu (x), pde (x)] − PD[pnu (u), pde (u)] =
2
≥ 0.

∫ (

pnu (x) pnu (u)−pde (x) pde (u))2pde (x)dx(7)

A proof of the above lemma (for a class of f -divergences) is provided in Appendix B. The right-hand side of Eq.(7) is non-negative, and it vanishes if and only
if pnu (v|u) = pde (v|u). Since PD[pnu (x), pde (x)] is a constant with respect to U , maximizing PD[pnu (u), pde (u)] with respect to U leads to pnu (v|u) = pde (v|u) (Figure 5).
That is, the hetero-distributional subspace can be characterized as the maximizer1 ofPD[pnu (u), pde (u)].
Although the hetero-distributional subspace can be characterized as the maximizer ofPD[pnu (u), pde (u)], we cannot directly ﬁnd the maximizer since pnu (u) and pde (u) are unknown. Here, we utilize a direct density-ratio estimator uLSIF (see Section 2.2) for approximating PD[pnu (u), pde (u)] from samples. Let us replace the density ratio pnu (u)/pde (u)in Eq.(6) by a density ratio estimator rb(u). Approximating the expectation over pnu (u)nnu
by an empirical average over {unui }i=1 , we have the following PD estimator.
c nu (u), pde (u)] :=
PD[p

nnu
1 ∑
1
rb(unui )− .
2nnu i=1
2

Since uLSIF was shown to be consistent (i.e., the solution converges to the optimalvalue) both in parametric and non-parametric cases (Kanamori et al., 2009a; Kanamori
c would be a consistent estimator of the true PD.
et al., 2009b), PD
1

As shown in Appendix B, the data processing inequality holds not only for PD, but also for any
f -divergences. Thus the characterization of the hetero-distributional subspace is not limited to PD, but
is applicable to all f -divergences.

Direct Density-ratio Estimation with Dimensionality Reduction

3.3

Least-squares(LHSS)

Hetero-distributional

Subspace

13

Searchc nu (u), pde (u)], our next task is to ﬁnd a maxGiven the uLSIF-based PD estimator PD[pc nu (u), pde (u)] with respect to U , and identify the hetero-distributional
imizer of PD[psubspace (cf. the data processing inequality given in Lemma 1). We call this procedure
Least-squares Hetero-distributional Subspace Search (LHSS).
We may employ various optimization techniques to ﬁnd a maximizer ofcPD[pnu (u), pde (u)]. Here we describe several possibilities.
3.3.1

Plain Gradient Algorithm

A gradient ascent algorithm would be a fundamental approach to non-linear smooth
optimization. We utilize the following lemma.
c nu (u), pde (u)] with respect to U is expressed as
Lemma 2 The gradient of PD[pbb∑cb ℓ,ℓ′∂ PD∂bhℓ 1 ∑
∂H
bℓ
bℓ α
α
α
b ℓ′
=
−
,
∂U
∂U
2 ℓ,ℓ′ =1
∂U
ℓ=1

(8)

b is given by Eq.(5) andwhere αnnu
∂ψℓ (unu
∂b
hℓ
1 ∑
i )
=
,
∂U
nnu i=1
∂U
)
(
nde
de
b ℓ,ℓ′′∂ψℓ (ude
∂ψ
)
(u
)
∂H
1 ∑
ℓ
j
j
de
=
ψℓ′ (ude,j ) + ψℓ (uj )∂Unde j=1∂U∂U∂ψℓ (u)1= − 2 (u − cℓ )(x − c′ℓ )⊤ ψℓ (u).
∂Uσ′dcℓ (∈ R ) is a pre-image of cℓ (∈ Rm ):

(9)
(10)
(11)

cℓ = U c′ℓ .
A proof of the above lemma is provided in Appendix C. Note that {bαℓ }bℓ=1 in Eq.(8)b in Eq.(5), which was taken into account when derivingb through Hc and hdepend on U
the gradient (see Appendix C). A plain gradient update rule is then given as
c
∂ PD
,
∂U
where t (> 0) is a learning rate. t may be chosen in practice by some approximate line
search method such as Armijo’s rule (Patriksson, 1999) or backtracking line search (Boyd& Vandenberghe, 2004).
A naive gradient update does not necessarily fulﬁll the orthonormality U U ⊤ = I m ,where I m is the m-dimensional identity matrix. Thus, after every gradient step, we
need to orthonormalize U by, e.g., the Gram-Schmidt process (Golub & Loan, 1996) toguarantee its orthonormality. However, this may be rather time-consuming.
U ←− U + t

Direct Density-ratio Estimation with Dimensionality Reduction

3.3.2

14

Natural Gradient AlgorithmcPD
In the Euclidean space, the ordinary gradient ∂∂Ugives the steepest direction. On theother hand, in the current setup, the matrix U is restricted to be a member of the Stiefel
manifold Sdm (R):Sdm (R) := {U ∈ Rm×d | U U ⊤ = I m }.
On a manifold, it is known that, not the ordinary gradient, but the natural gradientc(Amari, 1998) gives the steepest direction. The natural gradient ∇PD(U) at U is thec∂ PDdprojection of the ordinary gradient ∂U onto the tangent space of Sm (R) at U .
If the tangent space is equipped with the canonical metric, i.e., for any G and G′ in
the tangent space,)
1 (⟨G, G′ ⟩ = tr G⊤ G′ ,(12)2
the natural gradient is given by
(
)
c
c⊤
1
∂
PD
∂
PD
c
∇PD(U
)=
−U
U .
2 ∂U
∂U
c
Then the geodesic from U to the direction of the natural gradient ∇PD(U) over Sdm (R)can be expressed using t ∈ R as
)}
{ (
c
c⊤
∂
PD
∂
PD
−
U
,
U t := U exp t U ⊤∂U∂U
where ‘exp’ for a matrix denotes the matrix exponential, i.e., for a square matrix T ,exp(T ) :=

∞
∑
1

k!
k=0

T k.

(13)

Thus, line search along the geodesic in the natural gradient direction is equivalent toﬁnding a maximizer from
{U t | t ≥ 0}.
More details of geometric structure of the Stiefel manifold can be found in Nishimori and
Akaho (2005).
A natural gradient update rule is then given asU ←− U t ,where t (> 0) is the learning rate. Since the orthonormality of U is automatically satisﬁedin the natural gradient method, it would be computationally more eﬃcient than theplain gradient method. However, optimizing the m × d matrix U is still computationallyexpensive.

Direct Density-ratio Estimation with Dimensionality Reduction

15

Rotation acrossthe subspace

Rotation within
the subspace
Hetero-distributional
subspace

Figure 6: In the hetero-distributional subspace search, rotation which changes the subspace only matters (the solid arrow); rotation within the subspace (dotted arrow) can beignored 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).
3.3.3

Givens Rotation

Another simple strategy for optimizing U is to rotate the matrix in the plane spannedby two coordinate axes (which is called the Givens rotations; see Golub & Loan, 1996).
That is, we randomly choose a two-dimensional subspace spanned by the i-th and j-thvariables, and rotate the matrix U within this subspace:(i,j)

U ←− RθU,(i,j)

where Rθ is the rotation matrix by angle θ within the subspace spanned by the i-th
(i,j)
and j-th variables. Rθ is equal to the identity matrix except that its elements (i, i),(i, j), (j, i), and (j, j) form a two-dimensional rotation matrix:[ (i,j)] [](i,j)[Rθ ]i,i [Rθ ]i,jcos θ sin θ
=
.
(i,j)
(i,j)
− sin θ cos θ[Rθ ]j,i [Rθ ]j,j
The rotation angle θ (0 ≤ θ ≤ π) may be optimized by some secant method (Press et al.,1992).
As shown above, the update rule of the Givens rotations is computationally very
eﬃcient. However, since the update direction is not optimized as in the plain/natural
gradient methods, the Givens-rotation method could be potentially less eﬃcient as an
optimization strategy.
3.3.4

Subspace Rotation

Since we are searching for a subspace, rotation within the subspace does not have anyc (see Figure 6). This implies that the number ofinﬂuence on the objective value PD
parameters to be optimized in the gradient algorithm can be reduced.

Direct Density-ratio Estimation with Dimensionality Reduction

16

For a skew-symmetric matrix M (∈ Rd×d ), i.e., M ⊤ = −M , rotation of U can be
expressed as follows (Plumbley, 2005):
[ ]
U
[
]
I m O m,(d−m) exp(M ),Vwhere O d,d′ is the d × d′ matrix with all zeros, and exp(M ) is the matrix exponential of
M (see Eq.(13)). M = O d,d (i.e., exp(O d,d ) = I d ) corresponds to no rotation. Here we
update U through the matrix M .
Let us adopt Eq.(12) as the inner product in the space of skew-symmetric matrices.
Then we have the following lemma.
c with respect to M at M = O d,d is given by
Lemma 3 The derivative of PD
[
]
c
∂ PD
c
O m,m
V⊤
∂ PD
∂U
=
.
c
∂M
−( ∂ PD V ⊤ )⊤ O
M =O d,d(14)(d−m),(d−m)∂U

A proof of the above lemma is provided in Appendix D. The block structure of Eq.(14)has an intuitive explanation: the non-zero oﬀ-diagonal blocks correspond to the rotation
angles between the hetero-distributional subspace and its orthogonal complement whichc On the other hand, the derivative of rotation within
do aﬀect the objective function PD.
the two subspaces vanishes because this does not change the objective value. Thus thevariables to be optimized are only the angles corresponding to the non-zero oﬀ-diagonalcPD
blocks ∂∂UV ⊤ , which includes only m(d − m) variables. In contrast, the plain/natural
gradient algorithms optimize the matrix U , which contains md variables. Thus, when mis large, the subspace rotation approach may be computationally more eﬃcient than theplain/natural gradient algorithms.
The gradient ascent update rule of M is given byM ←− t

c
∂ PD
∂M

,
M =O d,dwhere t is a step-size. Then U is updated as

[ ]
U
[
]
U ←− I m O m,(d−m) exp(M ).
V

The conjugate gradient method (Golub & Loan, 1996) may be used for the update of M .
Following the update of U , its counterpart V also needs to be updated accordingly
since the hetero-distributional subspace and its complement speciﬁed by U and V should
be orthogonal to each other (see Figure 4). This can be achieved by setting[]⊤V ←− φ1 φ2 · · · φd−m ,where φ1 , φ2 , . . . , φd−m are orthonormal basis vectors in the orthogonal complement of
the hetero-distributional subspace.

Direct Density-ratio Estimation with Dimensionality Reduction

3.4

17

Proposed Algorithm: D3 -LHSS

Finally, we estimate the density ratio in the hetero-distributional subspace detected by
the above LHSS method.
A notable fact of the LHSS algorithm is that the density ratio estimator in the heterodistributional subspace has already been obtained during the hetero-distributional subspace search procedure. Thus, we do not need an additional estimation procedure—ourﬁnal solution is simply given by
rb(x) =

b
∑

b x),
α
bℓ ψℓ (Uℓ=1

b is a projection matrix obtained by the LHSS algorithm. {bwhere Uαℓ }bℓ=1 are the learnedb , which have been obtained and used when computing the gradient (seeparameters for U
Lemma 2).
This expression implies that if the dimensionality is not reduced (i.e., m = d), theproposed method agrees with the original uLSIF (see Section 2.2). Thus, the proposed
method could be regarded as a natural extension of uLSIF to high-dimensional data.
Given the true dimensionality m of the hetero-distributional subspace, we can estimate
the hetero-distributional subspace by the LHSS algorithm. When m is unknown, we may
choose the best dimensionality based on the CV score of the uLSIF estimator. We refer
to our proposed procedure D3 -LHSS (D-cube LHSS; Direct Density-ratio estimation with
Dimensionality reduction via Least-squares Hetero-distributional Subspace Search).
R
The complete procedure of D3 -LHSS is summarized in Figure 7. A MATLAB⃝
im3
plementation of D -LHSS is available from‘http://sugiyama-www.cs.titech.ac.jp/~sugi/software/D3LHSS/’.

4

Experiments

In this section, we investigate the experimental performance of the proposed method.
We employ the subspace rotation algorithm explained in Section 3.3.4 in our D3 -LHSS
implementation. In uLSIF, the number of parameters is ﬁxed to b = 100; the Gaussian
width σ and the regularization parameter λ are chosen based on cross-validation.

4.1

Illustrative Examples

First, we illustrate how the D3 -LHSS algorithm behaves.
As explained in Section 1, the previous D3 method, D3 -LFDA/uLSIF (Sugiyama et al.,2010a), has two potential weaknesses:• The component u inside the hetero-distributional subspace and its complementarycomponent v are assumed to be statistically independent (cf. Section 3.1).

Direct Density-ratio Estimation with Dimensionality Reduction

18

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 dimension m = 1, 2, . . . , dInitialize embedding matrix U m (∈ Rm×d );Repeat until U m converges
Choose Gaussian width σ and regularization parameter λ by CV;
Update U by some optimization 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 dimensionality mb that minimizes the CV score;
Set rb(x) = rbm
b (x);
Figure 7: Pseudo code of D3 -LHSS.
• Separability of samples drawn from two distributions implies that the two distributions are diﬀerent, but non-separability does not necessarily imply that the two
distributions are equivalent. Thus, D3 -LFDA/uLSIF may not be able to detect thesubspace in which the two distributions are diﬀerent, but samples are not reallyseparable.
Here, through numerical examples, we illustrate these weaknesses of D3 -LFDA/uLSIF,
and show these problems can be overcome by D3 -LHSS. Let us consider two-dimensional
examples (i.e., d = 2), and suppose that the two densities pnu (x) and pde (x) are diﬀerent
only in the one-dimensional subspace (i.e., m = 1) spanned by (1, 0)⊤ :x = (x(1) , x(2) )⊤ = (u, v)⊤ ,pnu (x) = p(v|u)pnu (u),pde (x) = p(v|u)pde (u).
Let nnu = nde = 1000. We use the following three datasets:“Rather-separate” dataset (Figure 8):p(v|u) = p(v) = N (v; 0, 12 ),pnu (u) = N (u; 0, 0.52 ),pde (u) = 0.5N (u; −1, 12 ) + 0.5N (u; 1, 12 ),where N (u; µ, σ 2 ) denotes the Gaussian density with mean µ and variance σ 2 withrespect to u. This is an easy and simple dataset for the purpose of illustrating theusefulness of dimensionality reduction in density ratio estimation.

Direct Density-ratio Estimation with Dimensionality Reduction

19“Highly-overlapped” dataset (Figure 9):p(v|u) = p(v) = N (v; 0, 12 ),pnu (u) = N (u; 0, 0.62 ),pde (u) = N (u; 0, 1.22 ).
Since v is independent of u, D3 -LFDA/uLSIF is still applicable in principle. However, unu and ude are highly overlapped and are not clearly separable. Thus thisdataset would be hard for D3 -LFDA/uLSIF.
“Dependent” dataset (Figure 10):p(v|u) = N (v; u, 12 ),pnu (u) = N (u; 0, 0.52 ),pde (u) = 0.5N (u; −1, 12 ) + 0.5N (u; 1, 12 ).
In this dataset, the conditional distribution p(v|u) is common, but the marginal
distributions pnu (v) and pde (v) are diﬀerent. Since v is not independent of u, thisdataset would be out of scope for D3 -LFDA/uLSIF.
The true hetero-distributional subspace for the “rather-separate” dataset is depictedby the dotted line in Figure 8(a); the solid line and the dashed line depict the heterodistributional subspace found by LHSS and LFDA with reduced dimensionality m = 1,
respectively. This graph shows that LHSS and LFDA both give very good estimates of
the true hetero-distributional subspace. In Figure 8(c), Figure 8(d), and Figure 8(e),density ratio functions estimated by the plain uLSIF without dimensionality reduction,
D3 -LFDA/uLSIF, and D3 -LHSS for the “rather-separate” dataset are depicted. Thesegraphs show that both D3 -LHSS and D3 -LFDA/uLSIF give much better estimates of thedensity ratio function (see Figure 8(b) for the proﬁle of the true density ratio function)
than the plain uLSIF without dimensionality reduction. Thus, the usefulness of dimensionality reduction in density ratio estimation was illustrated.
For the “highly-overlapped” dataset (Figure 9), LHSS gives a reasonable estimate of
the hetero-distributional subspace, while LFDA is highly erroneous due to less separability.
As a result, the density ratio function obtained by D3 -LFDA/uLSIF does not reﬂect thetrue redundant structure appropriately. On the other hand, D3 -LHSS still works well.
Finally, for the “dependent” dataset (Figure 10), LHSS gives an accurate estimate of
the hetero-distributional subspace. However, LFDA gives a highly biased solution since
the marginal distributions pnu (v) and pde (v) are no longer common in the “dependent”
dataset. Consequently, the density ratio function obtained by D3 -LFDA/uLSIF is highly
erroneous. In contrast, D3 -LHSS still works very well for the “dependent” dataset.
The experimental results for the “highly-overlapped” and “dependent” datasets illustrated typical failure modes of LFDA, and LHSS was shown to be able to successfully
overcome these weaknesses of LFDA.

Direct Density-ratio Estimation with Dimensionality Reduction

6

de

x
xnu
True Subspace
LHSS Subspace
LFDA Subspace

4
2
x(2)

20

0
−2
−4
−6
−6

−4

−2

0
x(1)

2

4

6

(a) Hetero-distributional subspace

3

3
2

5

2

5

1

1
0

0

0

0

−5

−5
0
5

−5

x

0

(2)

5

x (1)

x (1)

(b) r(x)(c) rb(x) by plain uLSIF.

−5

x

(2)

3

3
2

5

2

5

1

1
0

0

0

0
−5

−5
0
5

−5

x

(2)

0
5

x(1)

x

(d) rb(x) by D3 -LFDA/uLSIF.

(1)(e) rb(x) by D3 -LHSS.

Figure 8: “Rather-separate” dataset.

−5

x

(2)

Direct Density-ratio Estimation with Dimensionality Reduction

6

de

x
xnu
True Subspace
LHSS Subspace
LFDA Subspace

4
2
x(2)

21

0
−2
−4
−6
−6

−4

−2

0
x(1)

2

4

6

(a) Hetero-distributional subspace

2

2
5

1
0

0

5

1
0

0

−5

−5
0
5
x

−5

0

x(2)

5

(1)

x

−5

x(2)

(1)

(c) rb(x) by plain uLSIF.

(b) r(x)

2
1
5

0.5
0

0

5

1
0

0
−5

−5
0
5

−5

x(2)

0
5

x(1)

x(1)

(d) rb(x) by D3 -LFDA/uLSIF.

(e) rb(x) by D3 -LHSS.

Figure 9: “Highly-overlapped” dataset.

−5

x(2)

Direct Density-ratio Estimation with Dimensionality Reduction

6

de

x
xnu
True Subspace
LHSS Subspace
LFDA Subspace

4
2
x(2)

22

0
−2
−4
−6
−6

−4

−2

0
x(1)

2

4

6

(a) Hetero-distributional subspace

3

3
2

5

2

5

1

1
0

0

0

0

−5

−5
0
5
x

−5

0

(2)

x

5

(1)

x

−5

x(2)

(1)

(c) rb(x) by plain uLSIF.

(b) r(x)

2

3
5

1

2

5

1
0

0

0

0

−5

−5
0
5
x

−5

x

(2)

0
5

(1)

x

(d) rb(x) by D3 -LFDA/uLSIF.

(1)(e) rb(x) by D3 -LHSS.

Figure 10: “Dependent” dataset.

−5

x(2)

Direct Density-ratio Estimation with Dimensionality Reduction

4.2

23

Evaluation on Artificial Data

Next, we systematically compare the performance of the proposed D3 -LHSS with that ofthe plain uLSIF and D3 -LFDA/uLSIF for high-dimensional artiﬁcial data.
For the three datasets used in the previous experiments, we increase the entire dimensionality as d = 2, 3, . . . , 10 by adding dimensions consisting of standard normal noise.
The dimensionality of the hetero-distributional subspace is estimated based on the CVscore of uLSIF. We evaluate the error of a density ratio estimator rb(x) by∫
1
Error :=
(br(x) − r(x))2 pde (x)dx,(15)2
which uLSIF tries to minimize (see Section 2.2).
The left graphs in Figure 11 show the density-ratio estimation error averaged over
100 runs as functions of the entire input dimensionality d. The best method in termsof the mean error and comparable methods according to the t-test (Henkel, 1979) at thesigniﬁcance level 1% are speciﬁed by ‘◦’; otherwise methods are speciﬁed by ‘×’.
These plots show that, while the error of the plain uLSIF increases rapidly as the entire
dimensionality d increases, that of the proposed D3 -LHSS is kept moderate. Consequently,
the proposed method consistently outperforms the plain uLSIF. D3 -LHSS is comparable
to D3 -LFDA/uLSIF for the “rather-separate” dataset, and D3 -LHSS signiﬁcantly outperforms D3 -LFDA/uLSIF for the “highly-overlapped” and “dependent” datasets. Thus,
D3 -LHSS was overall shown to compare favorably with the other approaches.
The choice of the dimensionality of the hetero-distributional subspace in D3 -LHSS
and D3 -LFDA/uLSIF is illustrated in the middle and right columns of Figure 11; thedarker the color is, the more frequently the corresponding dimensionality is chosen. The
plots show that D3 -LHSS reasonably identiﬁes the true dimensionality (m = 1 in the current setup) for all the three datasets, while D3 -LFDA/uLSIF performs well only for the“rather-separate” dataset. This happened because D3 -LFDA/uLSIF cannot ﬁnd appropriate low-dimensional subspaces for the “highly-overlapped” and “dependent” datasets,
and therefore the CV scores misled the choice of subspace dimensionality.

4.3

Inlier-based Outlier Detection for Benchmark Data

Finally, we apply the proposed method to inlier-based outlier detection, i.e., ﬁnding outliers in an evaluation dataset based on another “model” dataset that only contains inliers(see Section A.2 for details).
We use the USPS hand-written digit dataset taken from the UCI Machine Learning
Repository (Asuncion & Newman, 2007). We regard samples in the class ‘1’ as inliers
and samples in other classes as outliers. We randomly take 500 samples from the class‘1’, and assign them to the model dataset. Then we randomly take 500 samples fromthe class ‘1’ without overlap, and 25 samples from one of the other classes. From thesesamples, density ratio estimation is performed and the outlier score is computed. Since
the USPS hand-written digit dataset contains 10 classes (i.e., from ‘0’ to ‘9’), we have 9

Direct Density-ratio Estimation with Dimensionality Reduction

24

100

90

0.4
0.3
0.2

9
80

8
70

7
60

6
50

5
40

4
30

3
20

Dimensionality chosen by CV

9

0.5

Error

10
90

Dimensionality chosen by CV

0.6

100

10

Plain uLSIF
D3−LFDA/uLSIF
D3−LHSS

2

0.1

80

8
70

7
60

6
50

5
40

4
30

3
20

2
10

10

1
0

1
0

2

3

4

5
6
7
8
Entire dimensionality d

9

2

10

3

4

5
6
7
8
Entire dimensionality d

9

0

10

2

3

4

5
6
7
8
Entire dimensionality d910(a) Density-ratio estimation error (b) Choice of Dimensionality by (c) Choice of Dimensionality by
D3 -LHSS
D3 -LFDA/uLSIF
100

0.25

10
90

90

Error

0.15

0.1

0.05

9
80

8
70

7
60

6
50

5
40

4
30

3
20

Dimensionality chosen by CV

9Dimensionality chosen by CV

0.2

100

10

Plain uLSIF
D3−LFDA/uLSIF
D3−LHSS

2

80

8
70

7
60

6
50

5
40

4
30

3
20

2
10

10

1
0

1
0

2

3

4

5
6
7
8
Entire dimensionality d

9

2

10

3

4

5
6
7
8
Entire dimensionality d

9

0

10

2

3

4

5
6
7
8
Entire dimensionality d

910(d) Density-ratio estimation error (e) Choice of Dimensionality by (f) Choice of Dimensionality by
D3 -LHSS
D3 -LFDA/uLSIF
100

90

0.3

0.2

0.1

9
80

8
70

7
60

6
50

5
40

4
30

3
20

2

Dimensionality chosen by CV

9

0.4
Error

10
90

Dimensionality chosen by CV

0.5

100

10

Plain uLSIF
3
D −LFDA/uLSIF
3
D −LHSS

80

8
70

7
60

6
50

5
40

4
30

3
20

2
10

1
0

10

1
0

2

3

4

5
6
7
8
Entire dimensionality d

9

10

2

3

4

5
6
7
8
Entire dimensionality d

9

10

0

2

3

4

5
6
7
8
Entire dimensionality d

910(g) Density-ratio estimation error (h) Choice of Dimensionality by (i) Choice of Dimensionality by
D3 -LHSS
D3 -LFDA/uLSIF

Figure 11: Top: “Rather-separate” dataset. Middle: “Highly-overlapped” dataset. Bottom: “Dependent” dataset. Left: Density-ratio estimation error (15) averaged over 100
runs as a function of the entire data dimensionality d. The best method in terms of themean 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 ‘×’. Center: The dimensionalityof the hetero-distributional subspace chosen by CV in LHSS. Right: The dimensionalityof the hetero-distributional subspace chosen by CV in LFDA.

Direct Density-ratio Estimation with Dimensionality Reduction25

diﬀerent tasks in total. The dimensionality of the samples is d = 256. For the D3 -LHSS
and D3 -LFDA/uLSIF methods, we choose the dimensionality of the hetero-distributional
subspace from m = 1, 2, . . . , 5 by cross-validation.
When evaluating the performance of outlier detection methods, it is important to take
into account both the detection rate (i.e., the amount of true outliers an outlier detection
algorithm can ﬁnd) and the detection accuracy (i.e., the amount of true inliers an outlier
detection algorithm misjudges as outliers). Since there is a trade-oﬀ between the detection
rate and the detection accuracy, we adopt the area under the ROC curve (AUC) as ourerror metric (Bradley, 1997).
The mean and standard deviation of AUC scores over 100 runs with diﬀerent random
seeds are summarized in Table 1, where the best method in terms of the mean AUC
and comparable methods according to the t-test at the signiﬁcance level 1% are speciﬁed by ‘◦ ’. The table shows that the proposed D3 -LHSS tends to outperform the plain
uLSIF and D3 -LFDA/uLSIF. It is also note worthy that D3 -LFDA/uLSIF is actuallyoutperformed by the plain uLSIF—the baseline method. This is perhaps because thenumerator and denominator datasets are highly overlapped in outlier detection scenarios,
so D3 -LFDA/uLSIF performs rather poorly (cf. Figure 9)
We also evaluate the performance of each method for an additional test dataset which
is not used for density ratio estimation. The test dataset consists of 100 randomly chosen
samples from the class ‘1’ and 5 randomly chosen samples from the outlier class (whichis the same as the evaluation dataset). The results are summarized in Table 2, showingthat the advantage of the proposed method is still valid in this more challenging scenario.

5

Conclusions

Density ratios are becoming quantities of interest in the machine learning and data mining
communities since it can be used for solving various important data processing tasks suchas non-stationarity adaptation, outlier detection, and feature selection (Sugiyama et al.,2009; Sugiyama et al., 2011). In this paper, we tackled a challenging problem of estimating
density ratios in high-dimensional spaces, and gave a new procedure in the framework of
Direct Density-ratio estimation with Dimensionality reduction (D3 ; D-cube). The basic
idea of D3 is to identify a subspace called the hetero-distributional subspace, in which two
distributions (corresponding to the numerator and denominator of the density ratio) are
diﬀerent.
In the existing approach of D3 (Sugiyama et al., 2010a), the hetero-distributional subspace is identiﬁed by ﬁnding a subspace in which samples drawn from the two distributions
are maximally separated from each other. To this end, supervised dimensionality reduction methods such as local Fisher discriminant analysis (LFDA) (Sugiyama, 2007) areutilized. This approach was shown to work well when the components inside and outside
the hetero-distributional subspace are statistically independent, and samples drawn fromthe two distributions are highly separable from each other in the hetero-distributional
subspace.

Direct Density-ratio Estimation with Dimensionality Reduction

26

Table 1: Outlier detection for the USPS hand-written digit dataset (d = 256). The
means (and standard deviations in the bracket) of AUC scores over 100 runs for theevaluation dataset are summarized. The best method in terms of the mean AUC value
and comparable methods according to the t-test at the signiﬁcance level 1% are speciﬁedby ‘◦ ’. The means (and standard deviations in the bracket) of the chosen dimensionality
by cross-validation are also included in the table.
Data
Digit 2
Digit 3
Digit 4
Digit 5
Digit 6
Digit 7
Digit 8
Digit 9
Digit 0
Average

D3 -LHSS
AUC
mb◦
0.956 (0.035) 4.3 (0.8)◦
0.967 (0.032) 4.4 (0.8)◦
0.907 (0.061) 4.4 (0.9)◦
0.965 (0.037) 4.3 (0.9)◦
0.974 (0.022) 4.4 (0.8)◦
0.924 (0.072) 4.4 (0.9)◦
0.929 (0.051) 4.2 (1.0)◦
0.942 (0.048) 4.6 (0.7)◦
0.986 (0.019) 4.2 (0.9)
0.950 (0.051) 4.4 (0.9)

D3 -LFDA/uLSIF
AUC
mb
0.889 (0.104) 1.7 (1.1)
0.868 (0.136) 1.8 (1.1)
0.825 (0.104) 1.4 (0.6)
0.882 (0.109) 1.6 (0.9)
0.891 (0.090) 1.7 (1.1)
0.642 (0.139) 2.3 (1.4)
0.804 (0.147) 1.8 (1.1)
0.790 (0.136) 1.8 (1.1)
0.920 (0.071) 1.9 (0.8)
0.835 (0.142) 1.8 (1.1)

Plain uLSIF
AUC
0.902 (0.038)
0.921 (0.039)
0.870 (0.036)
0.906 (0.037)
0.941 (0.029)
0.878 (0.035)
0.860 (0.033)
0.892 (0.035)◦
0.979 (0.019)
0.905 (0.049)

Table 2: Outlier detection for the USPS hand-written digit dataset (d = 256). The means(and standard deviations in the bracket) of AUC scores over 100 runs for unlearned test
dataset are summarized.
D3 -LHSS
D3 -LFDA/uLSIF
Plain uLSIF
Data
AUC
m
b
AUC
m
b
AUC
◦
0.946 (0.047) 4.3 (0.8) 0.817 (0.132) 1.7 (1.1) 0.905 (0.044)
Digit 2
Digit 3 ◦ 0.953 (0.061) 4.4 (0.8) 0.780 (0.161) 1.8 (1.1) 0.924 (0.045)
Digit 4 ◦ 0.880 (0.094) 4.4 (0.9) 0.767 (0.121) 1.4 (0.6) ◦ 0.870 (0.063)
Digit 5 ◦ 0.954 (0.057) 4.3 (0,9) 0.813 (0.142) 1.6 (0.9) 0.906 (0.047)
Digit 6 ◦ 0.959 (0.052) 4.4 (0.8) 0.806 (0.141) 1.7 (1.1) 0.939 (0.040)
Digit 7 ◦ 0.909 (0.079) 4.4 (0.9) 0.689 (0.173) 2.3 (1.4) 0.877 (0.056)
Digit 8 ◦ 0.903 (0.078) 4.2 (1.0) 0.741 (0.173) 1.8 (1.1) 0.861 (0.049)
Digit 9 ◦ 0.932 (0.072) 4.6 (0.7) 0.793 (0.128) 1.8 (1.1) 0.894 (0.054)
Digit 0 ◦ 0.982 (0.039) 4.2 (0.9) 0.859 (0.098) 1.9 (0.8) ◦ 0.982 (0.022)
Average
0.935 (0.073) 4.4 (0.9) 0.785 (0.150) 1.8 (1.1) 0.906 (0.060)

Direct Density-ratio Estimation with Dimensionality Reduction

27

However, as illustrated in Section 4.1, violation of these conditions can cause significant performance degradation. This problem can be overcome in principle by ﬁnding asubspace such that two conditional distributions are similar to each other in its complementary subspace. However, comparing conditional distributions is a cumbersome task.
To cope with this problem, we ﬁrst proved that the hetero-distributional subspace canbe characterized as the subspace in which two marginal distributions are maximally different under the Pearson divergence (Lemma 1). Based on this lemma, we proposed a
new algorithm for ﬁnding the hetero-distributional subspace called Least-squares Heterodistributional Subspace Search (LHSS). Since a density-ratio estimation method is utilized during hetero-distributional subspace search in the LHSS procedure, an additional
density-ratio estimation step is not needed after hetero-distributional subspace search.
Thus, two steps in the previous method (hetero-distributional subspace search followed
by density ratio estimation in the identiﬁed subspace) were merged into a single step (seeFigure 2). The proposed single-shot procedure, D3 -LHSS (D-cube LHSS), was shown tobe able to overcome the limitations of the D3 -LFDA/uLSIF approach through experiments.
In the experiments in Section 4, we employed the subspace rotation algorithm explained in Section 3.3.4 in our D3 -LHSS implementation. Although we experimentallyfound that the subspace rotation algorithm is useful, this does not necessarily mean thatsubspace rotation is always the best performing algorithm. Other approaches explained in
Section 3.3 may also be useful in some situations. Further investigating the optimization
issue is an important future work.
We gave a general proof of the data processing inequality (Lemma 1) for a class of f divergences (Ali & Silvey, 1966; Csiszár, 1967). Thus, the hetero-distributional subspaceis characterized not only by the Pearson divergence, but also by any f -divergences. Sincea framework of density ratio estimation for f -divergences has been provided in Nguyen
et al. (2010), an interesting future direction is to develop hetero-distributional subspace
search methods for general f -divergences.

Acknowledgments
MS was supported by SCAT, AOARD, and the JST PRESTO program. MY was supported by the JST PRESTO program. We thank Satoshi Hara for having performed
preliminary experiments using an earlier version of the proposed method. Our special
thanks also go to anonymous reviewers for their comments.

References
Akaike, H. (1974). A new look at the statistical model identiﬁcation. IEEE Transactions
on Automatic Control, AC-19, 716–723.
Akiyama, T., Hachiya, H., & Sugiyama, M. (2010). Eﬃcient exploration through active

Direct Density-ratio Estimation with Dimensionality Reduction

28

learning for value function approximation in reinforcement learning. Neural Networks,23, 639–648.
Ali, S. M., & Silvey, S. D. (1966). A general class of coeﬃcients of divergence of one
distribution from another. Journal of the Royal Statistical Society, Series B, 28, 131–
142.
Amari, S. (1998). Natural gradient works eﬃciently in learning. Neural Computation, 10,
251–276.
Aronszajn, N. (1950). Theory of reproducing kernels. Transactions of the American
Mathematical Society, 68, 337–404.
Asuncion, A., & Newman, D. (2007). UCI machine learning repository.
Bensaid, N., & Fabre, J. P. (2007). Optimal asymptotic quadratic error of kernel estimators of Radon-Nikodym derivatives for strong mixing data. Journal of Nonparametric
Statistics, 19, 77–88.
Bickel, S., Bogojeska, J., Lengauer, T., & Scheﬀer, T. (2008). Multi-task learning for HIV
therapy screening. Proceedings of 25th Annual International Conference on Machine
Learning (ICML2008) (pp. 56–63).
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).
Bishop, C. M. (2006). Pattern recognition and machine learning. New York, NY, USA:
Springer.
Boyd, S., & Vandenberghe, L. (2004). Convex optimization. Cambridge, UK: Cambridge
University Press.
Bradley, A. P. (1997). The use of the area under the ROC curve in the evaluation ofmachine learning algorithms. Pattern Recognition, 30, 1145–1159.
Chen, S.-M., Hsu, Y.-S., & Liaw, J.-T. (2009). On kernel estimators of density ratio.
Statistics, 43, 463–479.
Cheng, K. F., & Chu, C. K. (2004). Semiparametric density estimation under a twosample density ratio model. Bernoulli, 10, 583–604.
Cover, T. M., & Thomas, J. A. (2006). Elements of information theory. Hoboken, NJ,
USA: John Wiley & Sons, Inc. 2nd edition.
Csiszár, I. (1967). Information-type measures of diﬀerence of probability distributions
and indirect observation. Studia Scientiarum Mathematicarum Hungarica, 2, 229–318.

Direct Density-ratio Estimation with Dimensionality Reduction

29

Ćwik, J., & Mielniczuk, J. (1989). Estimating density ratio with application to discriminant analysis. Communications in Statistics: Theory and Methods, 18, 3057–3069.
Fishman, G. S. (1996). Monte Carlo: Concepts, algorithms, and applications. Berlin,
Germany: Springer-Verlag.
Gijbels, I., & Mielniczuk, J. (1995). Asymptotic properties of kernel estimators of theRadon-Nikodym derivative with applications to discriminant analysis. Statistica Sinica,
5, 261–278.
Golub, G. H., & Loan, C. F. V. (1996). Matrix computations. Baltimore, MD, USA:
Johns Hopkins University Press.
Hachiya, H., Akiyama, T., Sugiyama, M., & Peters, J. (2009a). Adaptive importancesampling for value function approximation in oﬀ-policy reinforcement learning. Neural
Networks, 22, 1399–1410.
Hachiya, H., Peters, J., & Sugiyama, M. (2009b). Eﬃcient sample reuse in EM-based
policy search. Machine Learning and Knowledge Discovery in Databases (pp. 469–484).
Berlin: Springer.
Härdle, W., Müller, M., Sperlich, S., & Werwatz, A. (2004). Nonparametric and semiparametric models. Berlin, Germany: Springer.
Henkel, R. E. (1979). Tests of significance. Beverly Hills, CA, USA.: SAGE Publication.
Hido, S., Tsuboi, Y., Kashima, H., Sugiyama, M., & Kanamori, T. (2008). Inlier-based
outlier detection via direct density ratio estimation. Proceedings of IEEE International
Conference on Data Mining (ICDM2008) (pp. 223–232). Pisa, Italy.
Hido, S., Tsuboi, Y., Kashima, H., Sugiyama, M., & Kanamori, T. (2010). Statisticaloutlier detection using direct density ratio estimation. Knowledge and Information
Systems. to appear.
Huang, J., Smola, A., Gretton, A., Borgwardt, K. M., & Schölkopf, B. (2007). Correcting
sample selection bias by unlabeled data. In B. Schölkopf, J. Platt and T. Hoﬀman(Eds.), Advances in neural information processing systems 19, 601–608. Cambridge,
MA, USA: MIT Press.
Hulle, M. M. V. (2005). Edgeworth approximation of multivariate diﬀerential entropy.
Neural Computation, 17, 1903–1910.
Jacob, P., & Oliveira, P. E. (1997). Kernel estimators of general Radon-Nikodym derivatives. Statistics, 30, 25–46.
Kanamori, T., Hido, S., & Sugiyama, M. (2009a). A least-squares approach to direct
importance estimation. Journal of Machine Learning Research, 10, 1391–1445.

Direct Density-ratio Estimation with Dimensionality Reduction

30

Kanamori, T., & Shimodaira, H. (2003). Active learning algorithm using the maximum
weighted log-likelihood estimator. Journal of Statistical Planning and Inference, 116,
149–162.
Kanamori, T., Suzuki, T., & Sugiyama, M. (2009b). Condition number analysis of kernelbased density ratio estimation (Technical Report). arXiv.
Kawahara, Y., & Sugiyama, M. (2009). Change-point detection in time-series data bydirect density-ratio estimation. Proceedings of 2009 SIAM International Conference on
Data Mining (SDM2009) (pp. 389–400). Sparks, Nevada, USA.
Kraskov, A., Stögbauer, H., & Grassberger, P. (2004). Estimating mutual information.
Physical Review E, 69, 066138.
Kullback, S., & Leibler, R. A. (1951). On information and suﬃciency. Annals of Mathematical Statistics, 22, 79–86.
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. to appear.
Nishimori, Y., & Akaho, S. (2005). Learning algorithms utilizing quasi-geodesic ﬂows on
the Stiefel manifold. Neurocomputing, 67, 106–135.
Patriksson, M. (1999). Nonlinear programming and variational inequality problems. Dordrecht, the Netherlands: Kluwer Academic.
Petersen, K. B., & Pedersen, M. S. (2008). The matrix cookbook (Technical Report).
Technical University of Denmark.
Plumbley, M. D. (2005). Geometrical methods for non-negative ICA: Manifolds, Liegroups and toral subalgebras. Neurocomputing, 67, 161–197.
Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. (1992). Numerical
recipes in C. Cambridge, UK: Cambridge University Press. 2nd edition.
Qin, J. (1998). Inferences for case-control and semiparametric two-sample density ratio
models. Biometrika, 85, 619–639.
Quiñonero-Candela, J., Sugiyama, M., Schwaighofer, A., & Lawrence, N. (Eds.). (2009).
Dataset shift in machine learning. Cambridge, Massachusetts, USA: MIT Press.
Shimodaira, H. (2000). Improving predictive inference under covariate shift by weighting
the log-likelihood function. Journal of Statistical Planning and Inference, 90, 227–244.
Silverman, B. W. (1978). Density ratios, empirical likelihood and cot death. Journal of
the Royal Statistical Society, Series C, 27, 26–33.

Direct Density-ratio Estimation with Dimensionality Reduction

31

Silverman, B. W. (1986). Density estimation for statistics and data analysis. London,
UK: Chapman and Hall.
Smola, A., Song, L., & Teo, C. H. (2009). Relative novelty detection. Proceedings of Twelfth International Conference on Artificial Intelligence and Statistics (AISTATS2009) (pp. 536–543). Clearwater Beach, FL, USA.
Steinwart, I. (2001). On the inﬂuence of the kernel on the consistency of support vector
machines. Journal of Machine Learning Research, 2, 67–93.
Stone, M. (1974). Cross-validatory choice and assessment of statistical predictions. Journal of the Royal Statistical Society, Series B, 36, 111–147.
Storkey, A., & Sugiyama, M. (2007). Mixture regression for covariate shift. Advances in
Neural Information Processing Systems 19 (pp. 1337–1344). Cambridge, Massachusetts,
USA: MIT Press.
Sugiyama, M. (2006). Active learning in approximately linear regression based on conditional expectation of generalization error. Journal of Machine Learning Research, 7,
141–166.
Sugiyama, M. (2007). Dimensionality reduction of multimodal labeled data by local Fisher
discriminant analysis. Journal of Machine Learning Research, 8, 1027–1061.
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., Kanamori, T., Suzuki, T., Hido, S., Sese, J., Takeuchi, I., & Wang, L.
(2009). A density-ratio framework for statistical data processing. IPSJ Transactions
on Computer Vision and Applications, 1, 183–208.
Sugiyama, M., & Kawanabe, M. (2010). Covariate shift adaptation: Towards machinelearning in non-stationary environment. Cambridge, Massachusetts, USA: MIT Press.
to appear.
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., & Müller, K.-R. (2005). Input-dependent estimation of generalization error
under covariate shift. Statistics & Decisions, 23, 249–279.
Sugiyama, M., & Nakajima, S. (2009). Pool-based active learning in approximate linear
regression. Machine Learning, 75, 249–274.

Direct Density-ratio Estimation with Dimensionality Reduction

32

Sugiyama, M., Suzuki, T., & Kanamori, T. (2010b). Density ratio estimation: A comprehensive review. RIMS Kokyuroku, 10–31.
Sugiyama, M., Suzuki, T., & Kanamori, T. (2011). Density ratio estimation in machinelearning: A versatile tool for statistical data processing. Cambridge, UK: Cambridge
University Press.
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.
(2010c). Least-squares conditional density estimation. IEICE Transactions on Information and Systems, E93-D, 583–594.
Sun, J., & Woodroofe, M. (1997). Semi-parametric estimates under biased sampling.
Statistica Sinica, 7, 545–575.
Suzuki, T., & Sugiyama, M. (2009). Estimating squared-loss mutual information for independent component analysis. Independent Component Analysis and Signal Separation(pp. 130–137). Berlin, Germany: Springer.
Suzuki, T., & Sugiyama, M. (2010). Suﬃcient dimension reduction via squared-loss mutual information estimation. Proceedings of the Thirteenth International Conference on
Artificial Intelligence and Statistics (AISTATS2010) (pp. 804–811). Sardinia, Italy.
Suzuki, T., Sugiyama, M., Kanamori, T., & Sese, J. (2009a). Mutual information estimation reveals global associations between stimuli and biological processes. BMC
Bioinformatics, 10, S52.
Suzuki, T., Sugiyama, M., Sese, J., & Kanamori, T. (2008). Approximating mutual
information by maximum likelihood density ratio estimation. Proceedings of ECMLPKDD2008 Workshop on New Challenges for Feature Selection in Data Mining andKnowledge Discovery 2008 (FSDM2008) (pp. 5–20). Antwerp, Belgium.
Suzuki, T., Sugiyama, M., & Tanaka, T. (2009b). Mutual information approximation via
maximum likelihood estimation of density ratio. Proceedings of 2009 IEEE International
Symposium on Information Theory (ISIT2009) (pp. 463–467). Seoul, Korea.
Takeuchi, I., Nomura, K., & Kanamori, T. (2009). Nonparametric conditional density
estimation using piecewise-linear solution path of kernel quantile regression. NeuralComputation, 21, 533–559.
Tsuboi, Y., Kashima, H., Hido, S., Bickel, S., & Sugiyama, M. (2009). Direct density ratio
estimation for large-scale covariate shift adaptation. Journal of Information Processing,17, 138–155.

Direct Density-ratio Estimation with Dimensionality Reduction

33

Ueki, K., Sugiyama, M., & Ihara, Y. (2010). Perceived age estimation under lighting condition change by covariate shift adaptation. 20th International Conference on Pattern
Recognition (ICPR2010) (pp. 3400–3403). Istanbul, Turkey.
Vapnik, V. N. (1998). Statistical learning theory. New York, NY, USA: Wiley.
Wahba, G. (1990). Spline models for observational data. Philadelphia, PA, USA: Society
for Industrial and Applied Mathematics.
Wiens, D. P. (2000). Robust weights and designs for biased regression models: Leastsquares and generalized M-estimation. Journal of Statistical Planning and Inference,
83, 395–412.
Yamada, M., & Sugiyama, M. (2010). Dependence minimizing regression with model
selection for non-linear causal inference under non-Gaussian noise. Proceedings of theTwenty-Fourth AAAI Conference on Artificial Intelligence (AAAI2010) (pp. 643–648).
Atlanta, Georgia, USA: The AAAI Press.
Yamada, M., Sugiyama, M., & Matsui, T. (2010). Semi-supervised speaker identiﬁcation
under covariate shift. Signal Processing, 90, 2353–2361.
Zadrozny, B. (2004).
Learning and evaluating classiﬁers under sample selection
bias. Proceedings of the Twenty-First International Conference on Machine Learning(ICML2004) (pp. 903–910). New York, NY, USA: ACM Press.

A

Usage of Density Ratios in Data Processing

We are interested in estimating density ratios since they are useful in various data processing tasks. Here, we brieﬂy review possible usage of density ratios (Sugiyama et al.,2009; Sugiyama et al., 2011).

A.1

Covariate Shift Adaptation

Covariate shift (Shimodaira, 2000) is a situation in supervised learning where input distributions change between the training and test phases, but the conditional distribution of
outputs given inputs remains unchanged. Extrapolation (i.e., prediction is made outside
the training region) would be a typical example of covariate shift. Standard learning techniques such as maximum likelihood estimation are biased under covariate shift; the bias
caused by covariate shift can be asymptotically canceled by weighting the loss function according to the importance 2 (Shimodaira, 2000; Zadrozny, 2004; Sugiyama & Müller, 2005;
Sugiyama et al., 2007; Quiñonero-Candela et al., 2009; Sugiyama & Kawanabe, 2010).

2

The test input density over the training input density is referred to as the importance in the context
of importance sampling (Fishman, 1996).

Direct Density-ratio Estimation with Dimensionality Reduction

34

The basic idea of covariate shift adaptation is summarized in the following importancesampling identity:∫
E [loss(x)] = loss(x)pnu (x)dx
pnu (x)∫= loss(x)r(x)pde (x)dx = E [loss(x)r(x)].
pde (x)That is, the expectation of a function loss(x) over pnu (x) can be computed by theimportance-weighted expectation of loss(x) over pde (x). Similarly, standard model selection criteria such as cross-validation (Stone, 1974; Wahba, 1990) or Akaike’s information
criterion (Akaike, 1974) lose their unbiasedness under covariate shift; proper unbiasedness
can be recovered by modifying the methods based on importance weighting (Shimodaira,
2000; Zadrozny, 2004; Sugiyama & Müller, 2005; Sugiyama et al., 2007). Furthermore,
the performance of active learning or the experiment design, i.e., the training input distribution is designed by the user to enhance the generalization performance, could also
be improved by the use of the importance (Wiens, 2000; Kanamori & Shimodaira, 2003;
Sugiyama, 2006; Sugiyama & Nakajima, 2009).
Thus, the importance plays a central role in covariate shift adaptation, and densityratio estimation methods could be used for reducing the estimation bias under covariate
shift. Examples of successful real-world applications include brain-computer interface(Sugiyama et al., 2007), robot control (Hachiya et al., 2009a; Akiyama et al., 2010;
Hachiya et al., 2009b), speaker identiﬁcation (Yamada et al., 2010), age prediction fromface images (Ueki et al., 2010), wafer alignment in semiconductor exposure apparatus(Sugiyama & Nakajima, 2009), and natural language processing (Tsuboi et al., 2009). A
similar importance-weighting idea also plays a central role in domain adaptation (Storkey& Sugiyama, 2007) and multi-task learning (Bickel et al., 2008).

A.2

Inlier-based Outlier Detection

Let us consider an outlier detection problem of ﬁnding irregular samples in a dataset(“evaluation dataset”) based on another dataset (“model dataset”) that only contains
regular samples. Deﬁning the density ratio over the two sets of samples, we can seethat the density ratio values for regular samples 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 (Hido et al., 2008; Smola et al., 2009; Hido
et al., 2010). Since the evaluation dataset usually has a wider support than the model
dataset, we regard the evaluation dataset as samples corresponding to pde (x) and themodel dataset as samples corresponding to pnu (x). Then outliers tend to have smaller
density-ratio values (i.e., close to zero). As such, density-ratio estimation methods could
be employed in outlier detection scenarios.
A similar idea could be used for change-point detection in time-series (Kawahara &Sugiyama, 2009).

Direct Density-ratio Estimation with Dimensionality Reduction

A.3

35

Conditional Density 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 density q(y|x). When thedomain of x is continuous, conditional density estimation is not straightforward since a
naive empirical approximation cannot be used (Bishop, 2006; Takeuchi et al., 2009).
In the context of density ratio estimation, 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 byr(x, y) :=q(x, y)= q(y|x),q(x)where q(x) is the marginal density of x. Then a density-ratio estimation method directly
gives an estimate of the conditional density (Sugiyama et al., 2010c).
When y is categorical, the same method can be used for probabilistic classiﬁcation(Sugiyama, 2010).

A.4

Mutual Information 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). Let us denote the marginal densities of x and y by q(x) and q(y),respectively. Then mutual information MI(X, Y ) between random variables X and Y isdeﬁned by∫∫q(x, y)MI(X, Y ) :=q(x, y) log
dxdy,
q(x)q(y)which plays a central role in information theory (Cover & Thomas, 2006).
Let us regard {(xk , y k )}nk=1 as samples corresponding to the numerator of the density
ratio and {(xk , y k′ )}nk,k′ =1 as samples corresponding to the denominator of the density
ratio, i.e.,r(x, y) :=q(x, y).
q(x)q(y)Then mutual information can be directly estimated using a density-ratio estimation
method (Suzuki et al., 2008; Suzuki et al., 2009b). General divergence functionals canalso be estimated in a similar way (Nguyen et al., 2010).
Mutual information can be used for measuring independence between random variables (Kraskov et al., 2004; Hulle, 2005) since it vanishes if and only if X and Y arestatistically independent. Thus density-ratio estimation methods are applicable, e.g.,
to variable selection (Suzuki et al., 2009a), independent component analysis (Suzuki &Sugiyama, 2009), supervised dimensionality reduction (Suzuki & Sugiyama, 2010), andcausal inference (Yamada & Sugiyama, 2010).

Direct Density-ratio Estimation with Dimensionality Reduction

B

36

Proof of Lemma 1

Here, let us consider the f -divergences (Ali & Silvey, 1966; Csiszár, 1967) and prove asimilar inequality for a broader class of divergences. An f -divergence is deﬁned using a
convex function f such that f (1) = 0 as()∫
pnu (x)If [pnu (x), pde (x)] := pde (x)fdx.
pde (x)
The f -divergence is reduced to the Kullback-Leibler divergence if
f (t) = − log t,
and the Pearson divergence if
1f (t) = (t − 1)2 .
2
Using Jensen’s inequality (Bishop, 2006), we have()∫∫
pnu (v|u)pnu (u)If [pnu (x), pde (x)] =pde (v|u)pde (u)fdudv
pde (v|u)pde (u))(∫∫
pnu (v|u)pnu (u)pde (v|u)
dv du≥ pde (u)fpde (v|u)pde (u))(∫∫
pnu (u)pnu (v|u)dv du= pde (u)fpde (u)()∫
pnu (u)= pde (u)fdupde (u)= If [pnu (u), pde (u)].
Thus, we have
If [pnu (x), pde (x)] − If [pnu (u), pde (u)] ≥ 0,
and the equality holds if and only if pnu (v|u) = pde (v|u).

C

Proof of Lemma 2

For
c + λI b )−1 ,F = (Hc nu (u), pde (u)] can be expressed as
PD[p
b
1
1∑ b
c
PD[pnu (u), pde (u)] =αb ℓ hℓ −
2 ℓ=1
2

1
1 ∑bb
=
hℓ hℓ′ Fℓ,ℓ′ − .
2 ℓ,ℓ′ =12bDirect Density-ratio Estimation with Dimensionality Reduction

37

Thus, its partial derivative with respect to U is given by
b
b
∑
c
∂ PD
∂b
hℓ 1 ∑ b b ∂Fℓ,ℓ′=
αbℓ+hℓ hℓ′
.
∂U
∂U
2
∂U
′
ℓ=1
ℓ,ℓ =1

(16)

Since
∂B −1
∂B −1
= −B −1B∂t∂tholds for a square invertible matrix B (Petersen & Pedersen, 2008), it holds thatc∂Fc + λI b )−1 ∂ H (Hc + λI b )−1 .
= −(H∂Uk,k′∂Uk,k′
Then we haveb∑ℓ,ℓ′ =1

b
hℓ b
hℓ′

[

∂F
∂Uk,k′]cbb ⊤ (Hc + λI b )−1 hc + λI b )−1 ∂ H (H
= −h
′
∂U
′
k,k
ℓ,ℓ
[
]
b
∑
c
∂H
α
bℓ α
b ℓ′
=−
.
∂Uk,k′′ℓ,ℓ′ =1ℓ,ℓ

Substituting this into Eq.(16), we obtain Eq.(8). Eqs.(9) and (10) are clear from Eqs.(3)and (2). Finally, we prove Eq.(11). The basis function ψℓ (u) can be expressed as()∥U (x − c′ℓ )∥2ψℓ (u) = ψℓ (U x) = exp −
.
2σ 2
⊤

⊤

A Aa
= 2Aa⊤ a (Petersen & Pedersen, 2008), we haveSince ∂a ∂A()∂ψℓ (u)1∥U (x − c′ℓ )∥2′′ ⊤= − 2 U (x − cℓ )(x − cℓ ) exp −
,
∂U
σ
2σ 2

from which we obtain Eq.(11).

D

Proof of Lemma 3

The proof we provide here essentially follows the argument in Plumbley (2005).
For
[ ]
[ ]
U
U0
,
W =
, W0 =V
V0

Direct Density-ratio Estimation with Dimensionality Reduction

38

rotation W from some W 0 can be expressed as follows (Plumbley, 2005):
W = exp(M )W 0 ,(17)where M is some skew-symmetric matrix. Let us consider the space of skew-symmetricmatrices, and let E be an element in that space with unit length. Then the gradient ofc∂ PDc, in this space is given as an element whose
a function PD(M)
withrespect
to
M
,
∂M
⟨
⟩
c
c
inner product ∂ PD , E is equal to the derivative of PD(M) in the direction E (Plumbley,∂M2005). Thus, for M = tE with t being a scalar, we have⟨⟩cc
∂ PD(tE)∂ PD(M
)
=
,E .
∂t
∂M
If we adopt Eq.(12) as the inner product of the space of skew-symmetric matrices, we
have
(
)
c
c
∂ PD(tE)
1∂ PD(M
) ⊤
= tr
E
.
(18)
∂t
2
∂M
can be expressed as follows
On the other hand, from Eq.(17) with M = tE, ∂W∂t(Petersen & Pedersen, 2008):∂W= E exp(tE)W 0 = EW .
∂tcThen ∂ PD(tE)
can be expressed as
∂t
(
)
(
)
⊤
c
c
c
∂ PD(tE)∂ PD ∂W∂ PD ⊤ ⊤
= tr
= tr
W E
.
∂t
∂W ∂t
∂W

(19)

Since E is skew-symmetric, it can be expressed as111
1
E = E + E = E − E⊤.
2222
Substituting this into Eq.(19), we have
(
)
(
)
c
c
c
1
∂ PD(tE)∂ PD ⊤ ⊤1∂ PD ⊤
= tr
W E
W E
− tr
∂t
2
∂W
2
∂W
(
(
)
)
c
c⊤
∂ PD
∂
PD
1
1
W ⊤ E ⊤ − tr W
E⊤
= tr
2
∂W
2
∂W
((
)
)
c
c⊤
∂ PD
∂
PD
1
W⊤ − W
E⊤ .
= tr
2
∂W
∂W

(20)

Direct Density-ratio Estimation with Dimensionality Reduction

39

Combining Eqs.(18) and (20), we have
c
c⊤
c
∂ PD
∂ PD
∂ PD
⊤
=
W −W
∂M
∂W
∂W
(
)
c
c ⊤
⊤
∂ PD
∂ PD
 ∂U U − U ∂U
=
(
)
c ⊤
c
⊤
∂ PD
∂ PD
U
−
V
∂V
∂U

c
∂ PD
V⊤−U
∂U
c
∂ PD
V⊤−V
∂V

(
(

c
∂ PD
∂V
c
∂ PD
∂V

)⊤ 

)⊤  .

(21)

b

hℓ
ℓ (u)
Eq.(11) implies that ∂ψ∂U
U ⊤ is symmetric. Then Eqs.(8) and (9) imply that ∂∂U
U⊤

and

b ′
∂H
c
ℓ,ℓ
PD
U ⊤ are also symmetric. Consequently, Eq.(10) imply that ∂∂U
U ⊤ is symmetric:
∂U

c
∂ PD
U⊤ =
∂U

(

c
∂ PD
U⊤
∂U

)⊤

⊤

c
∂ PD
.
=U
∂U

c
Since the range of V is assumed to be orthogonal to the range of U (see Section 3.1), PD
is independent of V , and thus we havec∂ PD
= O (d−m),d ,∂Vwhere O d,d′ is the d × d′ matrix with all zeros. Then Eq.(21) yields
[
]
c
∂ PD
c
V⊤
O m,m
∂ PD
∂U
=
,
c
PD
∂M
−( ∂∂UV ⊤ )⊤ O (d−m),(d−m)
which concludes the proof.