Journal of Machine Learning Research 10 (2009) 1391-1445

Submitted 3/08; Revised 4/09; Published 7/09

A Least-squares Approach to Direct Importance Estimation∗
Takafumi Kanamori

KANAMORI @ IS . NAGOYA - U . AC . JP

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

Shohei Hido

HIDO @ JP. IBM . COM

IBM Research
Tokyo Research Laboratory
1623-14 Shimotsuruma, Yamato-shi, Kanagawa 242-8502, Japan

Masashi Sugiyama

SUGI @ CS . TITECH . AC . JP

Department of Computer Science
Tokyo Institute of Technology
2-12-1 O-okayama, Meguro-ku, Tokyo 152-8552, Japan

Editor: Bianca Zadrozny

Abstract
We address the problem of estimating the ratio of two probability density functions, which is often
referred to as the importance. The importance values can be used for various succeeding tasks
such as covariate shift adaptation or outlier detection. In this paper, we propose a new importanceestimation method that has a closed-form solution; the leave-one-out cross-validation score can also
be computed analytically. Therefore, the proposed method is computationally highly efficient andsimple to implement. We also elucidate theoretical properties of the proposed method such as theconvergence rate and approximation error bounds. Numerical experiments show that the proposed
method is comparable to the best existing method in accuracy, while it is computationally more
efficient than competing approaches.
Keywords: importance sampling, covariate shift adaptation, novelty detection, regularization
path, leave-one-out cross validation

1. Introduction
In the context of importance sampling (Fishman, 1996), the ratio of two probability density functions is called the importance. The problem of estimating the importance is attracting a great deal ofattention these days since the importance can be used for various succeeding tasks such as covariate
shift adaptation or outlier detection.
Covariate Shift Adaptation: Covariate shift is a situation in supervised learning where
the distributions of inputs change between the training and test phases but the conditional distribution of outputs given inputs remains unchanged (Shimodaira, 2000;
Quiñonero-Candela et al., 2008). Covariate shift is conceivable in many real-world∗. A
MATLAB R
orR
implementationoftheproposedimportancerithm,
unconstrained
Least-Squares
Importance
Fitting(uLSIF),is
http://sugiyama-www.cs.titech.ac.jp/∼sugi/software/uLSIF/.

c 2009 Takafumi Kanamori, Shohei Hido and Masashi Sugiyama.

estimation
available

algofromK ANAMORI , H IDO AND S UGIYAMA

applications such as bioinformatics (Baldi and Brunak, 1998; Borgwardt et al., 2006),
brain-computer interfaces (Wolpaw et al., 2002; Sugiyama et al., 2007), robot control(Sutton and Barto, 1998; Hachiya et al., 2008), spam filtering (Bickel and Scheffer,2007), and econometrics (Heckman, 1979). Under covariate shift, standard learningtechniques such as maximum likelihood estimation or cross-validation are biased andtherefore unreliable—the bias caused by covariate shift can be compensated by weighting the loss function according to the importance (Shimodaira, 2000; Zadrozny, 2004;
Sugiyama and Müller, 2005; Sugiyama et al., 2007; Huang et al., 2007; Bickel et al.,2007).
Outlier Detection: The outlier detection task addressed here is to identify irregularsamples in a validation data set based on a model data set that only contains regularsamples (Schölkopf et al., 2001; Tax and Duin, 2004; Hodge and Austin, 2004; Hido
et al., 2008). The values of the importance for regular samples are close to one, while
those for outliers tend to be significantly deviated from one. Thus the values of theimportance could be used as an index of the degree of outlyingness.
Below, we refer to the two sets of samples as the training set and the test set.
A naive approach to estimating the importance is to first estimate the training and test density
functions from the sets of training and test samples separately, and then take the ratio of the estimated densities. However, density estimation is known to be a hard problem particularly in highdimensional cases if we do not have simple and good parametric density models (Vapnik, 1998;
Härdle et al., 2004). In practice, such an appropriate parametric model may not be available andtherefore this naive approach is not so effective.
To cope with this problem, direct importance estimation methods which do not involve density estimation have been developed recently. The kernel mean matching (KMM) method (Huang
et al., 2007) directly gives estimates of the importance at the training inputs by matching the two
distributions efficiently based on a special property of universal reproducing kernel Hilbert spaces(Steinwart, 2001). The optimization problem involved in KMM is a convex quadratic program, so
the unique global optimal solution can be obtained using a standard optimization software. However, the performance of KMM depends on the choice of tuning parameters such as the kernel parameter and the regularization parameter. For the kernel parameter, a popular heuristic of using themedian distance between samples as the Gaussian width could be useful in some cases (Schölkopf
and Smola, 2002; Song et al., 2007). However, there seems no strong justification for this heuristic
and the choice of other tuning parameters is still open.
A probabilistic classifier that separates training samples from test samples can be used for directly estimating the importance, for example, a logistic regression (LogReg) classifier (Qin, 1998;
Cheng and Chu, 2004; Bickel et al., 2007). Maximum likelihood estimation of LogReg models
can be formulated as a convex optimization problem, so the unique global optimal solution can beobtained. Furthermore, since the LogReg-based method only involves a standard supervised classification problem, the tuning parameters such as the kernel width and the regularization parametercan be optimized based on the standard cross-validation procedure. This is a very useful property
in practice.
The Kullback-Leibler importance estimation procedure (KLIEP) (Sugiyama et al., 2008b; Nguyen
et al., 2008) also directly gives an estimate of the importance function by matching the two distributions in terms of the Kullback-Leibler divergence (Kullback and Leibler, 1951). The optimization
1392

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

problem involved in KLIEP is convex, so the unique global optimal solution—which tends to besparse—can be obtained, when linear importance models are used. In addition, the tuning parameters in KLIEP can be optimized based on a variant of cross-validation.
As reviewed above, LogReg and KLIEP are more advantageous than KMM since the tuning
parameters can be objectively optimized based on cross-validation. However, optimization procedures of LogReg and KLIEP are less efficient in computation than KMM due to high non-linearityof the objective functions to be optimized—more specifically, exponential functions induced by the
LogReg model or the log function induced by the Kullback-Leibler divergence. The purpose ofthis paper is to develop a new importance estimation method that is equipped with a build-in model
selection procedure as LogReg and KLIEP and is computationally more efficient than LogReg and
KLIEP.
Our basic idea is to formulate the direct importance estimation problem as a least-squares
function fitting problem. This formulation allows us to cast the optimization problem as a convex quadratic program, which can be efficiently solved using a standard quadratic program solver.
Cross-validation can be used for optimizing the tuning parameters such as the kernel width or the
regularization parameter. We call the proposed method least-squares importance fitting (LSIF).
We further show that the solutions of LSIF is piecewise linear with respect to the ℓ1 -regularization
parameter and the entire regularization path (that is, all solutions for different regularization parameter values) can be computed efficiently based on the parametric optimization technique (Best,
1982; Efron et al., 2004; Hastie et al., 2004). Thanks to this regularization path tracking algorithm,
LSIF is computationally efficient in model selection scenarios. Note that in the regularization pathtracking algorithm, we can trace the solution path without a quadratic program solver—we just needto compute matrix inverses.
LSIF is shown to be efficient in computation, but it tends to share a common weakness of regularization path tracking algorithms, that is, accumulation of numerical errors (Scheinberg, 2006).
The numerical problem tends to be severe if there are many change points in the regularization
path. To cope with this problem, we develop an approximation algorithm in the same least-squares
framework. The approximation version of LSIF, which we call unconstrained LSIF (uLSIF), allows
us to obtain the closed-form solution that can be computed just by solving a system of linear equations. Thus uLSIF is numerically stable when regularized properly. Moreover, the leave-one-out
cross-validation score for uLSIF can also be computed analytically (cf. Wahba, 1990; Cawley and
Talbot, 2004), which significantly improves the computational efficiency in model selection scenarios. We experimentally show that the accuracy of uLSIF is comparable to the best existing method
while its computation is faster than other methods in covariate shift adaptation and outlier detection
scenarios.
Our contributions in this paper are summarized as follows. A proposed density-ratio estimation method, LSIF, is equipped with cross-validation (which is an advantage over KMM) and iscomputationally efficient thanks to regularization path tracking (which is an advantage over KLIEP
and LogReg). Furthermore, uLSIF is computationally even more efficient since its solution andleave-one-out cross-validation score can be computed analytically in a stable manner. The proposed
methods, LSIF and uLSIF, are similar in spirit to KLIEP, but the loss functions are different: KLIEP
uses the log loss while LSIF and uLSIF use the squared loss. The difference of the log functions
allows us to improve computational efficiency significantly.
The rest of this paper is organized as follows. In Section 2, we propose a new importanceestimation procedure based on least-squares fitting (LSIF) and show its theoretical properties. In
1393

K ANAMORI , H IDO AND S UGIYAMA

Section 3, we develop an approximation algorithm (uLSIF) which can be computed efficiently. In
Section 4, we illustrate how the proposed methods behave using a toy data set. In Section 5, we discuss the characteristics of existing approaches in comparison with the proposed methods and showthat uLSIF could be a useful alternative to the existing methods. In Section 6, we experimentallycompare the performance of uLSIF and existing methods. Finally in Section 7, we summarize ourcontributions and outline future prospects. Those who are interested in practical implementation
may skip the theoretical analyses in Sections 2.3, 3.2, and 3.3.

2. Direct Importance Estimation
In this section, we propose a new method of direct importance estimation.
2.1 Formulation and Notation
Let D ⊂ (Rd ) be the data domain and suppose we are given independent and identically distributed
tr(i.i.d.) training samples {xtri }ni=1from a training distribution with density ptr (x) and i.i.d. test samples
nte
te{x j } j=1 from a test distribution with density pte (x):i.i.d.

tr∼ ptr (x),{xtri }ni=1i.i.d.

te∼ pte (x).
{xtej }nj=1

We assume that the training density is strictly positive, that is,
ptr (x) > 0 for all x ∈ D .
tetr:and {xtej }nj=1
The goal of this paper is to estimate the importance w(x) from {xtri }ni=1w(x) =

pte (x),ptr (x)which is non-negative by definition. Our key restriction is that we want to avoid estimating densities
pte (x) and ptr (x) when estimating the importance w(x).
2.2 Least-squares Approach to Direct Importance Estimation
Let us model the importance w(x) by the following linear model:bb = ∑ αℓ ϕℓ (x),
w(x)

(1)

ℓ=1

where α = (α1 , α2 , . . . , αb )⊤ are parameters to be learned from data samples, ⊤ denotes the transposeof a matrix or a vector, and {ϕℓ (x)}bℓ=1 are basis functions such thatϕℓ (x) ≥ 0 for all x ∈ D and for ℓ = 1, 2, . . . , b.
tetr
, for example,and {xtej }nj=1
Note that b and {ϕℓ (x)}bℓ=1 could be dependent on the samples {xtri }ni=1bkernel models are also allowed. We explain how the basis functions {ϕℓ (x)}ℓ=1 are chosen in
Section 2.5.

