Extreme Dimension Reduction for Handling Covariate Shift

arXiv:1711.10938v2 [cs.LG] 12 Mar 2018

Fulton Wang
MIT
fultonw@mit.edu

Abstract
In the covariate shift learning scenario, the training and test covariate distributions differ, so that
a predictor’s average loss over the training and
test distributions also differ. In this work, we
explore the potential of extreme dimension reduction, i.e. to very low dimensions, in improving
the performance of importance weighting methods for handling covariate shift, which fail in high
dimensions due to potentially high train/test covariate divergence and the inability to accurately
estimate the requisite density ratios. We first formulate and solve a problem optimizing over linear
subspaces a combination of their predictive utility
and train/test divergence within. Applying it to
simulated and real data, we show extreme dimension reduction helps sometimes but not always,
due to a bias introduced by dimension reduction.

1

Introduction

Often, the population for which one has labelled training
data differs from the population one needs to make predictions for. In handling this discrepancy one often assumes
the training and test domain have identical conditional outcome distributions (PYtr|X = PYte|X ) but different marginal
tr
te
covariate distributions (PX
6= PX
). The popular importance weighting (IW) approach [26, 13, 29, 5] to solving
this covariate shift problem fails in high dimensions, in
tr
te
which the density ratios PX
(x)/PX
(x) used to reweight
the training loss are unreliably estimated and concentrated
on a small number of training samples [9], leading to high
variance estimates of a predictor’s expected test loss, the
quantity to be minimized over predictors. In this work, we
propose a extreme linear dimension reduction preprocessing
step, i.e. to very low dimensions (≤ 4 and often 1 in our
experiments), to improve IW methods’ performance on high
dimensional data.

Cynthia Rudin
Duke University
cynthia@cs.duke.edu

We propose this dimension reduction be extreme to maximally reduce the IW estimator variance due to density ratio
estimation error. We propose this dimension reduction be
linear so that the downstream IW method, which we will
assume to learn linear models, would still learn models linear in the original covariates, thus retaining interpretability.
We propose that some dimension reduction should be done,
so that we might uncover a subspace within which the train
and test covariate distributions were not too different. This
would afford the downstream IW estimator a larger effective
sample size and thus lower variance.
It is not clear a priori that extreme dimension reduction
would be computationally and experimentally favorable.
Firstly, given the extreme low dimension of the subspaces
we consider, in selecting one, we wish to estimate subspaces’
true predictive utility: the minimum achievable test domain
loss of a predictor acting on it. However, estimating this
utility for just a single candidate subspace requires solving
a convex optimization problem. We must confront this
potentically impractical computational problem; a cheaper
but less accurate proxy for that utility is unlikely to suffice
in extremely low dimensions. Secondly, as we point out,
dimension reduction introduces bias into the downstream
IW estimator, so that extreme dimension reduction might
add too much bias to actually help.
To test the utility of extreme dimension reduction, we first
develop a procedure that maximizes, over projection matrices, estimates of the resulting subspace’s aforementioned
true predictive utility, regularized by a penalty against subspaces where the downstream IW estimator would have high
variance. This procedure, using techniques from bilevel optimization, is verified to be computationally feasible. We
then study on real data whether it can identify subspaces
where the downstream IW method can reliably have low test
loss.
We find that on some, but not all real data, extreme dimension reduction can facilitate the learning of accurate
predictors for the test domain. We hypothesize that for such
success cases, our procedure returns sufficiently predictive
subspaces within which the reduction in downstream IW

estimator variance offsets the potential estimator bias. To
better understand this phenomena, we construct a simulated
data example where extreme dimension reduction helps, and
another one where the introduced bias is crippling. Finally,
reformulating the problem of learning a subgroup-specific
model as a covariate shift problem, we demonstrate the benefit of extreme dimension reduction in learning a depression
classifier for the college-aged population subgroup.

2
2.1

Background
Covariate Shift Problem

In the covariate shift problem, one is given N tr latr
tr tr
tr
belled training samples {xtr
i , yi }, with xi , yi ∼ PX,Y ,
te
the training distribution, and N unlabelled test samte
te
ples {xte
i }, with xi ∼ PX , the test distribution, with
tr
te
D
xi , xi ∈ R . Importantly, one makes the covariate
tr
te
shift assumption, that PYtr|X = PYte|X , but PX
6= PX
.
Given a model class F, the covariate shift problem seeks
te [l(f (X), Y )], where l is a loss function
argminf ∈F EPX,Y
that we assume to be convex in f (X).
2.2

Importance Weighted Loss Minimization

