Classical risk minimization (RM) minimizes the expected loss over the training distribution $p_{\mathrm{train}}(x)$,
$$\theta^*_{RM} = \arg \min_\theta E[\ell(x, \theta)]_{p_{\text{train}}}.$$
As the distribution $p_{\text{train}}$ is usually unknown, empirical risk minimization (ERM) replaces the analytical expectation with the sample expectation, i.e.,
$$\theta^*_{ERM} = \arg \min_\theta \frac{1}{N}\sum_{i=1}^N \ell(x_i, \theta).$$
Now suppose the training data are drawn from a distribution $p_{\mathrm{train}}(x)$, but you would like the model to perform well on data drawn from another distribution $p_{\mathrm{target}}(x)$. This is what's called "covariate shift", and in this case, we would like to minimize the expected loss over the target distribution, giving rise to importance-weighted risk minimization (IWRM)
$$ \theta^*_{IWRM} = \arg \min_\theta E[\ell(x, \theta)]_{p_{\text{target}}} \\ = \arg \min_\theta E\left[\frac{p_{\text{target}}(x)}{p_{\text{train}}(x)}\ell(x, \theta)\right]_{p_{\text{train}}} $$
and its finite-sample counterpart, importance-weighted empirical risk minimization (IWERM):
$$\theta^*_{IWERM}= \arg \min_\theta \frac{1}{N}\sum_{i=1}^N \underbrace{\frac{p_{\text{target}}(x_i)}{p_{\text{train}}(x_i)}}_{=w_i} \ell(x_i, \theta).$$
In practice, this amounts to simply weighting individual samples by their importance $w_i$.
Shimodaira (2000) proves that (assuming $p_{\text{train}}$ and $p_{\text{target}}$ to be known) the IWERM estimator is asymptotically unbiased, i.e., for $N\to\infty$ we have
$$\theta^*_{IWERM} = \arg \min_\theta E[\ell(x, \theta)]_{p_{\text{target}}}.$$
My question is: why does this not hold for finite sample numbers as well? I would think that the estimator should also be unbiased for finite $N$, but that does not appear to be the case. Can someone explain why?
Here's an example demonstrating the biasedness of IWERM estimation for $N<\infty$. $N=500$ samples $x_i$ are drawn from a Gamma distribution, truncated to [0,10] ($p_{\text{train}}$), and $p_{\text{target}}$ is the uniform distribution over the same interval. A line is fit to a 4th-order polynomial $f(x)$ (with additive Gaussian noise of constant variance), and the cost function is the squared error. The figure shows one realization of the data; the regression lines are averaged over 1000 realizations. OLS = (unweighted) ordinary least squares, IWLS = importance-weighted least squares, ideal = solution that minimizes the expected squared error over $p_{\text{train}}$, the uniform distribution. For larger values of $N$, the IWLS solution converges to the "ideal" solution, in accordance with the asymptotic unbiasedness shown by Shimodaira.
(The intro to this question is adapted from my answer here.)
The full R code to reproduce the above figure can be found here.