1394

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATIONb so that the following squared error J0
We determine the parameters {αℓ }bℓ=1 in the model w(x)is minimized:
1b − w(x))2 ptr (x)dx
(w(x)
2
Z
Z
Z
1
1
b 2 ptr (x)dx − w(x)w(x)pb=
w(x)2 ptr (x)dxw(x)
tr (x)dx +
2
2
Z
Z
Z
1
1
2
b
b
w(x)2 ptr (x)dx,w(x)ptr (x)dx − w(x)p=te (x)dx +22J0 (α) =Z

where in the second term the probability density ptr (x) is canceled with that included in w(x).
The squared loss J0 (α) is defined as the expectation under the probability of training samples. In
covariate shift adaptation (see Section 6.2) and outlier detection (see Section 6.3), the importancevalues on the training samples are used. Thus, the definition of J0 (α) well agrees with our goal.
The last term of J0 (α) is a constant and therefore can be safely ignored. Let us denote the first
two terms by J:1b 2 ptr (x)dx − w(x)pbJ(α) =w(x)te (x)dx
2

Z
 b

Z
1 b
αℓ αℓ′
ϕℓ (x)ϕℓ′ (x)ptr (x)dx − ∑ αℓ
ϕℓ (x)pte (x)dx=
2 ℓ,ℓ∑
′ =1
ℓ=1
Z

Z

1
= α⊤ Hα − h⊤ α,2(2)where H is the b × b matrix with the (ℓ, ℓ′ )-th element
Hℓ,ℓ′ =Z

ϕℓ (x)ϕℓ′ (x)ptr (x)dx,(3)and h is the b-dimensional vector with the ℓ-th elementhℓ =Zϕℓ (x)pte (x)dx.

Approximating the expectations in J by empirical averages, we obtain
ntr
1 nte
b = 1 ∑ w(xb tri )2 −b tej )J(α)∑ w(x2ntr i=1
nte j=11 b=αℓ αℓ′
2 ℓ,ℓ∑
′ =1

!
b
1 ntr
tr
tr
ϕℓ (xi )ϕℓ′ (xi ) − ∑ αℓ∑
ntr i=1
ℓ=1

!
1 nte
te
∑ ϕℓ (x j )nte j=1⊤
1 b= α⊤ Hα−bh α,
2b is the b × b matrix with the (ℓ, ℓ′ )-th element
where H
ntr

bℓ,ℓ′ = 1 ∑ ϕℓ (xtri )ϕℓ′ (xtri ),H
ntr i=1
1395(4)

K ANAMORI , H IDO AND S UGIYAMA

and b
h is the b-dimensional vector with the ℓ-th element
nte

1
b
hℓ =
∑ ϕℓ (xtej ).
nte j=1(5)

Taking into account the non-negativity of the importance function w(x), we can formulate our optimization problem as follows.


⊤
1 ⊤b
⊤
b
min
α Hα − h α + λ1b α
2α∈Rbsubject to α ≥ 0b ,(6)
where 0b and 1b are the b-dimensional vectors with all zeros and ones, respectively; the vector
inequality α ≥ 0b is applied in the element-wise manner, that is, αℓ ≥ 0 for ℓ = 1, 2, . . . , b. In Eq. (6),
we included a penalty term λ1⊤
b α for regularization purposes, where λ (≥ 0) is a regularization
parameter. The above is a convex quadratic programming problem and therefore the unique global
optimal solution can be computed efficiently by a standard optimization package. We call thismethod Least-Squares Importance Fitting (LSIF).
We can also use the ℓ2 -regularizer α⊤ α instead of the ℓ1 -regularizer 1⊤
b α without changing thecomputational property. However, using the ℓ1 -regularizer would be more advantageous since thesolution tends to be sparse (Williams, 1995; Tibshirani, 1996; Chen et al., 1998). Furthermore, as
shown in Section 2.6, the use of the ℓ1 -regularizer allows us to compute the entire regularization
path efficiently (Best, 1982; Efron et al., 2004; Hastie et al., 2004). The ℓ2 -regularization method
will be used for theoretical analysis in Section 3.3.
2.3 Convergence Analysis of LSIF
b of the LSIF algorithm;
Here, we theoretically analyze the convergence property of the solution α
practitioners may skip this theoretical analysis.
b (λ) be the solution of the LSIF algorithm with regularization parameter λ, and let α∗ (λ)
Let αbe the optimal solution of the ‘ideal’ problem:


1 ⊤
⊤
⊤
min
α Hα − h α + λ1b α
2α∈Rbsubject to α ≥ 0b .
(7)
Below, we theoretically investigate the learning curve (Amari et al., 1992) of LSIF, that is, web (λ)) and J(α∗ (λ)) in terms of the expectation over all possible
elucidate the relation between J(αtraining and test samples as a function of the number of samples.
Let E be the expectation over all possible training samples of size ntr and all possible test samples of size nte . Let A ⊂ {1, 2, . . . , b} be the set of active indices (Boyd and Vandenberghe, 2004),
that is,A = {ℓ | α∗ℓ (λ) = 0, ℓ = 1, 2, . . . , b}.
For the active set A = { j1 , j2 , . . . , j|A | } with j1 < j2 < · · · < j|A | , let E be the |A |×b indicator matrix
with the (i, j)-th element(
1 j = ji ,Ei, j =
0 otherwise.
1396

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATIONb be the active set of αb (λ):
Similarly, let A

b = {ℓ | αb ℓ (λ) = 0, ℓ = 1, 2, . . . , b}.
A

b = {bb |×b indicator matrix
b be the |A
For the active set A
j1 , bj2 , . . . , b
j|Ab | } with b
j1 < bj2 < · · · < bj|Ab | , let Ewith the (i, j)-th element similarly defined by
(
1 j=b
ji ,
Ebi, j =(8)0 otherwise.

First, we show the optimality condition of (6) which will be used in the following theoreticalanalyses. The Lagrangian of the optimization problem (6) is given as
⊤
1 b
⊤
−b
h α + λ1⊤
L(α, ξ) = α⊤ Hαb α − ξ α,2where ξ is the b-dimensional Lagrange multiplier vector. Then the Karush-Kuhn-Tucker (KKT)
conditions (Boyd and Vandenberghe, 2004) are expressed as follows:b −bHαh + λ1b − ξ = 0b ,(9)ξℓ αℓ = 0 for ℓ = 1, 2, . . . , b.

(10)

α ≥ 0b ,ξ ≥ 0b ,′b |-dimensional vector with the i-th element being the bLet bξ (λ) be the |Aji -th element of b
ξ(λ):
′

b
ξ′i (λ) = b
ξ bji (λ),

b |.
i = 1, . . . , |A

(11)

b be
We assume that b
ξ (λ) only contains non-zero elements of b
ξ(λ). Let G
!
b
b⊤
H
−
E
b=
G
b Ob b ,−E|A |×|A |b | × |Ab | matrix with all zeros. Then Eqs. (9) and (10) are together expressed
where O|Ab |×|Ab | is the |Ain a matrix form as
!
!
b
b (λ)
α
h
−
λ1
b
b ′
.
(12)
=
G
b
0|Ab |ξ (λ)
b we have the following lemma:
Regarding the matrix G,

b is invertible if H
b is invertible.
Lemma 1 The matrix G

b is invertible.
The proof of the above lemma is given in Appendix A. Below, we assume that H−1bb
Then the inverse of G exists and multiplying G from the left-hand side of Eq. (12) yields
!


b (λ)
α
−1 b
h
−
λ1
b
b
′
.
(13)
=G
b
0|A |
ξ (λ)
1397

K ANAMORI , H IDO AND S UGIYAMA

The following inversion formula holds for block matrices (Petersen and Pedersen, 2007):
−1  −1


M1 M2
M 1 + M −1
M 2 M −1
M 3 M 1−1 −M −1
M 2 M −1
1
0
1
0
=
,
(14)
−1
M3 M4
−M −1
M −1
0 M3M10where
M 0 = M 4 − M 3 M −1
1 M2.
Applying Eq. (14) to Eq. (13), we haveb is defined bywhere Abbb (λ) = A(αh − λ1b ),b=Hb −1 − Hb −1 Eb ⊤ (EbH
b −1 Eb ⊤ )−1 EbHb −1 .
A

(15)
(16)

When the Lagrange multiplier vector satisfiesξ∗ℓ (λ) > 0 for all ℓ ∈ A ,(17)

we say that the strict complementarity condition is satisfied (Bertsekas et al., 2003). An importantconsequence of strict complementarity is that the optimal solution and the Lagrange multipliers ofconvex quadratic problems are uniquely determined. Then we have the following theorem.
Theorem 2 Let P be the probability over all possible training samples of size ntr and test samples ofsize nte . Let ξ∗ (λ) be the Lagrange multiplier vector of the problem (7) and suppose ξ∗ (λ) satisfies
the strict complementarity condition (17). Then, there exists a positive constant c > 0 and a natural
number N such that for min{ntr , nte } ≥ N,
b 6= A ) < e−c min{ntr ,nte } .
P(AThe proof of the above theorem is given in Appendix B. Theorem 2 shows that the probabilityb of the empirical problem (6) is different from the active set A of the idealthat the active set A
b = A in practice.
problem (7) is exponentially small. Thus we may regard Ab
Let A be the ‘ideal’ counterpart of A:
A = H −1 − H −1 E ⊤ (EH −1 E ⊤ )−1 EH −1 ,and let Cw,w′ be the b × b covariance matrix with the (ℓ, ℓ′ )-th element being the covariance between
w(x)ϕℓ (x) and w′ (x)ϕℓ′ (x) under ptr (x). Letbw∗ (x) = ∑ α∗ℓ (λ)ϕℓ (x),ℓ=1bv(x) = ∑ [A1b ]ℓ ϕℓ (x).
ℓ=1Letf (n) = ω(g(n))
denote that f (n) asymptotically dominates g(n); more precisely, for all C > 0, there exists n0 suchthat
|Cg(n)| < | f (n)| for all n > n0 .
Then we have the following theorem.
1398

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

Theorem 3 Assume that(a) The optimal solution of the problem (7) satisfies the strict complementarity condition (17).
(b) ntr and nte satisfy

nte = ω(n2tr ).

(18)

Then, for any λ ≥ 0, we have 
11b (λ))] = J(α (λ)) +E[J(αtr(A(Cw∗ ,w∗ − 2λCw∗ ,v )) + o
.
2ntr
ntr
∗

(19)

The proof of the above theorem is given in Appendix C. This theorem elucidates the learning
curve of LSIF up to the order of n−1
tr . In Section 2.4.1, we discuss practical implications of thistheorem.
2.4 Model Selection for LSIFThe practical performance of LSIF depends on the choice of the regularization parameter λ andbasis functions {ϕℓ (x)}bℓ=1 (which we refer to as a model). Since our objective is to minimize thecost function J defined in Eq. (2), it is natural to determine the model such that J is minimized.
However, the value of the cost function J is inaccessible since it includes the expectation overunknown probability density functions ptr (x) and pte (x). The value of the empirical cost Jb may be
regarded as an estimate of J, but this is not useful for model selection purposes since it is heavily
biased—the bias is caused by the fact that the same samples are used twice for learning the parameter
α and estimating the value of J. Below, we give two practical methods of estimating the value of J
in more precise ways.
2.4.1 I NFORMATION C RITERION
Inhthe samei way as Theorem 3, we can obtain an asymptotic expansion of the empirical costbαb (λ)) as follows:
E J(
 
1
1
∗
b
b (λ))] = J(α (λ)) −tr(A(Cw∗ ,w∗ + 2λCw∗ ,v )) + o
.
E[J(α
2ntr
ntr

(20)

Combining Eqs. (19) and (20), we have

 
1
1
bα
b (λ))] = E[J(b (λ))] + tr(ACw∗ ,w∗ ) + o
.
E[J(α
ntr
ntr

From this, we can immediately obtain an information criterion (Akaike, 1974; Konishi and Kitagawa, 1996) for LSIF:
1 bbbαb (λ)) + tr(A
Jb(IC) = J(
Cw,
bw
b ),
ntr
b is defined by Eq. (16). E
b w,w′ is the b × b covariance matrix withb is defined by Eq. (8) and C
where A′tr
. Sincethe (ℓ, ℓ )-th element being the covariance between w(x)ϕℓ (x) and w′ (x)ϕℓ′ (x) over {xtri }ni=1bb
A and Cw,
bw
b are consistent estimators of A and C w∗ ,w∗ , the above information criterion is unbiased up
to the order of n−1
tr .
1399

K ANAMORI , H IDO AND S UGIYAMAbC
b w,
Note that the term tr(Abwb ) may be interpreted as the effective dimension of the model (Moody,
b w,
b =C
b = 1, we have H
1992). Indeed, when w(x)bwb and thus
−1

−1

⊤ −1
b
b w,
bC
b w,
b b wb E ⊤ (ECtr(Abwb ) = tr(I b ) − tr(E Cw,
bw
b E ) ) = b − |A |,b (λ) lies.
which is the dimension of the face on which α2.4.2 C ROSS - VALIDATION

Although the information criterion derived above is more accurate than just a naive empirical estimator, its accuracy is guaranteed only asymptotically. Here, we employ cross-validation for estimatingb ), which has an accuracy guarantee for finite samples.
J(α
tr
te
are divided into R disjoint subsets
First, the training samples {xtri }ni=1
and test samples {xtej }nj=1
te
R
tr
R
bXrtr ,Xrte (x) is obtained using{Xr }r=1 and {Xr }r=1 , respectively. Then an importance estimate w{X jtr } j6=r and {X jte } j6=r (that is, without Xrtr and Xrte ), and the cost J is approximated using theheld-out samples Xrtr and Xrte as(CV)
JbX tr ,X te =
r

r

1
1
bXrtr ,Xrte (xtr )2 − te ∑ w
bXrtr ,Xrte (xte ).
w
2|Xrtr | xtr∑|
X|r xte ∈X te
∈X tr
r

r

This procedure is repeated for r = 1, 2, . . . , R and its average Jb(CV) is used as an estimate of J:
R

1
(CV)
Jb(CV) = ∑ JbX tr ,X te .
R r=1 r r

We can show that Jb(CV) gives an almost unbiased estimate of the true cost J, where the ‘almost’-nesscomes from the fact that the number of samples is reduced in the cross-validation procedure due todata splitting (Luntz and Brailovsky, 1969; Wahba, 1990; Schölkopf and Smola, 2002).
Cross-validation would be more accurate than the information criterion for finite samples. However, it is computationally more expensive than the information criterion since the learning procedure should be repeated R times.
2.5 Heuristics of Basis Function Design for LSIF
A good model may be chosen by cross-validation or the information criterion, given that a family ofpromising model candidates is prepared. As model candidates, we propose using a Gaussian kerneltemodel centered at the test points {xtej }nj=1, that is,nte

b = ∑ αℓ Kσ (x, xte
w(x)
ℓ ),
ℓ=1

where Kσ (x, x′ ) is the Gaussian kernel with kernel width σ:kx − x′ k2
Kσ (x, x′ ) = exp −
.
2σ2

(21)

te
as the Gaussian centers, not the training points
The reason why we chose the test points {xtej }nj=1ntr
tr{xi }i=1 , is as follows (Sugiyama et al., 2008b). By definition, the importance w(x) tends to take1400

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

large values if the training density ptr (x) is small and the test density pte (x) is large; conversely,
w(x) tends to be small (that is, close to zero) if ptr (x) is large and pte (x) is small. When a function
is approximated by a Gaussian kernel model, many kernels may be needed in the region where theoutput of the target function is large; on the other hand, only a small number of kernels would beenough in the region where the output of the target function is close to zero. Following this heuristic,
we allocate many kernels at high test density regions, which can be achieved by setting the Gaussiante.
centers at the test points {xtej }nj=1
tetr. However,and {xtej }nj=1
Alternatively, we may locate (ntr +nte ) Gaussian kernels at both {xtri }ni=1
in our preliminary experiments, this did not further improve the performance, but just slightly inte
as Gaussiancreased the computational cost. When nte is large, just using all the test points {xtej }nj=1
centers is already computationally rather demanding. To ease this problem, we practically proposeteas Gaussian centers for computational efficiency, that is,using a subset of {xtej }nj=1bb = ∑ αℓ Kσ (x, cℓ ),
w(x)

(22)

ℓ=1

te
where cℓ , ℓ = 1, 2, . . . , b are template points randomly chosen from {xtej }nj=1
without replacement
and b (≤ nte ) is a prefixed number. In the rest of this paper, we usually fix the number of template
points at
b = min(100, nte ),and optimize the kernel width σ and the regularization parameter λ by cross-validation with gridsearch.
2.6 Entire Regularization Path for LSIF
b is piecewise linear with respect to the regularization parameter
We can show that the LSIF solution αλ (see Appendix D). Therefore, the regularization path (that is, solutions for all λ) can be computed
efficiently based on the parametric optimization technique (Best, 1982; Efron et al., 2004; Hastie
et al., 2004).
A basic idea of regularization path tracking is to check violation of the KKT conditions—which
are necessary and sufficient for optimality of convex programs—when the regularization parameter
λ is changed. The KKT conditions of LSIF are summarized in Section 2.3. The strict complementarity condition (17) assures the uniqueness of the optimal solution for a fixed λ, and thus theuniqueness of the regularization path. A pseudo code of the regularization path tracking algorithm
for LSIF is described in Figure 1—its detailed derivation is summarized in Appendix D. Thanks to
the regularization path algorithm, LSIF is computationally efficient in model selection scenarios.
The pseudo code implies that we no longer need a quadratic programming solver for obtaining
the solution of LSIF—just computing matrix inverses is enough. Furthermore, the regularization
path algorithm is computationally more efficient when the solution is sparse, that is, most of theelements are zero since the number of change points tends to be small for such sparse solutions.
Even though the regularization path tracking algorithm is computationally efficient, it tends tobe numerically unreliable, as we experimentally show in Section 4. This numerical instability isb When G
b is nearly singular, it is not easy to accurately
caused by near singularity of the matrix G.
obtain the solutions u, v in Figure 1, and therefore the change point λτ+1 cannot be accurately computed. As a result, we cannot accurately update the active set of the inequality constraints and thus1401

K ANAMORI , H IDO AND S UGIYAMAb and b
Input: Hh% see Eqs. (4) and (5) for the definitionsb (λ) for λ ≥ 0
Output: entire regularization path ατ ←− 0;k ←− argmaxi {bhi | i = 1, 2, . . . , b};b
λτ ←− hk ;b ←− {1, 2, . . . , b}\{k};Ab (λτ ) ←− 0b ;α
% the vector with all zeros
While λτ > 0b ←− O b ;E% the matrix with all zeros|A |×bb
For i = 1, 2, . . . , |A |b = {b
Eb b ←− 1; % A
j1 , bj2 , . . . , b
jb |b
j1 < bj2 < · · · < b
jb}
|A |

i, ji

end

|A |

!
b
b⊤
H
−
E
b ←−
G
b Ob b ;−E|!
A |×|A |bb −1 h ;u ←− G
0|Ab |!
1−1bb;v ←− G
0|Ab |

end

If v ≤ 0b+|Ab | % the final intervalλτ+1 ←− 0;b (λτ+1 ) ←− (u1 , u2 , . . . , ub )⊤ ;αelse % an intermediate intervalb |};k ←− argmaxi {ui /vi | vi > 0, i = 1, 2, . . . , b + |A
λτ+1 ←− max{0, uk /vk };b (λτ+1 ) ←− (u1 , u2 , . . . , ub )⊤ − λτ+1 (v1 , v2 , . . . , vb )⊤ ;α
If 1 ≤ k ≤ bb ←− Ab ∪ {k};Aelseb ←− Ab \{ bAjk−b };endendτ ←− τ + 1;

b (λ) ←−α(0b

λτ+1 −λλ−λτbb
λτ+1 −λτ α(λτ ) + λτ+1 −λτ α(λτ+1 )if λ ≥ λ0
if λτ+1 ≤ λ ≤ λτ

Figure 1: Pseudo code for computing the entire regularization path of LSIF. When the computationb −1 is numerically unstable, we may add small positive diagonals to H
b for stabilizationof G
purposes.

1402

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATIONb (λ) becomes unreliable; furthermore, such numerical error tends to be accuthe obtained solution αmulated through the path-tracking process. This instability issue seems to be a common pitfall ofsolution path tracking algorithms in general (see Scheinberg, 2006).
b tends to be nearly singular
When the Gaussian width σ is very small or very large, the matrix Hband thus the matrix G also becomes nearly singular. On the other hand, when the Gaussian width σ isb is well-conditioned
not too small or too large compared with the dispersion of samples, the matrix G
and therefore the path-following algorithm would be stable and reliable.

3. Approximation Algorithm
Within the quadratic programming formulation, we have proposed a new importance estimation
procedure LSIF and showed its theoretical properties. We also gave a regularization path tracking
algorithm that can be computed efficiently. However, as we experimentally show in Section 4, it
tends to suffer from a numerical problem and therefore is not practically reliable. In this section, wegive a practical alternative to LSIF which gives an approximate solution to LSIF in a computationally efficient and reliable manner.
3.1 Unconstrained Least-squares Formulation
The approximation idea we introduce here is very simple: we ignore the non-negativity constraintof the parameters in the optimization problem (6). This results in the following unconstrained
optimization problem.


⊤
1 ⊤b
λ ⊤
b
min β Hβ − h β + β β .
(23)
2β∈Rb 2In the above, we included a quadratic regularization term β⊤ β/2, instead of the linear one 1⊤b β since
the linear penalty term does not work as a regularizer without the non-negativity constraint. Eq. (23)
is an unconstrained convex quadratic program, so the solution can be analytically computed aseb + λI b )−1bβ(λ) = (Hh,where I b is the b-dimensional identity matrix. Since we dropped the non-negativity constraint β ≥0b , some of the learned parameters could be negative. To compensate for this approximation error,
we modify the solution bybβ(λ) = max(0b , eβ(λ)),

where the ‘max’ operation for a pair of vectors is applied in the element-wise manner. This is thesolution of the approximation method we propose in this section.
An advantage of the above unconstrained formulation is that the solution can be computed justby solving a system of linear equations. Therefore, its computation is stable when λ is not too small.
We call this method unconstrained LSIF (uLSIF). Due to the ℓ2 regularizer, the solution tends to
be close to 0b to some extent. Thus, the effect of ignoring the non-negativity constraint may not be
so strong—later, we analyze the approximation error both theoretically and experimentally in moredetail in Sections 3.3 and 4.5.
Note that LSIF and uLSIF differ only in parameter learning. Thus, the basis design heuristic of
LSIF given in Section 2.5 is also valid for uLSIF.
1403

K ANAMORI , H IDO AND S UGIYAMA

3.2 Convergence Analysis of uLSIF
Here, we theoretically analyze the convergence property of the solution bβ(λ) of the uLSIF algorithm; practitioners may skip Sections 3.2 and 3.3.
Let β◦ (λ) be the optimal solution of the ‘ideal’ version of the problem (23):


λ ⊤
1 ⊤
⊤
min β Hβ − h β + β β .
2β∈Rb 2
Then the ideal solution β∗ (λ) is given byβ∗ (λ) = max(0b , β◦ (λ)),β◦ (λ) = B−1
λ h,

(24)

Bλ = H + λI b .

Below, we theoretically investigate the learning curve of uLSIF.
Let B ⊂ {1, 2, . . . , b} be the set of negative indices of β◦ (λ), that is,

B = {ℓ | β◦ℓ (λ) < 0, ℓ = 1, 2, . . . , b},e ⊂ {1, 2, . . . , b} be the set of negative indices of eand Bβ(λ), that is,

e = {ℓ | eBβℓ (λ) < 0, ℓ = 1, 2, . . . , b}.

Then we have the following theorem.

Theorem 4 Assume that β◦ℓ (λ) 6= 0 for ℓ = 1, 2, . . . , b. Then, there exists a positive constant c anda natural number N such that for min{ntr , nte } ≥ N,e ) < e−c min{ntr ,nte } .
P(B 6= B

The proof of the above theorem is given in Appendix E. The assumption that β◦ℓ (λ) 6= 0 forℓ = 1, 2, . . . , b corresponds to the strict complementarity condition (17) in LSIF. Theorem 4 shows
e is different from B is exponentially small. Thus we may regard B
e = B inthat the probability that B
practice.
Let D be the b-dimensional diagonal matrix with the ℓ-th diagonal element(
0 ℓ ∈ B,Dℓ,ℓ =
1 otherwise.
Letbw◦ (x) = ∑ β◦ℓ (λ)ϕℓ (x),ℓ=1b∗u(x) = ∑ [B−1λ D(Hβ (λ) − h)]ℓ ϕℓ (x).
ℓ=1Then we have the following theorem.
1404

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

Theorem 5 Assume that(a) β◦ℓ (λ) 6= 0 for ℓ = 1, 2, . . . , b.
(b) ntr and nte satisfy Eq. (18).
Then, for any λ ≥ 0, we have

 
1
1
∗
−1
−1
−1
b
◦
◦
◦
E[J(β(λ))] = J(β (λ)) +.
tr(Bλ DHDBλ Cw ,w + 2Bλ Cw ,u ) + o
2ntr
ntr

(25)

The proof of the above theorem is given in Appendix F. Theorem 5 elucidates the learning curveof uLSIF up to the order of n−1
tr . An information criterion may be obtained in the same way as
Section 2.4.1. However, as shown in Section 3.4, we can have a closed-form expression of theleave-one-out cross-validation score for uLSIF, which would be practically more useful. For thisreason, we do not go into the detail of information criterion.
3.3 Approximation Error Bounds for uLSIFThe uLSIF method is introduced as an approximation of LSIF. Here, we theoretically evaluate theb (λ). More specifically, we use
difference between the uLSIF solution bβ(λ) and the LSIF solution α
the following normalized L2 -norm on the training samples as the difference measure and derive its
upper bounds:
r
2

tr
b (λ′ )) − w(xb tri ; αb tri ; b
infλ′ ≥0 n1tr ∑ni=1
w(x
β(λ))
diff(λ) =
,
(26)
ntr
b tri ; b
w(x
β(λ))
∑i=1

b α) is given by
where the importance function w(x;bb α) = ∑ αℓ ϕℓ (x).
w(x;
ℓ=1

In the theoretical analysis below, we assume
ntr

b tri ; bβ(λ)) 6= 0.
∑ w(x

i=1

For p ∈ N ∪ {∞}, let k · k p be the L p -norm, and let kαkHb be
kαkHb =pbα⊤ Hα,b is the b × b matrix defined by Eq. (4). Then we have the following theorem.
where H
Theorem 6 (Norm bound) Assume that all basis functions satisfy
0 < ϕℓ (x) ≤ 1.
1405

(27)

K ANAMORI , H IDO AND S UGIYAMA

Then we havekb
β(λ)kHbntr
b tri ; b
w(x
β(λ))
∑i=1


nte
1
b
·
,
≤ b2 1 +ntr
ntetr
λ minℓ ∑i=1 ϕℓ (xi ) minℓ ∑ j=1 ϕℓ (xtej )diff(λ) ≤

(28)
(29)

where b is the number of basis functions. The upper bound (29) is reduced as the regularization
parameter λ increases. For the Gaussian basis function model (22), the upper bound (29) is reduced
as the Gaussian width σ increases.
The proof of the above theorem is given in Appendix G. We call Eq. (28) the norm bound since
it is governed by the norm of bβ. Intuitively, the approximation error of uLSIF would small if λ iselarge since β ≥ 0 may not be severely violated due to the strong regularization effect. The upperbound (29) justifies this intuitive claim since the error bound tends to be small if the regularization
parameter λ is large. Furthermore, the upper bound (29) shows that for the Gaussian basis function
model (22), the error bound tends to be small if the Gaussian width σ is large. This is also intuitivesince the Gaussian basis functions are nearly flat when the Gaussian width σ is large—a differencebin parameters does not cause a significant change in the learned importance function w(x).
Fromthe above theorem, we expect that uLSIF is a nice approximation of LSIF when λ is large and σ islarge. In Section 4.5, we numerically investigate this issue.
Below, we give a more sophisticated bound on diff(λ). To this end, let us introduce an intermediate optimization problem defined by

⊤
1 ⊤b
λ ⊤
b
min γ Hγ − h γ + γ γ
2γ∈Rb 2subject to γ ≥ 0b ,(30)

which we refer to as LSIF with quadratic penalty (LSIFq). LSIFq bridges LSIF and uLSIF sincethe ‘goodness-of-fit’ part is the same as LSIF but the ‘regularization’ part is the same as uLSIF. Letbγ(λ) be the optimal solution of LSIFq (30). Based on the solution of LSIFq, we have the following
upper bound.
Theorem 7 (Bridge bound) For any λ ≥ 0, the following inequality holds:
diff(λ) ≤q
β(λ)kHb
λ kbγ(λ)k1 · kbγ(λ)k∞ − kbγ(λ)k22 + kbγ(λ) − b
ntr
b tri ; b
w(x
β(λ))
∑i=1

.

(31)

The proof of the above theorem is given in Appendix H. We call the above bound the bridgebound since the bridged estimator bγ(λ) plays a central role in the bound. Note that, in the bridgebound, the inside of the square root is assured to be non-negative due to Hölder’s inequality (seeAppendix H for detail). The bridge bound is generally much sharper than the norm bound (28), but
not always (see Section 4.5 for numerical examples).
1406

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

3.4 Efficient Computation of Leave-one-out Cross-validation Score for uLSIF
A practically important advantage of uLSIF over LSIF is that the score of leave-one-out crossvalidation (LOOCV) can be computed analytically—thanks to this property, the computational
complexity for performing LOOCV is the same order as just computing a single solution.
tetr, which generallyand {xtej }nj=1
In the current setup, we are given two sets of samples, {xtri }ni=1
have different sample size. For simplicity, we assume that ntr < nte and the i-th training sample xtri
te nteand the i-th test sample xte
i are held out at the same time; the test samples {x j } j=ntr +1 are always
used for importance estimation. Note that this assumption is only for the sake of simplicity; we canchange the order of test samples without sacrificing the computational advantages.
b(i) (x) be an estimate of the importance obtained without the i-th training sample xtri and the
Let w
i-th test sample xte
i . Then the LOOCV score is expressed as
1 ntr 1 (i) tr 2(i) teLOOCV =
∑ 2 (wb (xi )) − wb (xi ) .
ntr i=1(32)

Our approach to efficiently computing the LOOCV score is to use the Sherman-Woodbury-Morrison
formula (Golub and Loan, 1996) for computing matrix inverses: for an invertible square matrix A
and vectors u and v such that v⊤ A−1 u 6= −1,(A + uv⊤ )−1 = A−1 −A−1 uv⊤ A−1.
1 + v⊤ A−1 u(33)

Efficient approximation schemes of LOOCV have often been investigated under asymptotic
setups (Stone, 1974; Hansen and Larsen, 1996). On the other hand, we provide the exact LOOCV
score of uLSIF, which follows the same line as that of ridge regression (Hoerl and Kennard, 1970;
Wahba, 1990).
A pseudo code of uLSIF with LOOCV-based model selection is summarized in Figure 2—its
detailed derivation is described in Appendix I. Note that the basis-function design heuristic given
in Section 2.5 is used in the pseudo code, but the analytic form of the LOOCV score is available forany basis functions.

4. Illustrative Examples
In this section, we illustrate the behavior of LSIF and uLSIF using a toy data set.
4.1 Setup
Let the dimension of the domain be d = 1 and the training and test densities be
ptr (x) = N (x; 1, (1/2)2 ),pte (x) = N (x; 2, (1/4)2 ),where N (x; µ, σ2 ) denotes the Gaussian density with mean µ and variance σ2 . These densities are
depicted in Figure 3. The task is to estimate the importance w(x) = pte (x)/ptr (x).
1407

K ANAMORI , H IDO AND S UGIYAMAtr
te
Input: {xtri }ni=1
and {xtej }nj=1bOutput: w(x)b ←− min(100, nte ); n ←− min(ntr , nte );teRandomly choose b centers {cℓ }bℓ=1 from {xtej }nj=1
without replacement;
For each candidate of Gaussian
width
σ


kxtri − cℓ k2 + kxtri − cℓ′ k2
1 ntrb
for ℓ, ℓ′ = 1, 2, . . . , b;
Hℓ,ℓ′ ←−
∑ exp −
2
ntr i=1
2σ!
kxtej − cℓ k2
1 ntebfor ℓ = 1, 2, . . . , b;hℓ ←−∑ exp − 2σ2
nte j=1
kxtr − cℓ k2
tr
Xℓ,i←− exp − i 2
for i = 1, 2, . . . , n and ℓ = 1, 2, . . . , b;
2σ


te
2
kx − cℓ kteXℓ,i←− exp − i 2
for i = 1, 2, . . . , n and ℓ = 1, 2, . . . , b;
2σFor each candidate of regularization parameter λb ←− Hb + λ(ntr − 1) I b ;Bntr!
⊤ −1 tr
b
b X
h B
−1 ⊤
−1 tr
b
b
b
;
B0 ←− B h1n + B X diag⊤tr b −1 tr
ntr 1⊤
n − 1b (X ∗ B X )!
te b −1 tr
1⊤
−1 te
−1 tr
b (X ∗ B X )bb;
B1 ←− B X + B X diag⊤tr b −1 tr
ntr 1⊤
n − 1b (X∗ B X )ntr − 1B2 ←− max Ob×n ,(nte B0 − B1 ) ;ntr (nte − 1)
⊤
tr
⊤
te
⊤
wtr ←− (1b (X ∗ B2 )) ;wte ←− (1⊤b (X ∗ B2 )) ;w⊤ wtr 1⊤
wte
LOOCV(σ, λ) ←− tr
− n
;
2n
n
end
end
b, bλ) ←− argmin(σ,λ) LOOCV(σ, λ);
(σ

ntr
2
tr
2
tr
eℓ,ℓ′ ←− 1 ∑ exp − kxi − cℓ k + kxi − cℓ′ k
for ℓ, ℓ′ = 1, 2, . . . , b;Hb2
ntr i=12
σ!
nte
kxtej − cℓ k2
1e
for ℓ = 1, 2, . . . , b;hℓ ←−∑ exp − 2σb2
nte j=1
e + bλI b )−1eb ←− max(0b , (H
α
h);


b
kx − cℓ k2
b ℓ exp −b ←− ∑ α
w(x)
;
b2
2σ
ℓ=1

Figure 2: Pseudo code of uLSIF algorithm with LOOCV. B ∗ B′ denotes the element-wise multiplication of matrices B and B′ of the same size, thatis, the (i, j)-th element is given by
Bi, j B′i, j . For n-dimensional vectors b and b′ , diag bb′ denotes the n × n diagonal matrix
with i-th diagonal element bi /b′i . A MATLAB R or R implementation of uLSIF is available from http://sugiyama-www.cs.titech.ac.jp/∼sugi/software/uLSIF/.

1408

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

4.2 Importance Estimation
First, we illustrate the behavior of LSIF and uLSIF in importance estimation. We set the numberof training and test samples at ntr = 200 and nte = 1000, respectively. We use the Gaussian kernel
model (22), and the number of basis functions is set at b = 100. The centers of the kernel functionnte
are randomly chosen from the test points {xtei } j=1 without replacement (see Section 2.5).
We test different Gaussian widths σ and different regularization parameters λ. The following
two setups are examined:(A) λ is fixed at λ = 0.2 and σ is changed as 0.1 ≤ σ ≤ 1.0,
(B) σ is fixed at σ = 0.3 and λ is changed as 0 ≤ λ ≤ 0.5.
Figure 4 depicts the true importance andRits estimates obtained by LSIF and uLSIF, where allimportance functions are normalized so that w(x)dx = 1 for better comparison. Figures 4(a) andb
4(b) show that the estimated importance w(x)tends to be too peaky when the Gaussian width σ
is small, while it tends to be overly smoothed when σ is large. If the Gaussian width is chosen
appropriately, both LSIF and uLSIF seem to work reasonably well. As shown in Figures 4(c)and 4(d), the solutions of LSIF and uLSIF also significantly change when different regularization
parameters λ are used. Again, given that the regularization parameter is chosen appropriately, both
LSIF and uLSIF tend to perform well.
From the graphs, we also observe that model selection based on cross-validation works reasonably well for both LSIF (5-fold) and uLSIF (leave-one-out) to choose appropriate values of theGaussian width or the regularization parameter; this will be analyzed in more detail in Section 4.4.
4.3 Regularization Path
Next, we illustrate how the regularization path tracking algorithm for LSIF behaves. We set the
number of training and test samples at ntr = 50 and nte = 100, respectively. For better illustration,
we set the number of basis functions at a small value as b = 30 in the Gaussian kernel model (22)and use the Gaussian kernels centered at equidistant points in [0, 3] as basis functions.
We use the algorithm described in Figure 1 for regularization path tracking. Theoretically, theinequality λτ+1 < λτ is assured. In numerical computation, however, the inequality is occasionallyb for stabilization (seeviolated. In order to avoid this numerical problem, we slightly regularize Halso the caption of Figure 1).
Figure 5 depicts the values of the estimated coefficients {αℓ }bℓ=1 as functions of kαk1 for
σ = 0.1, 0.3, and 0.5. Note that small kαk1 corresponds to large λ. The figure indicates that theregularization parameter λ works as a sparseness controlling factor of the solution, that is, the larger(smaller) the value of λ (kαk1 ) is, the sparser the solution is.
The path following algorithm is computationally efficient and therefore practically very attractive. However, as the above experiments illustrate, the path following algorithm is numericallyb can ease to solve this problem, but this in turn results in accurather unstable. Modification of H
mulating numerical errors through the path tracking process. Consequently, the solutions for small
λ tend to be inaccurate. This problem becomes prominent if the number of change points in theregularization path is large.
1409

K ANAMORI , H IDO AND S UGIYAMA

1.5

ptr(x)
pte(x)

1.0

0.5

0.0

-1

0

1

2

3

x

Figure 3: The solid line is the probability density of training data, and the dashed line is the probability density of test data.

0.08

true importanceσ = 0.1
σ = 0.3 (chosen by cross validation)
σ = 1.0

0.00

0.04

w(x)

0.04
0.00

w(x)

0.08

true importanceσ = 0.1
σ = 0.4 (chosen by cross validation)
σ = 1.0

1.0

1.5

2.0

2.5

3.0

3.5

1.0

1.5

2.0

x

3.5

0.08

(b) uLSIF for λ = 0.2, σ = 0.1, 0.3, 1.0.

true importanceλ=0
λ = 0.09 (chosen by cross validation)
λ = 0.5

0.02
0.00

0.00

0.02

w(x)

0.04

0.06

true importanceλ=0
λ = 0.2 (chosen by cross validation)
λ = 0.5

0.06

3.0

0.04

0.08

(a) LSIF for λ = 0.2, σ = 0.1, 0.4, 1.0.

w(x)

2.5
x

1.0

1.5

2.0

2.5

3.0

3.5

1.0

x

1.5

2.0

2.5

3.0

3.5

x

(c) LSIF for σ = 0.3, λ = 0, 0.2, 0.5.

(d) uLSIF for σ = 0.3, λ = 0, 0.09, 0.5.

Figure 4: True and estimated importance functions obtained by LSIF and uLSIF for various different Gaussian widths σ and regularization parameters λ.
1410

0.0e+00

0

0.0e+00

2

4

coefficients !
5.0e+06
1.0e+07

coefficients !
6
8
10

12

1.5e+07

14

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

1.0e+072.0e+073.0e+0702||!||168||!||1101214(b) σ = 0.3.

0.0

0.2

coefficients !
0.4
0.6
0.8

1.0

(a) σ = 0.1.

4

0.2

0.4

0.6

0.8

||!||11.0

1.2(c) σ = 0.5.

Figure 5: Regularization path of LSIF: the values of the estimated coefficients {αℓ }bℓ=1 are depicted
as functions of the L1 -norm of the estimated parameter vector for σ = 0.1, 0.3, and 0.5.
Small kαk1 corresponds to large λ.

4.4 Cross-validation
Here we illustrate the behavior of the cross-validation scores of LSIF and uLSIF. We set the numberof training and test samples at ntr = 200 and nte = 1000, respectively. The number of template
points is b = 100 and the Gaussian kernel model (22) is used. The centers of the kernel functions
are randomly chosen from the test points as described in Section 4.2. The left column of Figure 6b ) over 50 runs for LSIF and its estimate by 5-fold CV (25,depicts the expectation of the true cost J(α50, and 75 percentiles are plotted in the figure) as functions of the Gaussian width σ for λ = 0.2, 0.5,
and 0.8. We used the regularization path tracking algorithm for computing the solutions of LSIF.
1411

5
J

0
−6

−5

−4

0

−2

J

2

10

4

6

15

K ANAMORI , H IDO AND S UGIYAMA

0.1

0.2

0.5

1.0

0.1

σ (Gaussian width)

0.2

0.5

1.0

σ (Gaussian width)(b) uLSIF with LOOCV (λ = 0.2).

J
−6

−3.0

−4

−2

−2.0

J

0

−1.0

2

4

0.0

(a) LSIF with 5CV (λ = 0.2).

0.1

0.2

0.5

1.0

0.1

σ (Gaussian width)

0.2

0.5

1.0

σ (Gaussian width)(d) uLSIF with LOOCV (λ = 0.5).

J
0.1

0.2

0.5

1.0

−5 −4 −3 −2 −1

−0.6

−0.4

J

−0.2

0

1

0.0

(c) LSIF with 5CV (λ = 0.5).

0.1

σ (Gaussian width)

0.2

0.5

1.0

σ (Gaussian width)(e) LSIF with 5CV (λ = 0.8).

(f) uLSIF with LOOCV (λ = 0.8).

Figure 6: The true cost J and its cross-validation estimate as functions of Gaussian width σ fordifferent values of λ. The solid line denotes the expectation of the true cost J over 50
runs, while ‘◦’ and error bars denote the 25, 50, and 75 percentiles of the cross-validation
score.
1412

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATIONThe right column shows the expected true cost and its LOOCV estimates for uLSIF in the samemanner.
The graphs show that overall CV gives reasonably good approximations of the expected cost,
although CV for LSIF with small λ and small σ is rather inaccurate due to numerical problems—the
solution path of LSIF is computed from λ = ∞ to λ = 0, and the numerical error is accumulated as
the tracking process approaches to λ = 0. This phenomenon seems problematic when σ is small.
4.5 Difference between LSIF and uLSIFIn Section 3.3, we analyzed the approximation error of uLSIF against LSIF. Here we numericallyinvestigate the behavior of the approximation error (26) as well as the norm bound (28) and thebridge bound (31). We set the number of training and test samples at ntr = 200 and nte = 1000,
respectively. The number of template points in the Gaussian kernel model (22) is set at b = 100.
The centers of the kernel functions are randomly chosen from the test points (see Section 4.2).
Figure 7 depicts the true approximation error as well as its upper bounds as functions of the regularization parameter λ; λ is varied from 0.001 to 10 and the three Gaussian widths σ = 0.1, 0.5, 1.0
are tested. The graphs show that when λ and σ are large, the approximation error tends to be small;this is in good agreement with the theoretical analysis given in Section 3.3. The bridge bound isfairly tight in the entire range and is sharper than the norm bound except when σ is small and λ islarge.
4.6 Summary
Through the numerical examples, we overall found that LSIF and uLSIF give qualitatively similar
results. However, the computation of the solution-path tracking algorithm for LSIF tends to be
numerically unstable, which can result in unreliable model selection performance. On the other
hand, only a system of linear equations needs to be solved in uLSIF, which turned out to be much
more stable than LSIF. Thus, uLSIF would be practically more reliable than LSIF.
Based on the above findings, we will focus on uLSIF in the rest of this paper.

5. Relation to Existing Methods
In this section, we discuss the characteristics of existing approaches in comparison with the proposed methods.
5.1 Kernel Density Estimator
The kernel density estimator (KDE) is a non-parametric technique to estimate a probability density
function p(x) from its i.i.d. samples {xk }nk=1 . For the Gaussian kernel (21), KDE is expressed as
pb(x) =n1
Kσ (x, xk ).
∑ntr (2πσ2 )d/2 k=1

The performance of KDE depends on the choice of the kernel width σ. The kernel width σ canbe optimized by likelihood cross-validation (LCV) as follows (Härdle et al., 2004): First, divide
the samples {xi }ni=1 into R disjoint subsets {Xr }Rr=1 . Then obtain a density estimate pbXk (x) from1413

K ANAMORI , H IDO AND S UGIYAMA

Difference measure
0.006
0.010

Norm boundBridge boundTrue difference

0.002

0.005

Difference measure
0.015
0.025

Norm boundBridge boundTrue difference

1e-03

1e-02

1e-01
!

1e+00

1e+01

1e-03

1e-02

1e+00

1e+01

(b) σ = 0.5.

0.008(a) σ = 0.1.

1e-01
!

0.002

Difference measure
0.004
0.006

Norm boundBridge boundTrue difference

1e-03

1e-02

1e-01
!

1e+00

1e+01

(c) σ = 1.0.

Figure 7: The approximation error of uLSIF against LSIF as functions of the regularization parameter λ for different Gaussian width σ. Its upper bounds are also plotted in the graphs.

{Xr }r6=k (i.e., without Xk ) and compute its log-likelihood for Xk :1∑
|Xk | x∈Xk

log pbXk (x).

Repeat this procedure for r = 1, 2, . . . , R and choose the value of σ such that the average of the above
hold-out log-likelihood over all r is maximized. Note that the average hold-out log-likelihood is analmost unbiased estimate of the Kullback-Leibler divergence from p(x) to pb(x), up to an irrelevant
constant.
KDE can be used for importance estimation by first obtaining density estimators pbtr (x) andtetrb, and then estimating the importance by w(x)=and {xtej }nj=1pbte (x) separately from {xtri }ni=1pbte (x)/ pbtr (x). A potential limitation of this approach is that KDE suffers from the curse of dimensionality (Vapnik, 1998; Härdle et al., 2004), that is, the number of samples needed to maintain the1414A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

same approximation quality grows exponentially as the dimension of the domain increases. This is
critical when the number of available samples is limited. Therefore, the KDE-based approach may
not be reliable in high-dimensional problems.
5.2 Kernel Mean Matching
The kernel mean matching (KMM) method allows us to directly obtain an estimate of the importance values at training points without going through density estimation (Huang et al., 2007). Theb such that the mean discrepancy between nonlinearly transformed
basic idea of KMM is to find w(x)samples drawn from pte (x) and ptr (x) is minimized in a universal reproducing kernel Hilbert space(Steinwart, 2001). The Gaussian kernel (21) is an example of kernels that induce universal reproducing kernel Hilbert spaces and it has been shown that the solution of the following optimization
problem agrees with the true importance:
min
w(x)

Z

Kσ (x, ·)pte (x)dx −subject to

Z

Z

2

Kσ (x, ·)w(x)ptr (x)dxH

w(x)ptr (x)dx = 1 and w(x) ≥ 0,

where k · kH denotes the norm in the Gaussian reproducing kernel Hilbert space and Kσ (x, x′ ) is theGaussian kernel (21).
An empirical version of the above problem is reduced to the following quadratic program:
"
#
ntr
1 ntr
tr
tr
min
wi wi′ Kσ (xi , xi′ ) − ∑ wi κi
ntr
2 i,i∑{wi }i=1
′ =1
i=1
ntr

subject to∑ wi − ntr ≤ ntr ε and 0 ≤ w1 , w2 , . . . , wn ≤ B,
tr

i=1

where
κi =

ntr nte∑ Kσ (xtri , xtej ).
nte j=1

B (≥ 0) and ε (≥ 0) are tuning parameters that control the regularization effects. The solutiontrtrbi }ni=1.
is an estimate of the importance at the training points {xtri }ni=1{w
Since KMM does not involve density estimation, it is expected to work well even in high dimensional cases. However, the performance is dependent on the tuning parameters B, ε, and σ, and theycannot be simply optimized, for example, by CV since estimates of the importance are available
only at the training points. A popular heuristic is to use the median distance between samples as
the Gaussian width σ, which is shown to be useful (Schölkopf and Smola, 2002; Song et al., 2007).
However, there seems no strong justification for this heuristic. For the choice of ε, a theoreticalresult given in Huang et al. (2007) could be used as guidance, although it is still hard to determine
the best value of ε in practice.
5.3 Logistic Regression
Another approach to directly estimating the importance is to use a probabilistic classifier. Let us
assign a selector variable η = −1 to training samples and η = 1 to test samples, that is, the training
1415

K ANAMORI , H IDO AND S UGIYAMA

and test densities are written as
ptr (x) = p(x|η = −1),pte (x) = p(x|η = 1).

Note that η is regarded as a random variable.
Application of the Bayes theorem yields that the importance can be expressed in terms of η asfollows (Qin, 1998; Cheng and Chu, 2004; Bickel et al., 2007):
w(x) =p(η = −1) p(η = 1|x).
p(η = 1) p(η = −1|x)

The probability ratio of test and training samples may be simply estimated by the ratio of the numbers of samples:p(η = −1) ntr≈.
p(η = 1)nte
The conditional probability p(η|x) could be approximated by discriminating test samples from training samples using a logistic regression (LogReg) classifier, where η plays the role of a class variable.
Below we briefly explain the LogReg method.
The LogReg classifier employs a parametric model of the following form for expressing theconditional probability p(η|x):pb(η|x) =1,
1 + exp (−η ∑mℓ=1 ζℓ φℓ (x))

where m is the number of basis functions and {φℓ (x)}mℓ=1 are fixed basis functions. The parameter ζis learned so that the negative regularized log-likelihood is minimized:bζ = argmin
ζ

"

ntr

!!

m

∑ log 1 + exp ∑ ζℓ φℓ (xtri )
ℓ=1

i=1

nte

m

j=1

ℓ=1

!!

+ ∑ log 1 + exp − ∑ ζℓ φℓ (x )
te

⊤

#

+ λζ ζ .

Since the above objective function is convex, the global optimal solution can be obtained by standard
nonlinear optimization methods such as Newton’s method, the conjugate gradient method, and the
BFGS method (Minka, 2007). Then the importance estimate is given by
ntr
b =
w(x)
exp
nte

m

!

∑ ζℓ φℓ (x) .

ℓ=1

(34)

An advantage of the LogReg method is that model selection (that is, the choice of the basis
functions {φℓ (x)}mℓ=1 as well as the regularization parameter λ) is possible by standard CV since thelearning problem involved above is a standard supervised classification problem.
1416

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

5.4 Kullback-Leibler Importance Estimation Procedure
The Kullback-Leibler importance estimation procedure (KLIEP) (Sugiyama et al., 2008a) also directly gives an estimate of the importance function without going through density estimation by
matching the two distributions in terms of the Kullback-Leibler divergence (Kullback and Leibler,1951).
Let us model the importance w(x) by the linear model (1). An estimate of the test density pte (x)b as
is given by using the model w(x)bpbte (x) = w(x)ptr (x).

In KLIEP, the parameters α are determined so that the Kullback-Leibler divergence from pte (x) topbte (x) is minimized:pte (x)
dx
b
w(x)p
D
tr (x)ZZ
pte (x)b
pte (x) log=
pte (x) log w(x)dx.
dx −ptr (x)D
D
Z

KL[pte (x)k pbte (x)] =pte (x) log(35)b
The first term is a constant, so it can be safely ignored. Since pbte (x) (= w(x)ptr (x)) is a probability
density function, it should satisfy
1=

Z

D

pbte (x)dx =

Z

D

b
w(x)p
tr (x)dx.

(36)

Then the KLIEP optimization problem is given by replacing the expectations in Eqs. (35) and (36)with empirical averages as follows:
!#
"
nte

max

{αℓ }bℓ=1

∑ log

j=1

b

∑ αℓ ϕℓ (xtej )ℓ=1b

subject to ∑ αℓ
ℓ=1

ntr

∑

i=1

!

ϕℓ (xtri )= ntr and α1 , α2 , . . . , αb ≥ 0.

This is a convex optimization problem and the global solution—which tends to be sparse (Boyd andVandenberghe, 2004)—can be obtained, for example, by simply performing gradient ascent andfeasibility satisfaction iteratively. Model selection of KLIEP is possible by LCV.
Properties of KLIEP-type algorithms have been theoretically investigated in Nguyen et al. (2008)
and Sugiyama et al. (2008b) (see also Qin, 1998; Cheng and Chu, 2004). Note that the importancemodel of KLIEP is the linear model (1), while that of LogReg is the log-linear model (34). A variantof KLIEP for log-linear models has been studied in Tsuboi et al. (2008).
5.5 Discussions
Table 1 summarizes properties of proposed and existing methods.
KDE is efficient in computation since no optimization is involved, and model selection is possible by LCV. However, KDE may suffer from the curse of dimensionality due to the difficulty ofdensity estimation in high dimensions.
1417

K ANAMORI , H IDO AND S UGIYAMA

Methods
KDE
KMM
LogReg
KLIEP
LSIF
uLSIF

Densityestimation
Necessary
Not necessary
Not necessary
Not necessary
Not necessary
Not necessaryModelselection
Available
Not available
Available
Available
Available
Available

Optimization
Analytic
Convex quadratic program
Convex non-linear
Convex non-linear
Convex quadratic program
Analytic

Out-of-sample
prediction
Possible
Not possible
Possible
Possible
Possible
Possible

Table 1: Relation between proposed and existing methods.

KMM can potentially overcome the curse of dimensionality by directly estimating the importance. However, there is no objective model selection method. Therefore, model parameters such as
the Gaussian width need to be determined by hand, which is highly unreliable unless we have strong
prior knowledge. Furthermore, the computation of KMM is rather demanding since a quadratic programming problem has to be solved.
LogReg and KLIEP also do not involve density estimation, but different from KMM, they givean estimate the entire importance function, not only the values of the importance at training points.
Therefore, the values of the importance at unseen points can be estimated by LogReg and KLIEP.
This feature is highly useful since it enables us to employ CV for model selection, which is a significant advantage over KMM. However, LogReg and KLIEP are computationally rather expensive
since non-linear optimization problems have to be solved. Note that the LogReg method is slightly
different in motivation from other methods, but has some similarity in computation and implementation, for example, the LogReg method also involves a kernel smoother.
The proposed LSIF method is qualitatively similar to LogReg and KLIEP, that is, it can avoid
density estimation, model selection is possible, and non-linear optimization is involved. LSIF isadvantageous over LogReg and KLIEP in that it is equipped with a regularization path tracking
algorithm. Thanks to this, model selection of LSIF is computationally much more efficient than
LogReg and KLIEP. However, the regularization path tracking algorithm tends to be numericallyunstable.
The proposed uLSIF method inherits good properties of existing methods such as no density
estimation involved and a build-in model selection method equipped. In addition to these preferable
properties, the solution of uLSIF can be computed in an efficient and numerically stable manner.
Furthermore, thanks to the availability of the closed-form solution of uLSIF, the LOOCV score canbe analytically computed without repeating hold-out loops, which highly contributes to reducing
the computation time in the model selection phase.
In the next section, we experimentally show that uLSIF is computationally more efficient thanexisting direct importance estimation methods, while its estimation accuracy is comparable to thebest existing methods.

6. Experiments
In this section, we compare the experimental performance of the proposed and existing methods.
1418

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

6.1 Importance Estimation
Let the dimension of the domain be d andptr (x) = N (x; (0, 0, . . . , 0)⊤ , I d ),pte (x) = N (x; (1, 0, . . . , 0)⊤ , I d ).
The task is to estimate the importance at training points:
wi = w(xtri ) =pte (xtri )ptr (xtri )for i = 1, 2, . . . , ntr .

We compare the following methods:
KDE(CV): The Gaussian kernel (21) is used, where the kernel widths of the training and test
densities are separately optimized based on 5-fold LCV.
KMM(med): The performance of KMM is dependent on B, ε, and σ. We set B = 1000 and ε =√√( ntr −1)/ ntr following the original paper (Huang et al., 2007), and the Gaussian width σ isset at the median distance between samples within the training set and the test set (Schölkopf
and Smola, 2002; Song et al., 2007).
LogReg(CV): The Gaussian kernel model (22) are used as basis functions. The kernel width σ andthe regularization parameter λ are chosen based on 5-fold CV.1
KLIEP(CV): The Gaussian kernel model (22) is used. The kernel width σ is selected based on
5-fold LCV.
uLSIF(CV): The Gaussian kernel model (22) is used. The kernel width σ and the regularization
parameter λ are determined based on LOOCV.
All the methods are implemented using the MATLAB R environment, where the CPLEX R optimizer is used for solving quadratic programs in KMM and the LIBLINEAR implementation is usedfor LogReg (Lin et al., 2007).
We fixed the number of test points at nte = 1000 and consider the following two setups for thenumber ntr of training samples and the input dimensionality d:(a) ntr is fixed at ntr = 100 and d is changed as d = 1, 2, . . . , 20,(b) d is fixed at d = 10 and ntr is changed as ntr = 50, 60, . . . , 150.
We run the experiments 100 times for each d, each ntr , and each method, and evaluate the quality oftrbi }ni=1
by the normalized mean squared error (NMSE):
the importance estimates {w
2

bi
1 ntr
wi
w
NMSE =∑ ∑n′tr wbi′ − ∑n′tr wi′ .
ntr i=1i =1i =11. In Sugiyama et al. (2008b) where KLIEP has been proposed, the performance of LogReg has been experimentallyinvestigated in the same setup. In that paper, however, LogReg was not regularized since KLIEP was not also
regularized. On the other hand, we use a regularized LogReg method and choose the regularization parameter in
addition to the Gaussian kernel width by CV here. Thanks to the regularization effect, the results of LogReg in thecurrent paper tends to be better than that reported in Sugiyama et al. (2008b).

1419

K ANAMORI , H IDO AND S UGIYAMAIn practice, the scale of the importance is not significant and the relative magnitude among wi is
important. Thus the above NMSE would be a suitable error metric for evaluating the performanceof each method.
NMSEs averaged over 100 trials (a) as a function of input dimensionality d and (b) as a functionof the training sample size ntr are plotted in log scale in Figure 8. Error bars are omitted for clear
visibility—instead, the best method in terms of the mean error and comparable ones based on thet-test at the significance level 1% are indicated by ‘◦’; the methods with significant difference fromthe best methods are indicated by ‘×’.
Figure 8(a) shows that the error of KDE(CV) sharply increases as the input dimensionality
grows, while LogReg, KLIEP, and uLSIF tend to give much smaller errors than KDE. This would
be the fruit of directly estimating the importance without going through density estimation. KMM
tends to perform poorly, which is caused by an inappropriate choice of the Gaussian kernel width.
On the other hand, model selection in LogReg, KLIEP, and uLSIF seems to work quite well. Figure 8(b) shows that the errors of all methods tend to decrease as the number of training samples
grows. Again LogReg, KLIEP, and uLSIF tend to give much smaller errors than KDE and KMM.
Next we investigate the computation time. Each method has a different model selection strategy,
that is, KMM does not involve CV, KDE and KLIEP involve CV over the kernel width, and LogReg
and uLSIF involve CV over both the kernel width and the regularization parameter. Thus the naive
comparison of the total computation time is not so meaningful. For this reason, we first investigate
the computation time of each importance estimation method after the model parameters are fixed.
The average CPU computation time over 100 trials are summarized in Figure 9. Figure 9(a)shows that the computation time of KDE, KLIEP, and uLSIF is almost independent of the input
dimensionality, while that of KMM and LogReg is rather dependent on the input dimensionality.
Note that LogReg for d ≤ 3 is slow due to some convergence problem of the LIBLINEAR package.
Among them, the proposed uLSIF is one of the fastest methods. Figure 9(b) shows that the computation time of LogReg, KLIEP, and uLSIF is nearly independent of the number of training samples,while that of KDE and KMM sharply increase as the number of training samples increases.
Both LogReg and uLSIF have high accuracy and their computation time after model selection
is comparable. Finally, we compare the entire computation time of LogReg and uLSIF including
CV, which is summarized in Figure 10. We note that the Gaussian width σ and the regularization
parameter λ are chosen over the 9×9 grid in this experiment for both LogReg and uLSIF. Therefore,
the comparison of the entire computation time is fair. Figures 10(a) and 10(b) show that uLSIF is
approximately 5 times faster than LogReg.
Overall, uLSIF is shown to be comparable to the best existing method (LogReg) in terms of theaccuracy, but is computationally more efficient than LogReg.
6.2 Covariate Shift Adaptation in Regression and Classification
Next, we illustrate how the importance estimation methods could be used in covariate shift adaptation (Shimodaira, 2000; Zadrozny, 2004; Sugiyama and Müller, 2005; Huang et al., 2007; Bickel
and Scheffer, 2007; Bickel et al., 2007; Sugiyama et al., 2007). Covariate shift is a situation insupervised learning where the input distributions change between the training and test phase but theconditional distribution of outputs given inputs remains unchanged. Under covariate shift, standard
learning techniques such as maximum likelihood estimation or cross-validation are biased—the bias
1420

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION−3

Average NMSE over 100 Trials (in Log Scale)10

KDE(CV)KMM(med)
LogReg(CV)KLIEP(CV)uLSIF(CV)

−4

10

−5

10

−6

10

5

10
d (Input Dimension)

1520(a) When input dimensionality is changed−3

Average NMSE over 100 Trials (in Log Scale)10

KDE(CV)KMM(med)
LogReg(CV)KLIEP(CV)uLSIF(CV)

−4

10

−5

10

−6

10

50

100
ntr (Number of Training Samples)

150(b) When training sample size is changed

Figure 8: NMSEs averaged over 100 trials in log scale for the artificial data set. Error bars areomitted for clear visibility. Instead, the best method in terms of the mean error andcomparable ones based on the t-test at the significance level 1% are indicated by ‘◦’; themethods with significant difference from the best methods are indicated by ‘×’.
1421

K ANAMORI , H IDO AND S UGIYAMAKDE(CV)KMM(med)
LogReg(CV)KLIEP(CV)uLSIF(CV)Average Computation TIME over 100 Trials [sec]

0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0

5

10
d (Input Dimension)

1520(a) When input dimensionality is changed

Average Computation TIME over 100 Trials [sec]0.15

KDE(CV)KMM(med)
LogReg(CV)KLIEP(CV)uLSIF(CV)

0.1

0.05

0
50

100
n (Number of Training Samples)

150

tr

(b) When training sample size is changed

Figure 9: Average computation time (after model selection) over 100 trials for the artificial data set.
1422

Average Total Computation TIME over 100 Trials [sec]A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

LogReg(CV)uLSIF(CV)

18
16
14
12
10
8
6
4
2
0

5

10
d (Input Dimension)

1520(a) When input dimensionality is changed

Average Total Computation TIME over 100 Trials [sec]12

LogReg(CV)uLSIF(CV)

10

8

6

4

2

0
50

100
ntr (Number of Training Samples)

150(b) When training sample size is changed

Figure 10: Average computation time over 100 trials for the artificial data set (including model
selection of the Gaussian width σ and the regularization parameter λ over the 9 × 9grid).
1423

K ANAMORI , H IDO AND S UGIYAMA

caused by covariate shift can be asymptotically canceled by weighting the loss function according
to the importance.
tr
In addition to training input samples {xtri }ni=1drawn from a training input density ptr (x) and testnte
teinput samples {x j } j=1 drawn from a test input density pte (x), suppose that we are given trainingtrtr. The task is to predict the outputs forat the training input points {xtri }ni=1output samples {ytri }ni=1
nte
te
tr
.
test inputs {x j } j=1 based on the input-output training samples {(xtri , ytri )}ni=1
We use the following kernel model for function learning:tfb(x; θ) = ∑ θℓ Kh (x, mℓ ),ℓ=1

tewhere Kh (x, x′ ) is the Gaussian kernel (21) and mℓ is a template point randomly chosen from {xtej }nj=1
without replacement. We set the number of kernels at t = 50. We learn the parameter θ by importance weighted regularized least-squares (IWRLS) (Evgeniou et al., 2000; Sugiyama and Müller,
2005):
"
#
2

ntr
2
tr
tr
tr
b
b i ) fb(xi ; θ) − yi + γkθk .
θIWRLS ≡ argmin ∑ w(x
(37)

θ

i=1

It is known that IWRLS is consistent when the true importance w(xtri ) is used as weights—unweighted
RLS is not consistent due to covariate shift, given that the true learning target function f (x) is notrealizable by the model fb(x) (Shimodaira, 2000).
The solution b
θIWRLS is analytically given bywherebb K + γI b )−1 K ⊤W
b ytr ,θIWRLS = (K ⊤W

Ki,ℓ = Kh (xtri , mℓ ),b = diag w(xb tr1 ), w(xb tr2 ), . . . , w(xb trn ) ,
W
tr

tr

y

= (ytr1 , ytr2 , . . . , ytrntr )⊤ .

diag (a, b, . . . , c) denotes the diagonal matrix with the diagonal elements a, b, . . . , c.
The kernel width h and the regularization parameter γ in IWRLS (37) are chosen by importanceweighted CV (IWCV) (Sugiyama et al., 2007). More specifically, we first divide the training samples
tr{ztri | ztri = (xtri , ytri )}ni=1into R disjoint subsets {Zrtr }Rr=1 . Then a function fbr (x) is learned usingtr{Z j } j6=r by IWRLS and its mean test error for the remaining samples Zrtr is computed:
1∑ tr|Zrtr | (x,y)∈
Z
r

where



b
fbr (x), y ,w(x)loss((by − y)2(Regression),loss (by, y) = 1yy}) (Classification).
2 (1 − sign{bWe repeat this procedure for r = 1, 2, . . . , R and choose the kernel width h and the regularization
parameter γ so that the average of the above mean test error over all r is minimized. We set the
number of folds in IWCV at R = 5. IWCV is shown to be an (almost) unbiased estimator of the1424

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

KDE
KMM
LogReg
KLIEP
uLSIF
(CV)
(med)
(CV)
(CV)
(CV)
kin-8fh1.00(0.34) 1.22(0.52) 1.55(0.39) 1.31(0.39) 0.95(0.31) 1.02(0.33)
kin-8fm 1.00(0.39) 1.12(0.57) 1.84(0.58) 1.38(0.57) 0.86(0.35) 0.88(0.39)
kin-8nh
1.00(0.26) 1.09(0.20) 1.19(0.29) 1.09(0.19) 0.99(0.22) 1.02(0.18)
kin-8nm 1.00(0.30) 1.14(0.26) 1.20(0.20) 1.12(0.21) 0.97(0.25) 1.04(0.25)
abalone1.00(0.50) 1.02(0.41) 0.91(0.38) 0.97(0.49) 0.94(0.67) 0.96(0.61)
image
1.00(0.51) 0.98(0.45) 1.08(0.54) 0.98(0.46) 0.94(0.44) 0.98(0.47)
ringnorm 1.00(0.04) 0.87(0.04) 0.87(0.04) 0.95(0.08) 0.99(0.06) 0.91(0.08)
twonorm 1.00(0.58) 1.16(0.71) 0.94(0.57) 0.91(0.61) 0.91(0.52) 0.88(0.57)
waveform 1.00(0.45) 1.05(0.47) 0.98(0.31) 0.93(0.32) 0.93(0.34) 0.92(0.32)
Average 1.00(0.38) 1.07(0.40) 1.17(0.37) 1.07(0.37) 0.94(0.35) 0.96(0.36)
Comp. time
—
0.82
3.50
3.27
2.23
1.00
Data

Uniform

Table 2: Mean test error averaged over 100 trials for covariate shift adaptation in regression andclassification. The numbers in the brackets are the standard deviation. All the error values
are normalized by that of ‘Uniform’ (uniform weighting, or equivalently no importanceweighting). For each data set, the best method in terms of the mean error and comparable
ones based on the Wilcoxon signed rank test at the significance level 1% are described inbold face. The upper half corresponds to regression data sets taken from DELVE (Rasmussen et al., 1996), while the lower half correspond to classification data sets taken from
IDA (Rätsch et al., 2001). All the methods are implemented using the MATLAB R environment, where the CPLEX R optimizer is used for solving quadratic programs in KMM andthe LIBLINEAR implementation is used for LogReg (Lin et al., 2007).

generalization error, while unweighted CV with misspecified models is biased due to covariate shift(Zadrozny, 2004; Sugiyama et al., 2007).
The data sets provided by DELVE (Rasmussen et al., 1996) and IDA (Rätsch et al., 2001)
are used for performance evaluation. Each data set consists of input/output samples {(xk , yk )}nk=1 .
te
We normalize all the input samples {xk }nk=1 into [0, 1]d and choose the test samples {(xtej , ytej )}nj=1nfrom the pool {(xk , yk )}k=1 as follows. We randomly choose one sample (xk , yk ) from the pool and(c)(c)
accept this with probability min(1, 4(xk )2 ), where xk is the c-th element of xk and c is randomly
determined and fixed in each trial of the experiments. Then we remove xk from the pool regardlessof its rejection or acceptance, and repeat this procedure until nte samples are accepted. We choosetrthe training samples {(xtri , ytri )}ni=1
uniformly from the rest. Thus, in this experiment, the test input(c)density tends to be lower than the training input density when xk is small. We set the number oftetrand {xtej }nj=1
samples at ntr = 100 and nte = 500 for all data sets. Note that we only use {(xtri , ytri )}ni=1ntete
for training regressors or classifiers; the test output values {y j } j=1 are used only for evaluating thegeneralization performance.
We run the experiments 100 times for each data set and evaluate the mean test error:


1 nte
te
te
b
f
(x
),
y
loss
∑
j
j .
nte j=1
1425

K ANAMORI , H IDO AND S UGIYAMA

The results are summarized in Table 2, where ‘Uniform’ denotes uniform weights (or equivalently,
no importance weight). The numbers in the brackets are the standard deviation. All the error values
are normalized so that the mean error of Uniform is one. For each data set, the best method in termsof the mean error and comparable ones based on the Wilcoxon signed rank test at the significancelevel 1% are described in bold face. The upper half of the table corresponds to regression data setstaken from DELVE (Rasmussen et al., 1996), while the lower half correspond to classification data
sets taken from IDA (Rätsch et al., 2001). All the methods are implemented using the MATLAB R
environment, where the CPLEX R optimizer is used for solving quadratic programs in KMM andthe LIBLINEAR implementation is used for LogReg (Lin et al., 2007).
The table shows that the generalization performance of uLSIF tends to be better than that ofUniform, KDE, KMM, and LogReg, while it is comparable to the best existing method (KLIEP).
The mean computation time over 100 trials is described in the bottom row of the table, where thevalue is normalized so that the computation time of uLSIF is one. This shows that the computation
time of uLSIF is much shorter than KLIEP. Thus, uLSIF is overall shown to be useful in covariate
shift adaptation.
6.3 Outlier Detection
Finally, we apply importance estimation methods in outlier detection.
Here, we consider an outlier detection problem of finding irregular samples in a data set (“evaluation data set”) based on another data set (“model data set”) that only contains regular samples.
Defining the importance over two sets of samples, we can see that the importance values for regularsamples are close to one, while those for outliers tend to be significantly deviated from one. Thus
the importance values could be used as an index of the degree of outlyingness in this scenario. Since
the evaluation data set has wider support than the model data set, we regard the evaluation data settr
as the training set {xtri }ni=1(that is, the denominator in the importance) and the model data set asntete
the test set {x j } j=1 (that is, the numerator in the importance). Then outliers tend to have smaller
importance values (that is, close to zero).
We again test KMM(med), LogReg(CV), KLIEP(CV), and uLSIF(CV) for importance estimation; in addition, we include native outlier detection methods for comparison purposes. The outlier
detection problem that the native methods used below solve is to find outliers in a single data set{xk }nk=1 —the native methods can be employed in the current scenario just by finding outliers fromall samples:tetr.
{xk }nk=1 = {xtri }ni=1∪ {xtej }nj=1

One-class support vector machine (OSVM): The support vector machine (SVM) (Vapnik, 1998;
Schölkopf and Smola, 2002) is one of the most successful classification algorithms in machinelearning. The core idea of SVM is to separate samples in different classes by the maximum
margin hyperplane in a kernel-induced feature space.
OSVM is an extension of SVM to outlier detection (Schölkopf et al., 2001). The basic ideaof OSVM is to separate data samples {xk }nk=1 into outliers and inliers by a hyperplane in a
Gaussian reproducing kernel Hilbert space. More specifically, the solution of OSVM is given1426

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

as the solution of the following convex quadratic programming problem:
1 n∑ wk wk′ Kσ (xk , xk′ ){wk }k=1 2 k,k′ =1
min
n

n

subject to ∑ wk = 1 and 0 ≤ w1 , w2 , . . . , wn ≤
k=1

1
,
νn

where ν (0 ≤ ν ≤ 1) is the maximum fraction of outliers.

We use the inverse distance of a sample from the separating hyperplane as an outlier score.
The OSVM solution is dependent on the outlier ratio ν and the Gaussian kernel width σ, andthere seems to be no systematic method to determine the values of these tuning parameters.
Here we use the median distance between samples as the Gaussian width, which is a popularheuristic (Schölkopf and Smola, 2002; Song et al., 2007). The value of ν is fixed at the true
output ratio, that is, the ideal optimal value. Thus the simulation results below should beslightly in favor of OSVM.
Local outlier factor (LOF): LOF is the score to detect a local outlier which lies relatively far fromthe nearest dense region (Breunig et al., 2000). For a prefixed natural number k, the LOF valueof a sample x is defined by
LOFR (x) =1 k imdk (nearesti (x))
,
∑
k i=1
imdk (x)where nearesti (x) denotes the i-th nearest neighbor of x and imdk (x) denotes the inverse mean
distance from x to its k nearest neighbors:
1

imdk (x) = 1kk ∑i=1 kx − nearesti (x)k.

If x alone is apart from a cloud of points, imdk (x) tends to become smaller than than imdk (nearesti (x))
for all i. Then the LOF value gets large and therefore such a point is regarded as an outlier.
The performance of LOF depends on the choice of the parameter k and there seems no systematic way to find an appropriate value of k. Here we test several different values of k.
Kernel density estimator (KDE’): A naive density estimation of all data samples {xk }nk=1 can also
be used for outlier detection since the density value itself could be regarded as an outlier score.
We use KDE with the Gaussian kernel (21) for density estimation, where the kernel width isdetermined based on 5-fold LCV.
All the methods are implemented using the R environment—we use the ksvm routine in thekernlab package for OSVM (Karatzoglou et al., 2004) and the lofactor routine in the dprep package
for LOF (Fernandez, 2005).
The data sets provided by IDA (Rätsch et al., 2001) are used for performance evaluation. These
data sets are binary classification data sets consisting of positive/negative and training/test samples.
We allocate all positive training samples for the “model” set, while all positive test samples and afraction ρ (= 0.01, 0.02, 0.05) of negative test samples are assigned in the “evaluation” set. Thus,
we regard the positive samples as regular and the negative samples as irregular.
1427

K ANAMORI , H IDO AND S UGIYAMA

In the evaluation of the performance of outlier detection methods, it is important to take into
account both the detection rate (the amount of true outliers an outlier detection algorithm can find)and the detection accuracy (the amount of true inliers that an outlier detection algorithm misjudges
as outliers). Since there is a trade-off between the detection rate and the detection accuracy, we
adopt the area under the ROC curve (AUC) as our error metric (Bradley, 1997).
The mean AUC values over 20 trials as well as the computation time are summarized in Table 3,
showing that uLSIF works fairly well. KLIEP works slightly better than uLSIF, but uLSIF is computationally much more efficient. LogReg overall works reasonably well, but it performs poorly forsome data sets and the average AUC performance is not as good as uLSIF or KLIEP. KMM andOSVM are not comparable to uLSIF in AUC and they are computationally inefficient. Note that we
also tested KMM and OSVM with several different Gaussian widths and experimentally found thatthe heuristic of using the median sample distance as the Gaussian kernel width works reasonably
well in this experiment. Thus the AUC values of KMM and OSVM are close to optimal. LOF withlarge k is shown to work well, although it is not clear whether the heuristic of simply using large k
is always appropriate or not. The computational cost of LOF is high since nearest neighbor searchis computationally expensive. KDE’ works reasonably well, but its performance is not as good as
uLSIF and KLIEP.
Overall, uLSIF is shown to work well with low computational costs.

7. Conclusions
The importance is useful in various machine learning scenarios such as covariate shift adaptationand outlier detection. In this paper, we proposed a new method of importance estimation that can
avoid solving a substantially more difficult task of density estimation. We formulated the importanceestimation problem as least-squares function fitting and casted the optimization problem as a convex
quadratic program (we referred to it as LSIF). We theoretically elucidated the convergence propertyof LSIF and showed that the entire regularization path of LSIF can be efficiently computed basedon a parametric optimization technique. We further developed an approximation algorithm (wecalled it uLSIF), which allows us to obtain the closed-form solution. We showed that the leave-oneout cross-validation score can be computed analytically for uLSIF—this makes the computation of
uLSIF highly efficient. We carried out extensive simulations in covariate shift adaptation and outlier
detection, and experimentally confirmed that the proposed uLSIF is computationally more efficientthan existing approaches, while the accuracy of uLSIF is comparable to the best existing methods.
Thanks to the low computational cost, uLSIF would be highly scalability to large data sets, which
is very important in practical applications.
We have given convergence proofs for LSIF and uLSIF. A possible future direction to pursue
along this line is to show the convergence of LSIF and uLSIF in non-parametric cases, for example,
following Nguyen et al. (2008) and Sugiyama et al. (2008b). We are currently exploring various
possible applications of important estimation methods beyond covariate shift adaptation or outlier
detection, for example, feature selection, conditional distribution estimation, independent component analysis, and dimensionality reduction—we believe that importance estimation could be used
as a new versatile tool in statistical machine learning.

1428

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

Data
Name
bananab-cancerdiabetes

f-solargermanheart

satimagesplicethyroidtitanic

twonorm

waveform

ρ
0.01
0.02
0.05
0.01
0.02
0.05
0.01
0.02
0.05
0.01
0.02
0.05
0.01
0.02
0.05
0.01
0.02
0.05
0.01
0.02
0.05
0.01
0.02
0.05
0.01
0.02
0.05
0.01
0.02
0.05
0.01
0.02
0.05
0.01
0.02
0.05

AverageComp. timeuLSIF KLIEP LogReg(CV) (CV) (CV)0.851 0.815 0.447
0.858 0.824 0.428
0.869 0.851 0.435
0.463 0.480 0.627
0.463 0.480 0.627
0.463 0.480 0.627
0.558 0.615 0.599
0.558 0.615 0.599
0.532 0.590 0.636
0.416 0.485 0.438
0.426 0.456 0.432
0.442 0.479 0.432
0.574 0.572 0.556
0.574 0.572 0.556
0.564 0.555 0.540
0.659 0.647 0.833
0.659 0.647 0.833
0.659 0.647 0.833
0.812 0.828 0.600
0.829 0.847 0.632
0.841 0.858 0.715
0.713 0.748 0.368
0.754 0.765 0.343
0.734 0.764 0.377
0.534 0.720 0.745
0.534 0.720 0.745
0.534 0.720 0.745
0.525 0.534 0.602
0.496 0.498 0.659
0.526 0.521 0.644
0.905 0.902 0.161
0.896 0.889 0.197
0.905 0.903 0.396
0.890 0.881 0.243
0.901 0.890 0.181
0.885 0.873 0.236
0.661 0.685 0.530
1.00 11.7
5.35

KMM
(med)
0.578
0.644
0.761
0.576
0.576
0.576
0.574
0.574
0.547
0.494
0.480
0.532
0.529
0.529
0.532
0.623
0.623
0.623
0.813
0.861
0.893
0.541
0.588
0.643
0.681
0.681
0.681
0.502
0.513
0.538
0.439
0.572
0.754
0.477
0.602
0.757
0.608
751

OSVM
(med)
0.360
0.412
0.467
0.508
0.506
0.498
0.563
0.563
0.545
0.522
0.550
0.576
0.535
0.535
0.530
0.681
0.678
0.681
0.540
0.548
0.536
0.737
0.744
0.723
0.504
0.505
0.485
0.456
0.526
0.505
0.846
0.821
0.781
0.861
0.817
0.798
0.596
12.4

LOF
k = 5 k = 30 k = 50
0.838 0.915 0.919
0.813 0.918 0.920
0.786 0.907 0.909
0.546 0.488 0.463
0.521 0.445 0.428
0.549 0.480 0.452
0.513 0.403 0.390
0.526 0.453 0.434
0.536 0.461 0.447
0.480 0.441 0.385
0.442 0.406 0.343
0.455 0.417 0.370
0.526 0.559 0.552
0.553 0.549 0.544
0.548 0.571 0.555
0.407 0.659 0.739
0.428 0.668 0.746
0.440 0.666 0.749
0.909 0.930 0.896
0.785 0.919 0.880
0.712 0.895 0.868
0.765 0.778 0.768
0.761 0.793 0.783
0.764 0.785 0.777
0.259 0.111 0.071
0.259 0.111 0.071
0.259 0.111 0.071
0.520 0.525 0.525
0.492 0.503 0.503
0.499 0.512 0.512
0.812 0.889 0.897
0.803 0.892 0.901
0.765 0.858 0.874
0.724 0.887 0.889
0.690 0.887 0.890
0.705 0.847 0.874
0.594 0.629 0.622
85.5

KDE’
(CV)
0.934
0.927
0.923
0.400
0.400
0.400
0.425
0.425
0.435
0.378
0.374
0.346
0.561
0.561
0.547
0.638
0.638
0.638
0.916
0.898
0.892
0.845
0.848
0.849
0.256
0.256
0.256
0.461
0.472
0.433
0.875
0.858
0.807
0.861
0.861
0.831
0.623
8.70

Table 3: Mean AUC values for outlier detection over 20 trials for the benchmark data sets. All themethods are implemented using the R environment, where quadratic programs in KMM
are solved by the ipop optimizer (Karatzoglou et al., 2004), the ksvm routine is used for
OSVM (Karatzoglou et al., 2004), and the lofactor routine is used for LOF (Fernandez,2005).

1429

K ANAMORI , H IDO AND S UGIYAMA

Acknowledgments
The authors wish to thank Issei Sato for fruitful discussion and helpful comments. The authors
would also like to thank the anonymous referees whose comments helped to improve the paper
further. This work was supported by MEXT (20680007), SCAT, and AOARD.

b
Appendix A. Existence of the Inverse Matrix of G

Here we prove Lemma 1.
Let us consider the following system of linear equations:
! 
!
⊤
b
b
0
x
H
−E
b
=
,
b Ob b0|Ab |y−E(38)|A |×|A |b |-dimensional vectors, respectively. From the upper half of Eq. (38), we
where x and y are b- and |Ahaveb −1 Eb ⊤ y.
x=H

Substituting this into the lower half of Eq. (38), we havebH
b −1 Eb ⊤y = 0 b .
E|A |b |, that is, E
b is |Ab is a row-full rank matrix. As a result,From the definition, the rank of the matrix E
−1 ⊤
bH
b E
b is invertible. Therefore, Eq. (38) has the unique solution x = 0b and y = 0 b .
the matrix E|A |b is invertible.
This implies that G

Appendix B. Active Set of LSIF

Here, we prove Theorem 2.
We prove that the active set A does not change under the infinitesimal shift of H and h if thestrict complementarity condition is satisfied. We regard the pair of a symmetric matrix and a vector(H ′ , h′ ) as an element in the ( b(b+1)+ b)-dimensional Euclidean space. We consider the following
2
linear equation:
 ′  ′−1  ′αh − λ1b
H
−E ⊤
=
,
0|A |
−E O|A |×|A |ξ′
where E is the |A | × b indicator matrix determined from the active set A (see Section 2.3 for thedetailed definition). If H ′ = H and h′ = h hold, the solution (α′ , ξ′ ) = (α∗ (λ), ξ∗ (λ)) satisfies
α′ℓ = 0, ξ′ℓ > 0, ∀ ℓ ∈ A ,α′ℓ > 0, ξ′ℓ = 0, ∀ ℓ 6∈ A ,(39)

because of the strict complementarity condition. On the other hand, if the norm of (H ′ , h′ ) − (H, h)is infinitesimal, the solution (α′ , ξ′ ) also satisfies Eq. (39) because of the continuity of the relation
between (H ′ , h′ ) and (α′ , ξ′ ).
1430

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATIONb(b+1)As a result, there exists an ε-ball Bε in R 2 +b such that the equality A = {ℓ | α′ℓ = 0} holdsb ) ≤ P((H,b bfor any (H ′ , h′ ) ∈ Bε . Therefore, we have P(A 6= Ah) 6∈ Bε ). Due to the large deviation
principle (Dembo and Zeitouni, 1998), there is a positive constant c such that−1b b
log P((H,h) 6∈ Bε ) > c > 0,min{ntr , nte }b 6= A ) < e−c min{ntr ,nte } holds.
if min{ntr , nte } is large enough. Thus, asymptotically P(AAppendix C. Learning Curve of LSIF

Here, we prove Theorem 3.
Let us consider the ideal problem (7). Let α∗ (λ) and ξ∗ (λ) be the optimal parameter and La′
grange multiplier (that is, the KKT conditions are fulfilled; see Section 2.3) and let ξ∗ (λ) be the′∗vector of non-zero elements of ξ (λ) defined in the same way as Eq. (11). Then α∗ (λ) and ξ∗ (λ)
satisfy ∗  α (λ)
h − λ1b
=
,
(40)
G ∗′
0|A |
ξ (λ)
where




H
−E ⊤
G=
.
−E O|A |×|A |From the central limit theorem and the assumption (18), we have

 

1
1
b
h = h + Op √= h + op
,
nte
ntr

(41)

where O p and o p denote the asymptotic order in probability. The assumption (a) implies that theequality
b=E
E
(42)

b is the same size as G if
holds with exponentially high probability due to Theorem 2. Note that Gb
E = E. Thus we haveb = G + δG,
G

where


δH
Ob×|A |,δG =O|A |×b O|A |×|A |b − H.
δH = H

Combining Eqs. (12), (40), (41), and (42), we have
!
 
 ∗ 
b (λ)
α
−1
1
α (λ)
b
′
+ op
= G G ∗′
.
b
ξ (λ)
ntr
ξ (λ)

(43)

(44)

The matrix Taylor expansion (Petersen and Pedersen, 2007) yieldsb −1 = G−1 − G−1 δGG−1 + G−1 δGG−1 δGG−1 − · · · ,
G
1431

(45)

K ANAMORI , H IDO AND S UGIYAMA

and the central limit theorem asserts that
δH = O p



1
√
ntr



.

(46)

Combining Eqs. (44), (45), (14), and (46), we haveb (λ) − α∗ (λ)δα = α= −AδHα∗ (λ) + AδHAδHα∗ (λ) + o



1
ntr



(47)
.

(48)

Through direct calculation, we can confirm that
AHA = A.

(49)α∗ (λ) = A(h − λ1b ).

(50)

A(Hα∗ (λ) − h) = −λA1b .

(51)

E[δH] = Ob×b .

(52)

Similar to Eq. (15), it holds that

From Eqs. (49) and (50), we have

Eqs. (43), (4), and (3) imply

From Eqs. (2) and (47), we have1b (λ)) = J(α∗ (λ)) + δα⊤ Hδα + (Hα∗ (λ) − h)⊤ δα.
J(α
2

(53)

From Eqs. (46), (48), and (49), we havehihiE δα⊤ Hδα = tr(H E δαδα⊤ )
 
h
i
1
∗
∗
⊤
= tr(AHA E (δHα (λ))(δHα (λ)) ) + o
ntr
 
h
i
1
= tr(A E (δHα∗ (λ))(δHα∗ (λ))⊤ ) + o
.
ntr

(54)

From Eqs. (48), (51), and (52), we have
 
h
i
h
i
1
⊤
∗
∗
∗
⊤
∗
E δα (Hα (λ) − h) = − E (δHα (λ) − δHAδHα (λ)) A(Hα (λ) − h) + o
ntr
 
h
i
1
∗
∗
⊤
=E (δHα (λ) − δHAδHα (λ)) λA1b + o
ntr
 
h
i
1
= − λtr(A E (δHα∗ (λ))(δHA1b )⊤ ) + o
.
(55)
ntr
1432

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

Combining Eqs. (53), (54), and (55), we haveih√√1tr(A E ( ntr δHα∗ (λ))( ntr δHα∗ (λ))⊤ )
2ntr
 
i
h√
√
λ
1
∗
⊤
− tr(A E ( ntr δHα (λ))( ntr δHA1b ) ) + o
.
ntr
ntr
√
According to the central limit theorem, ntr δHi, j asymptotically follows the normal distribution
with mean zero and varianceZ
ϕ2i (x)ϕ2j (x)ptr (x)dx − Hi,2 j ,b (λ))] =J(α∗ (λ)) +E [J(αand the asymptotic covariance betweenZ√

ntr δHi, j and√ntr δHi′ , j′ is given by

ϕi (x)ϕ j (x)ϕi′ (x)ϕ j′ (x)ptr (x)dx − Hi, j Hi′ , j′ .

Then we haveh√i√lim E ( ntr δHα∗ (λ))( ntr δHα∗ (λ))⊤ = Cw∗ ,w∗ ,ntr →∞
h√i√lim E ( ntr δHα∗ (λ))( ntr δHA1b )⊤ = Cw∗ ,v ,ntr →∞

where Cw,w′ is the b × b covariance matrix with the (ℓ, ℓ′ )-th element being the covariance between
w(x)ϕℓ (x) and w′ (x)ϕℓ′ (x) under ptr (x). Then we obtain Eq. (19).

Appendix D. Regularization Path of LSIF
Here, we derive the regularization path tracking algorithm given in Figure 1.
When λ is greater than or equal to maxk b
hk , the solution of the KKT conditions (9)–(10) is
hk , and the correprovided as α = 0b , ξ = λ1b − bh ≥ 0b . Therefore, the initial value of λ0 is maxk bb (λ0 ) = 0b .
sponding optimal solution is α′
Since bξ (λ) corresponds to non-zero elements of b
ξ(λ) as shown in Eq. (11), we have(bξ′i (λ) if j = b
ji ,bξ j (λ) =0otherwise.

(56)b (λ) and b
When λ is decreased, the solutions αξ(λ) still satisfy Eqs. (12) and (56) as long as thebactive set A remains unchanged. Change points of the active set can be found by examining the nonb (λ) and b
negativity conditions of αξ(λ) as follows. Suppose λ is decreased and the non-negativity
condition
!
b (λ)
α
≥ 02b
b
ξ(λ)

b (λ′ ) ≥ 0b and bb (λ′ − ε) ≥ 0b or
is violated at λ = λ′ . That is, both αξ(λ′ ) ≥ 0b hold, and either αbb , j should be added to the active set
b j (λ′ ) = 0 for j 6∈ A
ξ(λ′ − ε) ≥ 0b is violated for any ε > 0. If α
1433

K ANAMORI , H IDO AND S UGIYAMAb ; on the other hand, if bb, αb j (λ′ ) will take a positive value and thereforeAξ j (λ′ ) = 0 for some j ∈ A
b . Then, for the updated active set, we compute thej should be removed from the active set A

solutions by Eqs. (12) and (56). Iterating this procedure until λ reaches zero, we can obtain theentire regularization path.
Note that we omitted some minor exceptional cases for the sake of simplicity—treatments for all
possible exceptions and the rigorous convergence property are exhaustively studied in Best (1982).

Appendix E. Negative Index Set of β◦ (λ)
Here we prove Theorem 4.
As explained in Appendix B, we regard the pair of a symmetric matrix and a vector (H ′ , h′ ) as+ b)-dimensional Euclidean space.
an element in the ( b(b+1)
2
We consider the linear equationβ′ = (H ′ + λI b )−1 h′ .
Due to the assumption, for H ′ = H and h′ = h, we haveβ′ℓ 6= 0, ℓ = 1, 2, . . . , b.

(57)

On the other hand, if the norm of (H ′ , h′ ) − (H, h) is infinitesimal, the solution β′ also satisfies
Eq. (57), and the sign of β′ℓ is same as that of βℓ for ℓ = 1, 2, . . . , b, because of the continuity of therelation between (H ′ , h′ ) and β′ .
b(b+1)
e holds for any
As a result, there exists an ε-ball Bε in R 2 +b such that the equality B = B
′ ′
e
b b
(H , h ) ∈ Bε . Therefore, we have P(B 6= B ) ≤ P((H,h) 6∈ Bε ). Due to the large deviation principle(Dembo and Zeitouni, 1998), there is a positive constant c such that−1b b
log P((H,h) 6∈ Bε ) > c > 0,min{ntr , nte }e ) < e−c min{ntr ,nte } holds.
if min{ntr , nte } is large enough. Thus, asymptotically P(B 6= B

Appendix F. Learning Curve of uLSIF
Here, we prove Theorem 5.
Letbλ = Hb + λI b .
B

The matrix Taylor expansion (Petersen and Pedersen, 2007) yields
−1
−1
−1
−1
−1
−1
b−1
B
λ = Bλ − Bλ δHBλ + Bλ δHBλ δHBλ − · · · .

e ⊂ {1, 2, . . . , b} be the set of negative indices of eLet Bβ(λ), that is,
e = {ℓ | eBβℓ (λ) < 0, ℓ = 1, 2, . . . , b}.

b be the b-dimensional diagonal matrix with the ℓ-th diagonal element
Let D
(
e
b ℓ,ℓ = 0 ℓ ∈ B ,D
1 otherwise.
1434

(58)

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

The assumption (a) implies that the equality
b =D
D

(59)

holds with exponentially high probability due to Theorem 4. Combining Eqs. (59), (41), (58), and(24), we have
δβ = bβ(λ) − β∗ (λ)
−1

bB
bλ b
=D
h − DB−1
λ h

◦
◦
−1
−1
= −DB−1λ δHβ (λ) + DBλ δHBλ δHβ (λ) + o



1
ntr



.

From Eqs. (46) and (60), we have
 
h
i
h
i
1
⊤
◦
◦
−1
−1
⊤
.
E δβ Hδβ = tr(Bλ DHDBλ E (δHβ (λ))(δHβ (λ)) ) + o
ntr

From Eqs. (52) and (24), we have
h
i
h
i
◦
∗
⊤ −1
E δβ⊤ (Hβ∗ (λ) − h) =E (−δHβ◦ (λ) + δHB−1
δHβ
(λ))
B
D(Hβ
(λ)
−
h)
λ
λ
 
1
+o
ntr
h
i
◦
∗
−1
⊤
=E tr(B−1(δHβ(λ))(δHB
D(Hβ
(λ)
−
h))
)
λ
λ
 
1
+o
.
ntr
Combining Eqs. (53), (61), and (62), we have
h
i
√
√
1
E J(b
β(λ)) =J(β∗ (λ)) +tr(B−1
DHDB−1E[( ntr δHβ◦ (λ))( ntr δHβ◦ (λ))⊤ )
λ
λ
2ntr
 
√
√
1
1
◦
∗
−1
−1
⊤
+ tr(Bλ E[( ntr δHβ (λ))( ntr δHBλ D(Hβ (λ) − h)) ) + o
.
ntr
ntr
According to the central limit theorem, we have√√lim E[( ntr δHβ◦ (λ))( ntr δHβ◦ (λ))⊤ ] = Cw◦ ,w◦ ,ntr →∞√√∗⊤◦
lim E[( ntr δHβ◦ (λ))( ntr δHB−1
λ D(Hβ (λ) − h)) ] = Cw ,u .
ntr →∞

Then we obtain Eq. (25).

Appendix G. ‘Norm’ Upper Bound of Approximation Error for uLSIF
Here we prove Theorem 6.
Using the weighted norm (27), we can express diff(λ) asdiff(λ) =b (λ′ ) − b
infλ′ ≥0 kαβ(λ)kHb.
ntr
b tr ; b
β(λ))
∑ w(x
i=1

1435

i

(60)

(61)

(62)

K ANAMORI , H IDO AND S UGIYAMAb (λ′ ) = 0b holds for some large λ′ . Then we immediately have
As shown in Appendix D, αdiff(λ) ≤kb
β(λ)kHb,ntr
b tr ; b
β(λ))
∑ w(x
i

i=1

b Then kb
which proves Eq. (28). Let κmax be the largest eigenvalue of H.
β(λ)kHb can be upper bounded
as
√
√
kb
β(λ)kHb ≤ κmax kb
β(λ)k2 ≤ κmax ke
β(λ)k2 ,b and the second inequality iswhere the first inequality may be confirmed by eigen-decomposing Hbeb Then an upper
clear from the definitions of β(λ) and β(λ). Let κmin be the smallest eigenvalue of H.
bound of ke
β(λ)k22 is given as⊤b + λI b )−2bke
β(λ)k22 = bh (Hh≤where the last inequality follows from κmin > 0 .
Now we have11kb
hk22 ≤ 2 kb
hk22 ,(κmin + λ)2
λ√kb
β(λ)kHb
1
κmax kb
hk2
≤
ntr
ntr
tr
tr
b
b
λ
∑ w(x ; β(λ))∑ w(x ; β(λ))
i=1

i

i

i=1

√
κmax kb
hk2
.
= n
b
tr
tr
b
b
b
∑i=1 ∑ℓ=1 ϕℓ (xi )βℓ (λ)/kβ(λ)k1 λkβ(λ)k1
1

For the denominator of the above expression, we have
ntr

b

b
βℓ (λ)ntrbbβℓ (λ)ntr( ∑ ϕℓ (xtri )) · ∑= min ∑ ϕℓ (xtri ),∑ ∑ ϕℓ (xtri ) kbβ(λ)k ≥ min
ℓ
ℓ
kb
β(λ)k

i=1 ℓ=1

1

′

′

ℓ=1

i=1

i=1

1

where the last equality follows from the non-negativity of bβℓ (λ). The reciprocal of kb
hk2 /kbβ(λ)k1 islower bounded as follows:
(
)
)
(
e
e
e
β(λ)
βℓ (λ)kb
β(λ)k1
β(λ)
= max
,0
,0
,
≥ max
= max
b
b
b
ℓ kb
khk2
khk2
khk2
hk2
∞

1

where the last equality follows from the fact that there is an ℓ such that eβℓ (λ) > 0; otherwise, we
ntr
trbhave ∑i=1 w(xi ; β) = 0 which contradicts to the assumption of the theorem. Let us put
κe =

e
β(λ)
,
kb
hk2

where κ > 0 and e ∈ Rb such that kek2 = 1. Then we have(κmax + λ)−1 ≤ κ and e⊤bh > 0.
1436

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

Note that there exists an ℓ such that eℓ > 0. Then, we have
max
ℓ

e
1
βℓ (λ)= max κeℓ = κ max eℓ ≥max eℓbℓℓκ
max + λ ℓ
khk2
≥

1

κmax + λNow we prove the following lemma.

min{max eℓ | e⊤ e = 1, e⊤bh/kbhk1 > 0}.
eℓLemma 8 Let p1 , p2 , . . . , pb (b ≥ 2) be positive numbers such thatb

∑ pℓ = 1,ℓ=1

and let1pℓ
ε = √ min.
b ℓ 1 − pℓThen, there exists no e = (e1 , e2 , . . . , eb ) ∈ Rb such that the three conditions,
b

b

ℓ=1

ℓ=1

∑ e2ℓ = 1, ∑ pℓ eℓ > 0, and eℓ < ε for ℓ = 1, 2, . . . , b

are satisfied at the same time.
Proof We suppose that e ∈ Rb satisfies the three conditions. If minℓ pℓ /(1 − pℓ ) > 1, we havepℓ > 1/2 for all ℓ. However, this is contradictory to ∑bℓ=1 pℓ = 1. Therefore, we havemin pℓ /(1 − pℓ ) ≤ 1,ℓfrom which we have√ε ≤ 1/ b.

√
The equality constraint ∑bℓ=1 e2ℓ = 1 implies√the condition that there exists an ei such that |e√i | ≥ 1/ b.
Moreover, we have e1 , e2 , . . . , eb < ε ≤ 1/ b, and thus there is an ei such that ei ≤ −1/ b. Hence,
we have
1
1
pi
pk
pk
p
√i ≤ − pi ei < ∑ pℓ eℓ < ∑ pℓ √ min= (1 − pi ) √ min
≤ √ .
k
k
1
−
p
1
−
p
b
b
b
b
k
k
ℓ6=i
ℓ6=i
This results in contradiction.
Let pℓ = b
hℓ /kbhk1 and we use Lemma 8. Note that any element of b
h is positive. Then, we havepℓkb
β(λ)k1
1
1.
≥· √ minκmax + λb ℓ ∑i6=ℓ pikbhk2

Moreover, we haveteteϕℓ (xtej )minℓ ∑nj=1
ϕℓ (xtej )minℓ ∑nj=1
minℓ b
hℓ
= b
,
≥ b
≥
min
nte
ℓ ∑i6=ℓ pinte b
ϕℓ′ (xtej )∑ℓ′ =1 ∑ j=1hℓ′∑ℓ′ =1 bpℓ1437

K ANAMORI , H IDO AND S UGIYAMA

where the last inequality follows from the assumption 0 < ϕℓ (x) ≤ 1. Therefore, we have the inequality√κmax kb
hk2
1
λ
β(λ))
∑ni=1 w(xtri ; b


p
κmax
1
nte
·
.
≤ b b κmax 1 +
ntr
nte
tr
λ
minℓ ∑i=1 ϕℓ (xi ) minℓ′ ∑ j=1 ϕℓ′ (xtej )(63)

An upper bound of κmax is given as follows. For all a ∈ Rb , the inequality
b

b

b

ℓ=1

ℓ=1

ℓ=1

− ∑ |aℓ |ϕℓ (x) ≤ ∑ aℓ ϕℓ (x) ≤ ∑ |aℓ |ϕℓ (x)(64)

holds because of the positivity of ϕℓ (x). Let us define ā ∈ Rb for given a ∈ Rb asā = (|a1 |, |a2 |, . . . , |ab |)⊤ .

Note that kāk2 = kak2 holds. Then, using Eq. (64), we obtain the inequality
!2
!2
b
b
1 ntr
1 ntr
⊤b
tr
tr
a Ha =∑ ∑ aℓ ϕℓ (xi ) ≤ ntr ∑ ∑ |aℓ |ϕℓ (xi ) = ā⊤ Hb ā,ntr i=1i=1 ℓ=1ℓ=1

for any a ∈ Rb . Therefore, we obtainb ≤ max ā⊤ H
b ā =
max a⊤ Ha
kak2 =1

kak2 =1

max

kak2 =1, a≥0bwhere the last equality is derived from the relation,ba⊤ Ha,(65){ā | kak2 = 1, a ∈ Rb } = {a | kak2 = 1, a ≥ 0b , a ∈ Rb }.

On the other hand, due to the additional constraint a ≥ 0b , the inequalityb ≤ max a⊤ Hab
max a⊤ Ha
kak2 =1, a≥0b(66)

kak2 =1

holds. From Eqs. (65) and (66), we have

1 ntr
b = max a⊤ Hab =κmax = max a⊤ Ha
max∑
kak2 =1
kak2 =1,a≥0bkak2 =1, a≥0b ntr i=1

Using the assumption 0 < ϕℓ (x) ≤ 1, we have
1 ntr
κmax =
max∑
kak2 =1, a≥0b ntr i=1
b

=

max

kak2 =1, a≥0b

∑ aℓ

ℓ=1

b

∑

ℓ=1
!2

!2

aℓ ϕℓ (xtri )!2b∑ aℓ ϕℓ (xtri )

ℓ=1

1 ntr
≤
max
∑
kak2 =1, a≥0b ntr i=1

b

∑ aℓ

ℓ=1

.

!2

b

≤

max

kak2 =1, a≥0b

= b,b · ∑ a2ℓ
ℓ=1

(67)

where the Schwarz inequality for a and 1b is used in the last inequality. The inequalities (63) and(67) lead to the inequality (29).
It is clear that the upper bound (29) is a decreasing function of λ (> 0). For the Gaussian basis
function, ϕℓ (x) is an increasing function with respect to the Gaussian width σ. Thus, Eq. (29) is adecreasing function of σ.
1438

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

Appendix H. ‘Bridge’ Upper Bound of Approximation Error for uLSIF
Here we prove Theorem 7.
From the triangle inequality, we obtaindiff(λ) ≤b (λ′ ) − bγ(λ)kHb + kbγ(λ) − b
infλ′ ≥0 kαβ(λ)kHb.
ntrtrbb ; β(λ))
∑ w(x
i=1

(68)

i

We derive an upper bound of the first term.
First, we show that the LSIF optimization problem (6) is equivalently expressed as


⊤
1 ⊤b
α Hα − bh αminα∈Rb 2
subject to α ≥ 0b , 1⊤b α ≤ c,which we refer to as LSIF′ . The KKT conditions of LSIF (6) are given as

b −b
Hα
h + λ1b − µ = 0b ,α ≥ 0b , µ ≥ 0b , α⊤ µ = 0,

where µ is the Lagrange multiplier vector. Similarly, the KKT conditions of LSIF′ are given asb −b
h + µ0 1b − µ = 0b ,
 Hα
(69)
α ≥ 0 , µ ≥ 0b , α⊤ µ = 0, ⊤ b1b α − c ≤ 0, µ0 ≥ 0, (1⊤
α
−
c)µ
=
0,
0
b

b (λ),bwhere µ and µ0 are the Lagrange multipliers. Let (αµ(λ)) be the solution of the KKT conditionsb (λ).
bof LSIF. Then, we find that (α, µ, µ0 ) = (α(λ),bµ(λ), λ) is the solution of Eq. (69) with c = 1⊤bα′b (λ) is the unique optimal
Note that LSIF is a strictly convex optimization problem, and thus αb ,bsolution. Conversely, when the solution of Eq. (69) is provided as (αµ, µ0 ), LSIF with λ = µ0 hasb.
the same optimal solution αWhen the optimal solution of LSIFq is bγ(λ), the KKT conditions of LSIFq (30) are given asb bγ(λ) − b
b = 0b ,H
h + λbγ(λ) − ηb ≥ 0b , bγ(λ) η = 0,bγ(λ) ≥ 0b , η
⊤b

(70)

(71)

b is the Lagrange multiplier vector.
where ηb (λ1 ) be the optimal solution of LSIF′ with c = 1⊤b (λ1 )Let αγ(λ), and suppose that the solution αbbcoincides with that of LSIF with λ = λ1 . Then, from Eq. (69), we havebαb (λ1 ) − b
Hh + λ1 1b − bµ(λ1 ) = 0b ,(72)

b (λ1 ) ≥ 0b , bb (λ1 )⊤b
αµ(λ1 ) ≥ 0b , αµ(λ1 ) = 0,b (λ1 ) − 1⊤1⊤γ(λ) ≤ 0,bαbb

λ1 ≥ 0,From Eqs. (70) and (72), we obtainb (λ1 ) − 1⊤(1⊤γ(λ))λ1 = 0.
bα
bb

b α
b.
b (λ1 ) − bγ(λ)) = −λ1 1b + λbγ(λ) + bH(µ(λ1 ) − η
1439

(73)
(74)

(75)

K ANAMORI , H IDO AND S UGIYAMA

Applying Eqs. (71), (73), (74), and (75), we haveb αb (λ′ ) − bγ(λ)k2Hb ≤ (αb (λ1 ) − bγ(λ))⊤ H(b (λ1 ) − bγ(λ))
inf kαλ′ ≥0b (λ1 ) − bγ(λ))⊤ 1b + λ(αb (λ1 ) − bγ(λ))⊤bγ(λ)= − λ1 (αb)
b (λ1 ) − bγ(λ))⊤ (b+ (αµ(λ1 ) − ηb − bγ(λ)⊤bb (λ1 )⊤bγ(λ) − kbγ(λ)k22 ) − αb (λ1 )⊤ η
= λ(αµ(λ1 )b (λ1 )⊤bγ(λ) − kbγ(λ)k22 ).
≤ λ(α

(76)

b (λ1 ) ≥ 0b , bγ(λ) ≥ 0b , and 1⊤
b (λ1 ) ≤ 1⊤
From αγ(λ), we havebαbbb (λ1 ) ≤ 1⊤b (λ1 )k1 = 1⊤γ(λ) ≤ kbγ(λ)k1 .
kα
bb
bα

Then we have the following inequality:b (λ1 )⊤bγ(λ) ≤ αb (λ1 )⊤ (kbγ(λ)k∞ 1b )αb (λ1 )k1 · kbγ(λ)k∞ ≤ kbγ(λ)k1 · kbγ(λ)k∞ .
= kα

(77)

For p and q such that 1/p + 1/q = 1 and 1 ≤ p, q ≤ ∞, Hölder’s inequality states thatkα ∗ βk1 ≤ kαk p · kβkq ,where α ∗ β denotes the element-wise product of α and β. Setting p = 1, q = ∞, and α = β = bγ(λ)
in Hölder’s inequality, we havekbγ(λ)k1 · kbγ(λ)k∞ − kbγ(λ)k22 ≥ 0.

(78)

Combining Eqs. (68), (76), (77), and (78), we obtainqλ(kbγ(λ)k1 · kbγ(λ)k∞ − kbγ(λ)k22 ) + kbγ(λ) − bβ(λ)kHb.
diff(λ) ≤ntr
b tri ; b
w(x
β(λ))
∑i=1

Appendix I. Closed Form of LOOCV Score for uLSIF

Here we derive a closed form expression of the LOOCV score for uLSIF (see Figure 2 for thepseudo code).
Letϕ(x) = (ϕ1 (x), ϕ2 (x), . . . , ϕb (x))⊤ .
b and the vector b
Then the matrix H
h are expressed asntr

b = 1 ∑ ϕ(xtri )ϕ(xtri )⊤ ,H
ntr i=1nte

1bh=∑ ϕ(xtej ),nte j=1
1440

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

and the coefficients eβ(λ) can be computed byebb−1β(λ) = B
λ h.

(i)

Let b
β be the estimator obtained without the i-th training sample xtri and the i-th test sample xte
i .
Then the estimator has the following closed form:
(i)
(i)
b
β (λ) = max(0b , eβ (λ)),
−1

(i)
1
1
tr
tr ⊤
e
b
(ntr H − ϕ(xi )ϕ(xi ) ) + λI b
(nteb
h − ϕ(xtej )).
β (λ) =
ntr − 1nte − 1b=H
b + λ(ntr −1) I b and e
b−1b
Let B
β=B
h in the following calculation. Using the Sherman-Woodburyntr(i)Morrison formula (33), we can simplify the expression of eβ (λ) as follows:

−1 

(i)
ntr − 1 b 1
1
nte b
tr
tr ⊤
te
e
β (λ) =h−ϕ(xi )
B − ϕ(xi )ϕ(xi )ntr
ntr
nte − 1nte − 1ntr − 1 b−1
1−1
trtr ⊤ b −1b=B +B ϕ(xi )ϕ(xi ) B
ntr
b−1 ϕ(xtr )ntr − ϕ(xtri )⊤ B
i


nte b
1
te
×
h−
ϕ(xi )
nte − 1nte − 1!
βϕ(xtri )⊤ e(ntr − 1)nte e−1trb ϕ(xi )β+B=ntr (nte − 1)b−1 ϕ(xtr )ntr − ϕ(xtri )⊤ B
i
!
−1
tr
⊤
te )
b
ϕ(x
)
B
ϕ(x
(ntr − 1)−1−1iib ϕ(xteb ϕ(xtri ) .
−
B
B
i )+
−1
tr
tr
ntr (nte − 1)⊤bntr − ϕ(xi ) B ϕ(xi )(i)b Applying
Thus the matrix inversion required for computing eβ (λ) for all i = 1, 2, . . . . , ntr is only B.
this to Eq. (32) and rearrange the formula, we can compute the LOOCV score analytically.

References
H. Akaike. A new look at the statistical model identification. IEEE Transactions on AutomaticControl, AC-19(6):716–723, 1974.
S. Amari, N. Fujita, and S. Shinomoto. Four types of learning curves. Neural Computation, 4(4):605–618, 1992.
P. Baldi and S. Brunak. Bioinformatics: The Machine Learning Approach. MIT Press, Cambridge,
1998.
D. Bertsekas, A. Nedic, and A. Ozdaglar. Convex Analysis and Optimization. Athena Scientific,
Belmont, MA, 2003.
M. J. Best. An algorithm for the solution of the parametric quadratic programming problem. CORR
Report 82-24, Faculty of Mathematics, University of Waterloo, 1982.
1441

K ANAMORI , H IDO AND S UGIYAMA

S. Bickel and T. Scheffer. Dirichlet-enhanced spam filtering based on biased samples. In
B. Schölkopf, J. Platt, and T. Hoffman, editors, Advances in Neural Information Processing Systems 19. MIT Press, Cambridge, MA, 2007.
S. Bickel, M. Brückner, and T. Scheffer. Discriminative learning for differing training and test
distributions. In Proceedings of the 24th International Conference on Machine Learning, 2007.
K. M. Borgwardt, A. Gretton, M. J. Rasch, H.-P. Kriegel, B. Schölkopf, and A. J. Smola. Integrating
structured biological data by kernel maximum mean discrepancy. Bioinformatics, 22(14):e49–
e57, 2006.
S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, Cambridge,
2004.
A. P. Bradley. The use of the area under the ROC curve in the evaluation of machine learningalgorithms. Pattern Recognition, 30(7):1145–1159, 1997.
M. M. Breunig, H.-P. Kriegel, R. T. Ng, and J. Sander. LOF: Identifying density-based local outliers.
In Proceedings of the ACM SIGMOD International Conference on Management of Data, 2000.
G. C. Cawley and N. L. C. Talbot. Fast exact leave-one-out cross-validation of sparse least-squares
support vector machines. Neural Networks, 17(10):1467–75, 2004.
S. S. Chen, D. L. Donoho, and M. A. Saunders. Atomic decomposition by basis pursuit. SIAM
Journal on Scientific Computing, 20(1):33–61, 1998.
K. F. Cheng and C. K. Chu. Semiparametric density estimation under a two-sample density ratio
model. Bernoulli, 10(4):583–604, 2004.
A. Dembo and O. Zeitouni. Large Deviations Techniques and Applications. Springer-Verlag, New
York, 1998.
B. Efron, T. Hastie, R. Tibshirani, and I. Johnstone. Least angle regression. The Annals of Statistics,
32(2):407–499, 2004.
T. Evgeniou, M. Pontil, and T. Poggio. Regularization networks and support vector machines.
Advances in Computational Mathematics, 13(1):1–50, 2000.
E. A. Fernandez. The dprep Package, 2005. URL http://math.uprm.edu/ edgar/dprep.pdf.
G. S. Fishman. Monte Carlo: Concepts, Algorithms, and Applications. Springer-Verlag, Berlin,
1996.
G. H. Golub and C. F. Van Loan. Matrix Computations. Johns Hopkins University Press, Baltimore,
MD, 1996.
H. Hachiya, T. Akiyama, M. Sugiyama, and J. Peters. Adaptive importance sampling with automatic model selection in value function approximation. In Proceedings of the Twenty-Third AAAI
Conference on Artificial Intelligence (AAAI2008), pages 1351–1356, Chicago, USA, Jul. 13–17
2008. The AAAI Press.
1442

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

L. K. Hansen and J. Larsen. Linear unlearning for crossvalidation. Advances in Computational
Mathematics, 5:269–280, 1996.
W. Härdle, M. Müller, S. Sperlich, and A. Werwatz. Nonparametric and Semiparametric Models.
Springer Series in Statistics. Springer, Berlin, 2004.
T. Hastie, S. Rosset, R. Tibshirani, and J. ZHu. The entire regularization path for the support vector
machine. Journal of Machine Learning Research, 5:1391–1415, 2004.
J. J. Heckman. Sample selection bias as a specification error. Econometrica, 47(1):153–162, 1979.
S. Hido, Y. Tsuboi, H. Kashima, M. Sugiyama, and T. Kanamori. Inlier-based outlier detection via
direct density ratio estimation. In Proceedings of IEEE International Conference on Data Mining(ICDM2008), Pisa, Italy, Dec. 15–19 2008.
V. Hodge and J. Austin. A survey of outlier detection methodologies. Artificial Intelligence Review,22(2):85–126, 2004.
A. E. Hoerl and R. W. Kennard. Ridge regression: Biased estimation for nonorthogonal problems.
Technometrics, 12(3):55–67, 1970.
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, and T. Hoffman, editors, Advances in Neural
Information Processing Systems 19, pages 601–608. MIT Press, Cambridge, MA, 2007.
A. Karatzoglou, A. Smola, K. Hornik, and A. Zeileis. An S4 package for kernel methods in R.
Journal of Statistical Planning and Inference, 11(9):1–20, 2004.
S. Konishi and G. Kitagawa. Generalized information criteria in model selection. Biometrika, 83:
875–890, 1996.
S. Kullback and R. A. Leibler. On information and sufficiency. Annals of Mathematical Statistics,
22:79–86, 1951.
C.-J. Lin, R. C. Weng, and S. S. Keerthi. Trust region Newton method for large-scale logistic
regression. Technical report, Department of Computer Science, National Taiwan University,
2007. URL http://www.csie.ntu.edu.tw/˜cjlin/liblinear/.
A. Luntz and V. Brailovsky. On estimation of characters obtained in statistical procedure of recognition. Technicheskaya Kibernetica, 3, 1969. in Russian.
T.

P. Minka.
A comparison of numerical optimizers for logistic regression.
Technical
report,Microsoft
Research,
2007.
URL
http://research.microsoft.com/˜minka/papers/logreg/minka-logreg.pdf.

J. E. Moody. The effective number of parameters: An analysis of generalization and regularization in nonlinear learning systems. In J. E. Moody, S. J. Hanson, and R. P. Lippmann, editors,
Advances in Neural Information Processing Systems 4, pages 847–854. Morgan Kaufmann Publishers, Inc., 1992.
1443

K ANAMORI , H IDO AND S UGIYAMA

X. Nguyen, M. J. Wainwright, and M. I. Jordan. Estimating divergence functions and the likelihood
ratio by penalized convex risk minimization. In Advances in Neural Information Processing
Systems 20, Cambridge, MA, 2008. MIT Press.
K. B. Petersen and M. S. Pedersen. The matrix cookbook. Technical report, Technical Universityof Denmark, 2007. URL http://www2.imm.dtu.dk/pubdb/p.php?3274.
J. Qin. Inferences for case-control and semiparametric two-sample density ratio models. Biometrika,
85(3):619–639, 1998.
J. Quiñonero-Candela, M. Sugiyama, A. Schwaighofer, and N. Lawrence, editors. Dataset Shift in
Machine Learning. MIT Press, Cambridge, MA, 2008.
C. E. Rasmussen, R. M. Neal, G. E. Hinton, D. van Camp, M. Revow, Z. Ghahramani, R. Kustra,
and R. Tibshirani. The DELVE manual, 1996. URL http://www.cs.toronto.edu/˜delve/.
G. Rätsch, T. Onoda, and K.-R. Müller. Soft margins for adaboost. Machine Learning, 42(3):287–320, 2001.
K. Scheinberg. An efficient implementation of an active set method for SVMs. Journal of Machine
Learning Research, 7:2237–2257, 2006.
B. Schölkopf and A. J. Smola. Learning with Kernels. MIT Press, Cambridge, MA, 2002.
B. Schölkopf, J. C. Platt, J. Shawe-Taylor, A. J. Smola, and R. C. Williamson. Estimating the
support of a high-dimensional distribution. Neural Computation, 13(7):1443–1471, 2001.
H. Shimodaira. Improving predictive inference under covariate shift by weighting the log-likelihood
function. Journal of Statistical Planning and Inference, 90(2):227–244, 2000.
L. Song, A. Smola, A. Gretton, K. M. Borgwardt, and J. Bedo. Supervised feature selection via dependence estimation. In Proceedings of the 24th International Conference on Machine learning,
pages 823–830, New York, NY, USA, 2007. ACM.
I. Steinwart. On the influence of the kernel on the consistency of support vector machines. Journalof Machine Learning Research, 2:67–93, 2001.
M. Stone. Cross-validatory choice and assessment of statistical predictors. Journal of the Royal
Statistical Society B, 32(2):111–147, 1974.
M. Sugiyama and K.-R. Müller. Input-dependent estimation of generalization error under covariate
shift. Statistics & Decisions, 23(4):249–279, 2005.
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.
M. Sugiyama, S. Nakajima, H. Kashima, P. von Bünau, and M. Kawanabe. Direct importanceestimation with model selection and its application to covariate shift adaptation. In J. C. Platt,
D. Koller, Y. Singer, and S. Roweis, editors, Advances in Neural Information Processing Systems
20, pages 1433–1440, Cambridge, MA, 2008a. MIT Press.
1444

A L EAST- SQUARES A PPROACH TO D IRECT I MPORTANCE E STIMATION

M. Sugiyama, T. Suzuki, S. Nakajima, H. Kashima, P. von Bünau, and M. Kawanabe. Direct importance estimation for covariate shift adaptation. Annals of the Institute of Statistical Mathematics,
60(4):699–746, 2008b.
R. S. Sutton and G. A. Barto. Reinforcement Learning: An Introduction. MIT Press, Cambridge,
MA, 1998.
D. M. J. Tax and R. P. W. Duin. Support vector data description. Machine Learning, 54(1):45–66,
2004.
R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal StatisticalSociety, Series B, 58(1):267–288, 1996.
Y. Tsuboi, H. Kashima, S. Hido, S. Bickel, and M. Sugiyama. Direct density ratio estimation forlarge-scale covariate shift adaptation. In Proceedings of 2008 SIAM International Conference on
Data Mining (SDM2008), Atlanta, Georgia, USA, 2008.
V. N. Vapnik. Statistical Learning Theory. Wiley, New York, 1998.
G. Wahba. Spline Model for Observational Data. Society for Industrial and Applied Mathematics,
Philadelphia and Pennsylvania, 1990.
P. M. Williams. Bayesian regularization and pruning using a Laplace prior. Neural Computation, 7(1):117–143, 1995.
J. R. Wolpaw, N. Birbaumer, D. J. McFarland, G. Pfurtscheller, and T. M. Vaughan. Brain-computerinterfaces for communication and control. Clinical Neurophysiology, 113(6):767–791, 2002.
B. Zadrozny. Learning and evaluating classifiers under sample selection bias. In Proceedings of theTwenty-First International Conference on Machine Learning, New York, NY, 2004. ACM Press.

1445