To minimize expected test loss over predictors, noting that
te
PX
(X)
te [l(f (X), Y )]
tr [ tr
EPX,Y
= EPX,Y
l(f (X), Y )], past
PX (X)
work constructs an unbiased estimator of test loss by forming the empirical expectation version of the latter expectation, and then minimizes it over predictors, adding some
regularization:
X
tr
argminf ∈F
ŵi l(f (xtr
(1)
i , yi ) + Ω(f )
i

where ŵi estimates wi :=

te
PX
(xtr
i )
tr (xtr ) .
PX
i

We note that the “effecP 2
i ŵi [10].

tive sample size” of the estimator is N tr /
2.3

Density Ratio Estimation

To carry out importance weighting requires estimating the
P te (xtr )
density ratios PXtr (xitr ) . A variety of methods [13, 5] exi
X
ist for doing so, but one, least squares importance fitting
(LSIF)[17] and its computationally more efficient variant,
unconstrained LSIF (uLSIF) stands out because it involves
only a (possibly constrained) quadratic program, and admits
a cross validation scheme.
LSIF assumes the estimated ratio function can be written as αT φ(x), where φ(x) = {φm (x)} is a set of M
basis functions they choose to be gaussian kernels centered at test points and α ∈ RK is a parameter they aspire to identify by minimizing the expected squared erP te (X)
ror to the true ratios: EPXtr [(αT φ(X) − PXtr (X) )2 ] =
X

EPXtr [(αT φ(X))2 ] − 2EPXte [αT φ(X)] + C, where C is

a constant. Substituting empirical expectations, dividing by
2, and regularizing, they obtain the fit value of α as
α∗ = argminα 12 αT Hα − hT α + γ1T α subject to α ≥ 0,
P
P
1
tr T
te
with H = N1tr k φ(xtr
i )φ(xi ) , h = N te
k φ(xi ),
γ is the regularization constant, and 1 is the vector of all 1’s.
The estimated ratio function is then ŵ(x) = α∗T φ(x), with
the element-wise positive constraint α ≥ 0 ensuring ŵ(x)
is positive. uLSIF simply removes the positivity constraint
on α in obtaining α∗ . With αT φ(x) no longer guranteed
to be positive, under uLSIF, the learned ratio function is
ŵ(x) = max(α∗T φ(x), 0).

3
3.1

Extreme Dimension Reduction for
Importance Weighting
Motivation

Our dimension reduction approach is motivated by the following: suppose we will apply a linear projection prior to an
IW method for handling covariate shift. Furthermore, suppose (as we will in this work) that the model class of the IW
method is the set of linear models. That is, we first choose
subspace dimension K, a projection matrix A ∈ RK×D ,
and then minimize over b ∈ RK (with regularization)
X

tr
ŵiA l bT AT xtr
, where (2)
L̂(b; A, ŵA ) := N1tr
i , yi
i

ŵiA estimates wiA :=

te
T tr
PA
T X (A xi )
P trT (AT xtr
i )
A X

and ŵA := {ŵiA }.

This projection by A can potentially reduce the variance of
the downstream IW loss estimator if we choose A to satisfy
two criteria. Firstly, if the dimension of the subspace A
projects to is greatly reduced, the density ratio estimates
themselves will be of lower variance, being those in a lower
dimensional space. Secondly, we can choose A so the effective sample
P size of the IW estimator following projection by
A, N tr / i (ŵiA )2 , would be above some threshold.
Of course, projection by A might hurt performance if other
criteria are not satisfied. We do not want projection by A to
introduce too much bias in the downstream IW estimator,
so that we may not want A to be too low-dimensional. And
of course we want A to have high predictive utility;
U (A) := argminb∈RK EP te [l(bT AT X)]

(3)

should be low.
Constructively, we can obtain an A satisfying a combination
of those criteria by minimizing over projection
P matrices an
estimator of U (A) subject to a regularizer i (ŵiA )2 that
discourages subspaces where the downstream IW estimator
would have high variance. We can estimate U (A) with
Û (A) := L̂(b∗ ; A, ŵA ), where
∗

(4)
A

2

b = argminb∈RK L̂(b; A; ŵ ) + ckbk ,

(5)

i.e. b∗ is a linear predictor with low out-of-sample test
loss, obtained via IW loss minimization regularized by a
tradeoff constant c. To complete the problem, we still need
to construct the projected density ratio estimates ŵiA , which
notably depend on A. We can do so by constraining that the
ŵiA are equal to the projected density ratio estimates that
would be returned by a ratio estimation method, i.e. LSIF or
T te
uLSIF, run on the projected covariates {AT xtr
i }, {A xi }.
3.2

Formulation

Here, we describe an optimization problem for finding a
projection matrix A ∈ {RK×D , AT A = I} with which
to project the covariates prior to the application of an IW
method for tackling the covariate shift problem, as described
in Section 3.1. K < D is the dimension of the desired
subspace. In particular, we assume given A, a linear model
will be fit using importance weighting, i.e. minimizing over
b ∈ RD the estimator of Equation 2 with regularization,
and that uLSIF will be used for density ratio estimation.
We propose finding A by solving:
X
X
min
L̂(b∗ ; A, ŵA ) + λ
(ŵiA )2
A

i

(6)

3.3

Analysis of Dimension Reduction for Importance
Weighting

We first show that in expectation, following dimenK×D
sion reduction
,
P A 2 by any projection matrix A ∈ R
tr
N / i (wi ) ), the effective sample size of the downstream IW loss estimator of Equation 2, where the density ratio estimates are replaced by the true density ratios, is larger
than the effective sample size of the IW loss estimator
had no
P
2
dimension reduction been performed, N tr / i (w(xtr
i ) ):
K×D
T
Lemma 3.1.
s.t.
P ForA 2A ∈ R
P A trA2 =
tr
tr
EPXtr [N / i (wi ) )] ≥ EPXtr [N / i (w(xi ) )].

I,

Pf. Let C be a matrix whose columns span the orthogonal
complement to the subspace spanned by the columns of A,
and let (U, V ) = (AT X, C T X). Then
X
P te (X) 2
2
tr [( tr
(w(xtr
(13)
EPXtr [ N1tr
i ) )] = EPX
P (X) ) ]
i

te

te

P (U )P (V |U ) 2
tr [( tr
= EPU,V
) ]
tr
i
h Pte (U )P (V |U ) te
P (V |U ) 2
P (U ) 2
= EPUtr ( P tr (U ) ) EPVtr|U [( P
tr (V |U ) ) ]
X
P te (U )
1
(wiA )2 ]
≥ EPUtr [( P
tr (U ) )] = EP tr [ N tr
X

(14)

(15)
(16)

i

i

subject to

b∗ = argminb∈RK L̂(b; A, ŵA ) + ckbk2
ŵiA = max(α∗T φ(utr
i ), 0)
∗
1 T
α = argminα 2 α Hα0 − hT α + γ1T α
X
T tr T
H = N1tr
φ(AT xtr
i )φ(A xi ) ,
h=

1
N te

k
X

AT A = I,

φ(AT xte
i )

(7)
(8)
(9)
(10)
(11)

k

(12)

where L̂(b∗ ; A, ŵA ) is as defined in Equation 2 and λ, γ, c
are hyperparameters described shortly.
The first term of the objective is the predictive utility estimate for A from Section 3.1, Û (A). It has such an
interpretation firstly due to the constraint of Equation 7,
and secondly because by Equations 8-11, the ŵiA are the
te
te
PA
T X (xi )
P trT (xte
i )
A X
T te
{AT xtr
},
{A
x
},
the
i
i

estimates of

returned by running uLSIF on

training and test covariates, projected by A. Again, note that the ŵiA thus depend on A
through an optimization problem: that of running uLSIF.
We use uLSIF instead of LSIF for ratio estimation due to its
computational benefits and relatively strong performance.
Due to the interpretation of the ŵiA , the second term of
the objective is a regularizer that encourages
the effective
P
sample size of that estimator, N tr / i (ŵiA )2 , to be high.
Regarding hyperparameter selection, we update c and γ
throughout individual runs of the gradient descent procedure we use with “in-line” cross-validation, and we choose
λ, controlling the lower bound on effective sample size, using weighted cross-validation [28] (details in supplement).

where we have used the fact that for any fixed U , 1 ≤
P te (V |U )

EPVtr|U [( PVtr|U (V |U ) )2 ] := P E(PVte|U kPVte|U ) + 1 to justify
V |U

Equation 16, which follows because P E(·k·), denoting the
Pearson Divergence between distributions, is always at least
0. This verifies our intuition that projecting onto a subspace
can, under idealized circumstances, increase effective sample size. Though, as the effective sample size following
projection is only larger in expectation, the realized effective sample size might still be small for some projections,
thus the need to explicitly regularize against such cases.

While dimension reduction reduces downstream estimator
variance (in expectation), it introduces bias in the estimator.
Under the covariate shift assumption, the estimator of Equation 2, with true density ratios replacing estimates thereof,
does not unbiasedly estimate test loss.

A
T T
tr [L̂(b; A, w )] 6= EP te [L b A X, Y ]
Lemma 3.2. EPX,Y
X,Y
for b ∈ RD , A as in Lemma 3.1.
tr
tr
Pf. Let U, V be as in Lemma 3.1, (xtr
i , yi ) ∼ PX,Y . Then,
X

A
A
T T tr tr
1
tr [L̂(b; A, w )] := EP tr [
EPX,Y
w
L
b
A
x
,
y
]
tr
i
i
i
N
U,Y
tr
= EPU,Y

P te (U )
[ PUtr (U ) L(f (U ), Y )]
U
P te (U )
[ PUtr (U ) EPY |U,V [L(f (U ), Y

i

tr
= EPU,V
)]]
U
Z
te
P (U )
= PUtr (U )PVtr|U (U ) PUtr (U ) EPY |U,V [L(f (U ), Y )]]dU dV
U
Z
te
tr
= PU (U )PV |U (U )EPY |U,V [L(f (U ), Y )]]dU dV

T T
te [L b A X, Y ].
= EPUte PVtr|U PY |U,V [L(f (U ), Y )] 6= EPX,Y

To understand the source of this bias, we can upper bound
the bias with an interpretable quantity (proof in supplement):
Lemma 3.3. Given A, b, let U, V be as defined in Lemma
tr
tr
3.1, and (xtr
i , yi ) ∼ PX,Y . Then
A

T

T

tr [L̂(b; A, w )] − EP te [l(b A X, Y )]
EPX,Y
X,Y

1

1

≤ EPUtr [P E(PVte|U kPVtr|U )] 2 EPUte [V arPYtr|U (l(bT U, Y ))] 2 .

The LHS is the absolute bias of our estimator. The first term
of the bound measures, roughly speaking, how much the
training and test distributions differ in the subspace orthogonal to the subspace A parameterizes. The second measures
how much variance there is in l(bT U, Y ), with U fixed, averaged over the test distribution of U , i.e. it represents how
well the loss can be predicted, given only U . The bias upper
bound is high if both terms are simultaneously high.
This analysis studies the bias and variance of the downstream IW estimator L̂(b; A, ŵA ), under the assumption that
A is fixed. But of course A is not fixed, and furthermore,
its selection uses that same IW estimator when obtaining a
subspace’s estimated predictive utility Û (A) of Equation 4.
Therefore, the bias and variance in L̂(b; A, wA ) also affects
the selection of A. In particular, the total variance of our
procedure contains the variability in selecting A. This is
why in experiments, the loss variance of our two-step procedure can be sometimes higher than IW, even if conditional
on A being chosen, we would expect the downstream IW
method to have lower variance. Also, of keen interest would
be a uniform bound on bias, over all A and b.
3.4

Solving the Optimization Problem

The objective of the optimization problem, which we will
refer to as G, depends only on projection matrix A, through
several intermediate variables. Thus, we solve the optimization problem via gradient descent over projection matrices
using the Pymanopt [30] package. The challenge is in cal∗
culating the gradient dG
dA , as intermediate variables α and
∗
b depend on other variables not analytically, but as the
solution to convex optimization problems parameterized by
the dependent variables. We call them argmin variables.
Fortunately we can apply existing work [11, 22] to efficiently calculate dG
dA . We used Autograd[23] to calculate
gradients not involving argmin variables. As the problem
is not convex in A, we utilize multiple random restarts.
Differentiation with argmin variables: In reverse mode
differentiation, the generic task is to recursively, given an
dv dG
objective of the form G(v(u)) and dG
dv , to compute du dv . In
our problem, we use reverse mode differentiation to compute
db∗ dG
dG
N tr vector of
dw db∗ given db∗ , where w is the length
dα∗ dG
dG
A
weights {ŵi } of Equation 8, and dA dα∗ given dα
∗ . For
concreteness, we will illustrate the latter calculation; the

same technique suffices for both, as both argmin variables
are the solutions to unconstrained problems. Please see [2]
∗
dG
∗
for how to calculate dα
dA dα∗ if LSIF is used, where α is
given by a constrained optimization problem.
∗

Calculating db
to explicitly form
dw : We first describe how
db∗
db∗ dG
dw , as the naive way to calculate dw db∗ is simply to form
db∗
∗
dw and then matrix multiply. Since b minimizes the convex
function f (b) := L̂(b∗ ; A, ŵA ) + ckbk2 , b∗ satisfies the
df ∗
stationarity condition db
(b (w), w) = 0, where 0 is a length
D vector of zeros, and the notation suggests f depends on
both b∗ and w, and b∗ further depends on w. Differentiating
∗
d df
d df
db∗
with respect to w gives db
dw ( db db ) + dw db = 0, where dw
d df
is the desired N × K Jacobian matrix, db
db is the K × K
df
d
Hessian matrix of f , and dw db is N × K. Rearranging,
∗
we obtain a multiple linear system we can solve for db
dw :
d df db∗ T
db ( db ) dw

d df T
= −( dw
db ) . However solving this multiple
∗
linear system is in general not feasible, as calculating db
dw
involves solving N separate linear systems, with w being
N tr -dimensional.
∗

dG
Efficient Calculation of db
dw db∗ : Fortunately, two computational tricks can be applied. Firstly, we can solve a single
linear system instead of N of them. We can express the
∗
dG
d df d df −1 dG
dG
desired gradient as db
dw db∗ = − dw db ( db db )
db∗ . As db∗
dG
d df
is assumed available, we can first solve ( db db )v = db∗ for
d df
db∗ dG
v, and then left multiply by − dw
db to get dw db∗ .

Moreover, we can avoid explicitly constructing the Hessian matrix in this linear system. Matrix free linear solvers
such as conjugate gradient, when solving a system Cx = d,
do not explicitly require C, but only that the matrix vector product Cu be able to be calculated for any vector u.
Furthermore, Hessian-vector products can be calculated efficiently. Returning to the example, note that for any u,
d df
d
T df
( db
db )u = db (u db ). This means the requisite Hessian
vector can be computed by first analytically forming the
df
scalar valued function uT db
and then calculating its gradient with respect to b analytically or via finite differences.

4

Simulation Study

We now apply extreme dimension reduction to two synthetic
examples in order to give insight into when the bias in the
IW estimator (c.f. Lemma 3.2) is detrimental. In the first,
this bias is nonexistent. In the second, this bias is crippling.
4.1

Example with no estimator bias

We generated 12-dimensional covariates and real valued
labels. The labels only depend on the first two covariates,
X1 and X2 . In particular, we let Y |X ∼ N (0.2kX1 k +
tr
kX2 k, 0.01). Figure 1a shows samples from PY,X
and
1
te
tr
te
PY,X1 , and Figure 1b shows samples from PY,X2 and PY,X
.
2
Training samples are in red, test samples are in test. Thus

Method
JP(1)
UW
IW
SIR(1)
RP(1)

50
.16(.08)
.26(.04)
.26(.04)
.26(.03)
.26(.04)

Data size
100
.13(.04)
.26(.03)
.26(.02)
.26(.03)
.26(.03)

150
.13(.03)
.26(.02)
.25(.02)
.26(.03)
.25(.02)

200
.14(.03)
.25(.02)
.25(.02)
.26(.02)
.25(.02)

Table 1: Example 1 - test loss over 50 replicates.

(a) X1 is moderately predic- (b) X2 is very predictive in
tive in both the test and train- the test distribution, but not
ing distribution.
very predictive in the training
distribution.

Figure 1: Example 1 - the relative importance between the
only 2 predictive features is reversed between training &
test distributions. Training points in red, test in blue.
for X1 and X2 , the test distribution contains only positive
values. The distributions of each covariate are independent
of each other, generated as follows:
X1tr , X2tr ∼ Uniform(−1, 1)
X1te , X2te ∼ Uniform(0, 1)
tr
X3tr , . . . , X12
∼ 0.9 Uniform(−1, 0) + 0.1 Uniform(0, 1)
te
te
X3 , . . . , X12 ∼ 0.1 Uniform(−1, 0) + 0.9 Uniform(0, 1)
Note that as our model class contains only linear models,
when trained and tested on P tr , a model that only uses X1
will have comparable performance compared to a model that
only uses X2 . This is because even though X2 appears more
informative, the model class does not contain the “v-shaped”
predictor that the mean of PY |X2 follows. Contrarily, when
trained and tested on P te , a model that uses only X2 will
have far superior performance to one that uses only X1 ,
because due to the covariate shift, the test covariates only
have support under one of the two “arms” of the “v”.
We compare our method (dimension reduction, then applying uLSIF), denoted JP(K), K the subspace dimension, to:
• Unweighted (UW): performing no covariate shift correction, minimizing (regularized) unweighted loss over
the training data
• Naive Importance Weighting (IW): Not applying any
dimension reduction, and then applying uLSIF.
• Random Projection to K dimensions (RP(K)): Apply-

ing a random projection to generate K-dimensional
covariates, then applying uLSIF to the projected data.
• Sliced Inverse Regression (SIR(K)): Applying sliced
inverse regression [19] to project to K-dimensions, and
then applying uLSIF.
Table 1 shows the mean and standard deviation of out-ofsample test loss (measured by absolute prediction error)
over 50 replicates for the methods, as N , the number of
generated training and test samples changes. To elaborate,
tr
te
we generate N samples each from PX,Y
and PX
, set aside
1
of
the
test
data,
fit
a
model
using
the
N
labelled
training
3
samples and remaining 2N
unlabelled
test
samples,
and
3
evaluate the predictions on the set aside test data.
We note the performance of IW is often worse than that
of UW, due to the difficulty in estimating the weights and
small sample size. Furthermore, the sliced inverse regression methods, which do not account for covariate shift in
estimating the projection, do not find any subspace particularly informative and thus do about the same as the random projection methods. On the other hand, our method
is able to find that X2 is most useful in the test domain, so
that the downstream IW method is applied to the correct
1-dimensional subspace. We note that the variance in the
loss for our method is slightly increased, due to the variance
in the selection of the subspace within which to apply IW.
Our dimension reduction method projects the covariates
to retain only X2 , because it finds Û (X2 )  Û (X1 ) (by
abuse of notation we will use X2 also to refer to the 1dimensional projection matrix that retains only X2 ). Why
is it able to accurately estimate Û (X1 ) and Û (X2 ) despite
the fact that their estimation uses estimators of a predictor b’s test loss, L̂(b; X1 , ŵX1 ) and L̂(b; X2 , ŵX2 ), that are
potentially biased, by Lemma 3.2? The reason is that we
have designed the data distributions so that L̂(b; X1 , wX1 )
and L̂(b; X2 , wX2 ), the same estimators but with true density ratios replacing their estimates, are actually unbiased.
We show this in the case of X1 . Since L̂(b; X1 , wX1 ) is a
function of just X1 , X2 , Y , in calculating its expectation
with respect to any distribution, we can marginalize over
variables other than X1 , X2 , Y . Specifically, in calculating EP tr [L̂(b; X1 , wX1 )], we can invoke Lemma 3.2 with
U = X1 , V = X2 , to obtain
EP tr [L̂(b; X1 , wX1 )] = EP 1\2 [L̂(b; X1 , wX1 )],
1\2
te
where PX1 ,X2 ,Y := PX
P tr
P
.
1 X2 |X1 Y |X1 ,X2

(17)
(18)

tr
te
However, actually PX
P
= PY,X
(due
2 |X1 Y |X1 ,X2
2 |X1
1\2

to the symmetry in kX2 k). This implies PX1 ,X2 ,Y =
te
by Equation 18. Specifically, expectations with
PX
1 ,X2 ,Y
respect to these 2 distributions are equal, so that by Equation 17, EP tr [L̂(b; X1 , wX1 )] = EP te [L̂(b; X1 , wX1 )], i.e.
L̂(b; X1 , wX1 ) unbiasedly estimates the test loss of b.

Dataset
Example 1
Example 2

50
.16(.08)
.26(.10)

Data size
100
150
.13(.04) .13(.03)
.20(.11) .23(.11)

200
.14(.03)
.22(.13)

Table 2: Comparison of our method’s performance on the 2
simulated examples.
tr
te
(a) PX
and PX
1 ,Y
1 ,Y

tr
te
(b) PX
and PX
2 ,Y
2 ,Y

Figure 2: Example 2 - P te is unchanged from before, so that
X2 is still most predictive in the test distribution. Training
points in red, test in blue.

1\2

(a) PX1 ,Y

2\1

(b) PX2 ,Y

Figure 3: Example 2 - Our method erroneously finds X1
and X2 to be equally predictive in the test distribution.
4.2

5

Example with crippling estimator bias

Now, we slightly modify the data distribution from the previous example, and obtain an example where the previously
unbiased estimators suffer from crippling bias. In particular,
we make a single change, letting
(
X1tr wpr 0.5
tr
tr
X2 |X1 =
(19)
−X1tr otherwise.
For this second example, Figure 2a shows samples from
tr
te
tr
PY,X
and PY,X
, and Figure 2b shows samples from PY,X
1
1
2
te
te
and PY,X2 . P is unchanged, and X1 and X2 remain not
predictive in P tr . However, the important difference in this
1\2
2\1
example is that in PX1 ,X2 ,Y and PX1 ,X2 ,Y , where by the
same reasoning of Equations 17 and 18,
EP tr [L̂(b; X2 , wX2 )] = EP 2\1 [L̂(b; X2 , wX2 )],
2\1
te
where PX1 ,X2 ,Y := PX
P tr
P
.
2 X1 |X2 Y |X1 ,X2

1-dimensional projection will contain equal contributions
from X1 and X2 . As P te is unchanged between the two examples, it remains the fact that the ideal projection to apply
before learning a model for P te would retain only X2 . Any
other projection would be suboptimal. This is reflected in
Table 2, where the performance of our method is degraded,
so that they approach that of the other baseline methods in
the first example. We can verify the reason for this degradation by examining the component of the 1-dimensional
projection vector corresponding to X1 and X2 , for both
examples. We do so when N tr = N te = 200, i.e. fairly
large. The average value of those 2 components for the
1st example are 0.01, 0.99, respectively. The corresponding
values for the 2nd example are 0.49, 0.53, respectively, reflecting that in the 2nd example, our method finds X1 and
X2 equally predictive.

(20)
(21)

Therefore, to understand the bias in L̂(b; X1 , wX1 ) and
1\2
L̂(b; X2 , wX2 ) for any b, we only need to examine PX1 ,Y

Experiments with Real Data

5.1
5.1.1

Covariate Shift Experiments
Data Preprocessing

We took several datasets from [20, 1], and for each, introduced covariate shift by creating a sampling scheme that
can repeatedly generate training and test data samples. For
each dataset, we first identify a single predictive vector in
covariate space by generating 100 random vectors, and measuring the predictive utility of a vector by projecting the
covariates onto it, running kernel density regression, and
examining the in-sample squared error. We first sample the
dataset to form the training data. Then to generate the test
data, we subsample from the remaining dataset according
to the projections along the vector. Let t1 , t0 denote the
max and min projected covariate value in the dataset, and
σ be the standard deviation of the projected values. We
select data with probability proportional to the density of
a N (α(t1 − to ), cσ 2 ) distribution based on their projected
value. Values of α close to 0 or 1, and small values of c will
lead to a small effective sample size; we choose them so
this is achieved. Given that vector, multiple training and test
datasets can then be sampled.

2\1

and PX2 ,Y (note the estimators are functions of X1 , Y and
1\2

X2 , Y , respectively). Figure 3a shows samples from PX1 ,Y ,
2\1

and Figure 3b shows samples from PX2 ,Y . What we see
is that those two distributions are identical. This means
that our method will find Û (X1 ) ≈ Û (X2 ), so that the

5.1.2

Experimental Results

For each dataset, we generated 50 training and test data pairs.
For each pair, we set aside 13 of the test data to evaluate outof-sample performance on the test distribution, and train

Figure 4: For several datasets, the mean loss (“loss”), downstream effective sample size (“N_eff”), and loss standard
deviation (“loss std”) of our method is shown in black as the subspace dimension varies. The respective values of the
unweighted and importance weighting baseline are shown in blue and red, respectively. Dataset dimension in parentheses.
the model with the remaining data. Experimental results
are in Figure 4, where for each dataset, in 3 subplots we
plot in black the respective quantities for our method as the
dimension of the projection changes: “loss” - the average
loss over the 50 replicates, “N_eff” - the average effective
sample size enjoyed by the downstream IW estimator, “loss
std” - the standard deviation in loss over the replicates. For
comparison, we also indicate with horizontal lines the values
of those 3 quantities for the unweighted (UW) baseline in
blue, and those for the naive importance weighting (IW)
baseline in red. We normalize the results so that UW has a
average loss and average effective sample size of 1. Thus for
the loss and N_eff subplots, there is always a horizontal blue
line at a height of 1. In parentheses are dataset dimension
and whether it is classification (C) or regression (R). We
use absolute prediction error loss for regression problems
as before, and 0-1 loss for classification problems.
We see that extreme dimension reduction does not always
help; whether it does depends on the dataset and the dimension of the subspace. However, for many datasets, there is
some subspace dimension for which our dimension reduction procedure does help. For a given dataset, the loss can
both increase (i.e. communities dataset) or decrease (i.e.
facebook) with a reduction in subspace dimension, depending whether the reduction in the variance of the estimator
L̂(b; A, wA ), used both in selecting A and the downstream

IW, is enough to offset a potential increase in estimator bias.
What dimension reduction offers to the covariate shift problem is a way to navigate that tradeoff. We also note that
as alluded to in Section 3.3 and the simulation study, the
variance of our two-step procedure can be higher than IW.
We also caution that the 3 quantities we report for our
method are not directly comparable between different subspace dimensions, because the value of hyperparameter λ
chosen by cross-validation also differs between dimensions.
This means we cannot make claims like “the loss standard
deviation is lowered when reducing the subspace dimension
solely due to decreased density ratio estimation variance”.
Despite this caveat, we can still identify a phenomena that
is paradoxical at first: when going from 2 dimensions to 1,
the utilized effective sample size, N_eff, tends to increase.
Coupled with the accompanying decrease in density ratio
estimation variance, we would expect the loss standard deviation to decrease under this change. However, we observe
exactly the opposite. This phenomena is likely due to the
increased number of bad local optima in the non-convex
optimization problem when the subspace dimension is extremely low, i.e. 1. We anecdotally observed the same
phenomena in the simulation study, where loss standard
deviation was sensitive to how many random restarts we
used for optimization.

(a) Our method has better rel- (b) Our method is able to use
ative performance when la- more data to learn the subbelled data is scarce.
group model than IW & SUB.

Figure 5: Our method performs relative better when labelled
data is scarce due to its ability to utilize more data for fitting.

5.2

Case Study on Learning Subgroup Models

Recent work has shown that some real world models[6]
have lower prediction performance on minority subgroups
than for the majority group. In theory, this problem could
be alleviated by learning a subgroup model by only using
data from the subgroup. However, as data are scarce, the
performance of the subgroup model could be improved by
utilizing the appropriate data from the majority group. In
correcting for covariate shift, a model is learned to minimize average prediction error with respect to one population
given labelled samples drawn from another population. This
motivates us to recast the subgroup model learning problem
as a covariate shift correction problem. We will describe this
formulation, then use our method to learn a classification
model for depression for the age 18-24 subgroup.
Formally, given (xi , yi , zi ) ∼ PX,Y,Z , where feature vectors xi ∈ RD and zi is an indicator denoting membership in the subgroup of interest, given loss function l and
model class F, the subgroup model learning problem seeks
argminf ∈F EPX,Y |Z=1 [l(f (X), Y )]. We reformulate it as
te
a covariate shift problem where PX,Y
= PX,Y |Z=1 and
tr
PX,Y = PX,Y . The covariate shift assumption is then that
PY |X,Z=1 = PY |X (which might not actually hold).
We apply our covariate shift correction method to the National Health and Nutrition Examination Survey (NHANES)
[16] dataset, downloaded via the eponymous R package
[24]. This dataset contains demographic, physical, health,
and lifestyle variables from 10000 individuals in the United
States, obtained via surveys between 1999 and 2004. Importantly, these individuals are sampled to be representative of
the United States general population as a whole, so that subgroups underrepresented in the population will also be so in
the dataset. Given the recent rise in mental health issues of
college students [3], we use our method to learn a classifier
to predict the presence of depression for the subgroup of
persons of college age, i.e. those between 18 and 24 years
of age. Such persons, according to the NHANES dataset,
are on average different from the general population. For
example, they tend to be less likely to smoke, and have sleep

troubles, more likely to be single and physically active, and
on average have lower weight and poverty levels. These differences open up the possibility that the best (linear) model
for this subgroup may be different than that of the general
population. Furthermore, only about 10% of the general
population (and thus NHANES dataset) falls into this subgroup, so that approaches beyond simply fitting a model to
data from the subgroup are needed to handle data scarcity.
Regarding features, we only use those that can be easily
obtained, such as those aforementioned as well as income,
pulse, years of education, excluding those require medical tests such as urine sugar level, and those too closely
correlated with depression, i.e. how often the individual
considers their own general mental health to be bad. The
utility of the classifier would be to identify those at risk
of depression, using easily obtainable features that are not
proxies for depression. Following this guideline, we end up
using 27 features, and after removing those for which not all
features were available, retain data from 4517 individuals,
of which 21% were labelled as being depressed (proportion
of days feeling depressed not equal to “None”).
We compare the performance of our method (projecting
to 1-dimension, denoted JP(1)) to the naive importance
weighting and unweighted methods (denoted IW and UW
in figures) in learning a classifier for depression for the age
18 to 24 subgroup using the covariate shift reformulation,
as well as the baseline of fitting the classifier using only
data from the subgroup (denoted SUB). We perform this
comparison as the number of labelled data varies - for various N , we uniformly subsample a subset of size N , train
a model using a given method, and obtain predictions for
the remaining unsampled data. We do this 100 times for
each N and method, and report the average 1-AUC for them
in Figure 5a as well as the average effective sample size
(N_eff) in fitting the models in Figure 5b.
We see that our method, for all dataset sizes, it is able to
consider only (1-dimensional) subspaces in which the effective sample size of the loss estimator is close to that of
the entire labelled data without hurting performance (recall
the tradeoff controlling the lower bound on effective sample
size is chosen via cross validation). When the labelled data
is scarce, the resulting low variance of the loss estimator
offsets the bias that results from only estimating density
ratios within a subspace, when compared to IW, whose effective sample size is always smaller. When data is plentiful,
i.e. at N = 1500, the relative advantage of our method in
utilizing more data is diminished compared to the baselines,
who utilize less biased estimators. A similar relationship
holds between IW and the more biased UW, though the
dataset size at which their relative performances change order occurs at N = 1000. Finally, SUB, whose estimator is
unbiased but enjoys far lower effective sample sizes, has far
worse performance than all of the covariate shift approaches,
which utilize labelled data from outside the subgroup.

6

Related Work

Our work studies the benefit of dimension reduction prior
to applying importance weighting approaches [13, 29, 5]
for solving the covariate shift problem [26]. In particular,
we use dimension reduction as a way to potentially reduce
the otherwise high variance in importance weighting approaches pointed out by [9]. Past work addressing this
issue [26, 9, 32, 14] has used various forms of regularization to discourage exceedingly large weights, at the cost
of increased bias. Another line of work eschews weight
estimation and builds predictors that are robust to potential
shifts in the covariate [31] or conditional outcome [21, 8]
distribution. A separate work [25] reduces variance by using the predictor minimizing training loss as a prior when
minimizing the reweighted loss. One work [27] does apply
dimension reduction before estimating density ratios, by
searching for a subspace in which the training and test densities are maximally different, and estimating the ratios within
it. However, their chosen subspace, by construction, would
result in a small underlying effective sample size relative to
other subspaces.
Our approach is similar to supervised representation learning methods that find representations that are both predictive
and similar in some sense between the two domains; penalizing subspaces with small effective sample sizes is equivalent
to penalizing subspaces in which the Pearson divergence is
high. However, none of those approaches perform importance weighting when evaluating the predictive utility of a
representation, which is advantageous in many situations
when the downstream model (linear model in our case) is
misspecified so that importance weighting is actually necessary to begin with. For example, [7] finds a linear subspace
minimizing an additive combination of maximum mean discrepancy, and the lowest possible unweighted squared loss
of a linear model acting on it. Recent work has used neural
networks to find representations with which low unweighted
training domain loss can be achieved, and which are similar
between domains, as accomplished through using closed
form discrepancy measures [15] or adversarial training [12].
However, the learned model is not linear and interpretable.

7

Conclusion

To address the high variance of importance weighting (IW)
approaches to handling covariate shift, we have explored the
benefit of extreme dimension reduction, i.e. to very low dimensions, prior to applying such approaches. We presented
and solved a bilevel optimization problem that searches for
a subspace that is predictive, and for which the downstream
IW procedure would enjoy a large effective sample size.
We illustrated through lemmas and simulated and real data
that dimension reduction helps sometimes but not always,
due to bias introduced into the IW loss estimator, and the
fact that we incur additional variance in the selection of the

subspace. Study into how to alleviate these issues, as well as
how to avoid bad local optima in our optimization problem,
are potential future work. We also formulated the problem of learning a subgroup-specific model as a covariate
shift problem, and demonstrate the advantage of dimension reduction in learning a classifier for depression for the
college-aged subgroup. Finally, we believe our work has
implications for causal inference, which can be formulated
as a covariate shift problem [15], and where the propensity
score, essentially a density ratio, is used widely [4], yet
unreliably estimated in higher dimensions [18].

References
[1] Libsvm datasets. http://www.csie.ntu.edu.
tw/~cjlin/libsvmtools/datasets/.
[2] Brandon Amos and J Zico Kolter. Optnet: Differentiable optimization as a layer in neural networks. arXiv
preprint arXiv:1703.00443, 2017.
[3] American Psychological Association. The Crisis on
Campus, 2011.
[4] Peter C Austin. An introduction to propensity score
methods for reducing the effects of confounding in observational studies. Multivariate behavioral research,
46(3):399–424, 2011.
[5] Steffen Bickel, Michael Brückner, and Tobias Scheffer.
Discriminative learning under covariate shift. Journal
of Machine Learning Research, 10(Sep):2137–2155,
2009.
[6] Joy Buolamwini and Timnit Gebru. Gender shades:
Intersectional accuracy disparities in commercial gender classification. In Sorelle A. Friedler and Christo
Wilson, editors, Proceedings of the 1st Conference on
Fairness, Accountability and Transparency, volume 81
of Proceedings of Machine Learning Research, pages
77–91, New York, NY, USA, 23–24 Feb 2018. PMLR.
[7] Bo Chen, Wai Lam, Ivor Tsang, and Tak-Lam Wong.
Extracting discriminative concepts for domain adaptation in text mining. In Proceedings of the 15th ACM
SIGKDD international conference on Knowledge discovery and data mining, pages 179–188. ACM, 2009.
[8] Xiangli Chen, Mathew Monfort, Anqi Liu, and
Brian D Ziebart. Robust covariate shift regression.
In Artificial Intelligence and Statistics, pages 1270–
1279, 2016.
[9] Corinna Cortes, Yishay Mansour, and Mehryar Mohri.
Learning bounds for importance weighting. In Advances in neural information processing systems,
pages 442–450, 2010.

[10] Arnaud Doucet, Nando De Freitas, and NJ Gordon. Sequential monte carlo methods in practice. series statistics for engineering and information science, 2001.
[11] Chuan-sheng Foo, Chuong B Do, and Andrew Y Ng.
Efficient multiple hyperparameter learning for loglinear models. In Advances in neural information
processing systems, pages 377–384, 2008.
[12] Yaroslav Ganin, Evgeniya Ustinova, Hana Ajakan,
Pascal Germain, Hugo Larochelle, François Laviolette,
Mario Marchand, and Victor Lempitsky. Domainadversarial training of neural networks. Journal of
Machine Learning Research, 17(59):1–35, 2016.
[13] Arthur Gretton, Alex Smola, Jiayuan Huang, Marcel Schmittfull, Karsten Borgwardt, and Bernhard
Schölkopf. Covariate shift by kernel mean matching.
2009.
[14] Weihua Hu, Issei Sato, and Masashi Sugiyama. Robust
supervised learning under uncertainty in dataset shift.
arXiv preprint arXiv:1611.02041, 2016.
[15] Fredrik Johansson, Uri Shalit, and David Sontag.
Learning representations for counterfactual inference.
In International Conference on Machine Learning,
pages 3020–3029, 2016.
[16] Clifford L Johnson, Ryne Paulose-Ram, Cynthia L
Ogden, Margaret D Carroll, Deanna Kruszan-Moran,
Sylvia M Dohrmann, and Lester R Curtin. National health and nutrition examination survey. analytic
guidelines, 1999-2010. 2013.
[17] Takafumi Kanamori, Shohei Hido, and Masashi
Sugiyama. A least-squares approach to direct importance estimation. Journal of Machine Learning
Research, 10(Jul):1391–1445, 2009.
[18] Joseph DY Kang and Joseph L Schafer. Demystifying
double robustness: A comparison of alternative strategies for estimating a population mean from incomplete
data. Statistical science, pages 523–539, 2007.
[19] Ker-Chau Li. Sliced inverse regression for dimension
reduction. Journal of the American Statistical Association, 86(414):316–327, 1991.
[20] M. Lichman. UCI machine learning repository, 2013.
[21] Anqi Liu and Brian Ziebart. Robust classification
under sample selection bias. In Advances in neural
information processing systems, pages 37–45, 2014.
[22] Dougal Maclaurin, David Duvenaud, and Ryan Adams.
Gradient-based hyperparameter optimization through
reversible learning. In International Conference on
Machine Learning, pages 2113–2122, 2015.

[23] Dougal Maclaurin, David Duvenaud, and Ryan P
Adams. Autograd: Reverse-mode differentiation of native python. In ICML workshop on Automatic Machine
Learning, 2015.
[24] Randall Pruim. NHANES: Data from the US National
Health and Nutrition Examination Study, 2015. R
package version 2.1.0.
[25] Sashank Jakkam Reddi, Barnabas Poczos, and Alexander J Smola. Doubly robust covariate shift correction.
2015.
[26] Hidetoshi Shimodaira. Improving predictive inference
under covariate shift by weighting the log-likelihood
function. Journal of statistical planning and inference,
90(2):227–244, 2000.
[27] Masashi Sugiyama, Satoshi Hara, Paul Von Bünau, Taiji Suzuki, Takafumi Kanamori, and Motoaki
Kawanabe. Direct density ratio estimation with dimensionality reduction. In Proceedings of the 2010
SIAM International Conference on Data Mining, pages
595–606. SIAM, 2010.
[28] Masashi Sugiyama, Matthias Krauledat, and KlausRobert Muller. Covariate shift adaptation by importance weighted cross validation. Journal of Machine
Learning Research, 8(May):985–1005, 2007.
[29] Masashi Sugiyama, Shinichi Nakajima, Hisashi
Kashima, Paul V Buenau, and Motoaki Kawanabe.
Direct importance estimation with model selection
and its application to covariate shift adaptation. In
Advances in neural information processing systems,
pages 1433–1440, 2008.
[30] James Townsend, Niklas Koep, and Sebastian Weichwald. Pymanopt: A python toolbox for optimization
on manifolds using automatic differentiation. Journal
of Machine Learning Research, 17(137):1–5, 2016.
[31] Junfeng Wen, Chun-nam Yu, and Russell Greiner. Robust learning under uncertain test distributions: Relating covariate shift to model misspecification. In
Proceedings of the 31st International Conference on
Machine Learning (ICML-14), pages 631–639, 2014.
[32] Makoto Yamada, Taiji Suzuki, Takafumi Kanamori,
Hirotaka Hachiya, and Masashi Sugiyama. Relative
density-ratio estimation for robust distribution comparison. In Advances in neural information processing
systems, pages 594–602, 2011.

Supplementary Material

Choosing Hyperparameters
The hyperparameters in the formulation fall into two groups: those whose optimal value does not change with A,
namely λ, and those whose optimal value does change with A: c, which prevents overfitting when evaluating the
predictive utility of A, and those required for density ratio estimation, namely γ as well as the bandwidth of the
Gaussian kernel basis functions, which we refer to as σ. We do not consider the latter group of hyperparameters
fixed, but instead update them throughout the gradient descent over A. We also note that by updating these
hyperparameters as the optimization runs, they are no longer hyperparameters in the conventional sense. The
reduced number of hyperparameters is thus computationally convenient.
The optimal values of γ and σ change with A because they are used to obtain estimates ŵiA of which changes
with A. Certainly, the hyperparameters of a density ratio estimation procedure should change as the true density
ratios change. To choose γ, σ for a given A, we use exactly the cross validation procedure [1] suggest, with 3-fold
cross-validation.
The optimal value of c is presumed to change with A for similar reasons. We would like to choose c such that given
A, b∗ , obtained in Equation 7 (i.e. IW loss minimization, regularized by c), equals argminb∈RK EP te [l(bT AT X)],
i.e. the best linear model for the test domain when covariates are first projected by A. Since c is the regularization
needed to solve a A-dependent learning problem, c should depend on A. To evaluate out-of-sample test domain
performance for a given A and linear predictor b, we use 3-fold weighted cross-validation as described in [2], where
a weighted average of out-of-sample labelled training data samples is used to approximate test domain loss. A
key difference is that given A, we estimate density ratios of the covariates, projected by A, as the density ratios in
the full space would not be estimated reliably.
While σ, γ are updated throughout a single optimization run, λ is held fixed, and thus chosen by “traditional”
out-of-sample validation. In particular, to evaluate a given λ, we use 5-fold weighted cross-validation, similar
to the previous section; using in-sample labelled training domain and unlabelled test domain data, we obtain
an A using the gradient descent procedure, with c, γ, σ updated through out it. The procedure returns an A,
and we use that A to calculate the density ratios of projected out-of-sample training data, which we then use to
approximate test domain loss.
In summary, in λ, the only “conventional” hyperparameter, is chosen to be the value whose resultant predictors
have the lowest average value of a loss estimate computed using the validation sets from 5-fold cross-validation
and estimates of projected density ratios. The remaining hyperparameters γ, σ, c are updated throughout the
optimization procedure for a fixed λ using grid search 3-fold cross-validation.

Proof of Lemma 3.1.2
Lemma. Given A, b, let C be a matrix whose columns span the orthogonal complement to the subspace spanned
tr
tr
by the columns of A, and let (U, V ) = (AT X, C T X) and (xtr
i , yi ) ∼ PX,Y . Then
1

1

A
T T
te
tr
T
2
tr [L̂(b; A, w )] − EP te [l(b A X, Y )] ≤ EP tr [P E(P
te [V arP tr (l(b U, Y ))] 2 ,
EPX,Y
V |U kPV |U )] EPU
X,Y
U
Y |U
1

where
L̂(b; A, wA ) :=

1
N tr

X

tr
wiA l bT AT xtr
i , yi

i

wiA :=

te
T tr
PA
T X (A xi )
P trT (AT xtr
i )
A



(1)
(2)

X

Pf. We first prove an upper bound of the bias of our loss estimator.
X P teT (AT xtr )
te
,V )
i
tr
T
A X
tr
EPXtr,Y [ N1tr
l(bT AT xtr
[ P tr (U
i , yi )] − EPU
(U ,V ) l(b U , Y )]
P tr (AT xtr )
,V ,Y P
AT X

i

P te

(3)

i

(AT X )

te

P (U ,V )
T
= EPXtr,Y [ PAtrT X (AT X ) l(bT AT X , Y )] − EPUtr,V ,Y [ P
tr (U ,V ) l(b U , Y )]
AT X
te

(4)

te

P (U )
,V )
T
T
= EPUtr,V,Y [ P
[ P tr (U
tr (U ) L(b (U ), Y )] − EP tr
(U ,V ) l(b U , Y )]
U ,V ,Y P

(5)

P te (U ,V )
T
[P
tr (U ,V ) l(b U , Y

P te (U )
T
[P
tr (U ) l(b U , Y

)] − EPUtr,V ,Y
)]
(6)
= EPUtr,Y
h
i
(7)
= EPUte EPVtr|U [EPY |U,V [l(bT U, Y )]] − EPVte|U [EPY |U,V [l(bT U, Y )]]
hZ
i

= EPUte
(PVtr|U (V ) − PVte|U (V ) EPY |U,V [l(bT U, Y )]dV
(8)
V
Z
Z
h
i



= EPUte
PVtr|U (V ) − PVte|U (V ) EPY |U,V [l(bT U, Y )] − µU dV +
PVtr|U (V ) − PVte|U (V ) µU dV
(9)
V
V
Z
hZ
 Z
i

P te (V ) 
= EPUte
1 − PVtr|U (V ) EPY |U,V [l(bT U, Y )] − µU PVtr|U (V )dV + µU −
PVtr|U (V )dV −
PVte|U (V )dV
V |U

V

h Z

PVte|U (V
PVtr|U (V

) 2 tr
) PV |U (V
)

1  Z
2

V

(EPY |U,V [l(bT U, Y )] − µU )2 PVtr|U (V )dV
V
V
h
1 i
1 
varPVtr|U (EPY |U,V [l(bT U, Y )]) 2
= EPUte PE(PVte|U ||PVtr|U ) 2
≤ EPUte

(1 −

)dV

V

1
2

+ µU 1 − 1

i


h
h
i 1 
i 1
2
2
≤ EPUte PE(PVte|U ||PVtr|U )
EPUte varPVtr|U (EPY |U,V [l(bT U, Y )])

(10)
(11)
(12)

where µU := EPVtr|U [EPY |U,V [l(bT U, Y )]] in Equation 10 and 12 I used Cauchy-Schwartz and PE(pX ||qX ) :=
X (X)
2
EqX [( pqX
(X) − 1) ], is the Pearson Divergence between distributions. We can now upper bound the second term
in Equation 12.

By the law of total variance, for a given U ,
varPVtr|U (EPY |V,U [l(bT U, Y )]) = varPYtr|U (l(bT U, Y )) − EPVtr|U [varPY |U,V (l(bT U, Y ))].
Taking expectations of this equality over PUte gives
EPUte [varPVtr|U (EPY |U,V [l(bT U, Y )])] = EPUte [varPYtr|U (l(bT U, Y ))] − EPUte [EPVtr|U [varPY |U,V (l(bT U, Y ))]]
≤ EPUte [varPYtr|U (l(bT U, Y ))]

(13)
(14)

due to the second term on the RHS of Equation 13 always being at least 0. Combining Equation 14 and 12 gives
the upper bound. The proof of the lower bound of the bias is entire analogous, except it uses the same steps to
upper bound the negative of the bias, and we omit it for brevity.

References
[1] Takafumi Kanamori, Shohei Hido, and Masashi Sugiyama. A least-squares approach to direct importance
estimation. Journal of Machine Learning Research, 10(Jul):1391–1445, 2009.
[2] Masashi Sugiyama, Matthias Krauledat, and Klaus-Robert Muller. Covariate shift adaptation by importance
weighted cross validation. Journal of Machine Learning Research, 8(May):985–1005, 2007.

