A Unified Framework for Multi-distribution Density Ratio Estimation

In multi-class experiments, we have a pair of random variables $(X,Y)\in{\mathcal{X}}\times{\mathcal{Y}}$ with joint distribution $D(X,Y)$, where ${\mathcal{X}}$ is the sample space and ${\mathcal{Y}}=[k]\mathrel{\mathop{:}}=\{1,\ldots,k\}$ is the finite label space. Define the probability simplex as $\Delta_{k}\mathrel{\mathop{:}}=\{{\bm{p}}\in{\mathbb{R}}^{k}_{\geq 0}|\mathbf{1}^{\top}{\bm{p}}=1\}$. According to chain rule of probability, any joint distribution $D(X,Y)$ can be decomposed into class priors $\pi_{i}\mathrel{\mathop{:}}={\mathbb{P}}(Y=i)$ and class conditionals $P_{i}(x)\mathrel{\mathop{:}}={\mathbb{P}}(X=x|Y=i)$ for $i\in[k]$, or into sample marginal $M(x)\mathrel{\mathop{:}}={\mathbb{P}}(X=x)$ and class probability function ${\bm{\eta}}:{\mathcal{X}}\to\Delta_{k}$ (i.e., $\eta_{i}(x)={\mathbb{P}}(Y=i|X=x)$). We write ${\bm{\eta}}(x)$ as a vector ${\bm{\eta}}$ and omit $x$ when it is clear from context. Thus we can also represent the joint distribution as $D=(\bm{\pi},P_{1},\ldots,P_{k})$ (where ${\bm{\pi}}\in\Delta_{k}$) or $(M,{\bm{\eta}})$. For any $i\in[k]$, we assume $P_{i}$ has density $p_{i}$ with respect to the Lebesgue measure.

Remark on notations. To avoid confusion, we would like to emphasize that the class probability is denoted as $\eta_{i}(x)={\mathbb{P}}(Y=i|X=x)$ and the class conditional is denoted as $P_{i}(x)={\mathbb{P}}(X=x|Y=i)$ with density $p_{i}(x)$. The former further satisfies the normalization constraint: $\forall x\in{\mathcal{X}},\sum_{i=1}^{k}\eta_{i}(x)=1$, while $i$ in the latter one only serves as the index for $k$ different distributions.

In multi-class classification, given independent and identically distributed (i.i.d.) samples from the joint distribution $D(X,Y)$, we want to learn a probabilistic classifier $\hat{{\bm{\eta}}}:{\mathcal{X}}\to\Delta_{k}$ to approximate the true class probability function ${\bm{\eta}}$ by minimizing the following $\ell$-risk:

$${\mathcal{L}}_{\text{CPE}}(\hat{{\bm{\eta}}};D)=\mathbb{E}_{D(x,y)}[\ell(y,\hat{{\bm{\eta}}}(x))]=\mathbb{E}_{x\sim M}[\mathbb{E}_{y\sim{\bm{\eta}}(x)}[\ell(y,\hat{{\bm{\eta}}}(x))]]=\mathbb{E}_{x\sim M}[L({\bm{\eta}}(x),\hat{{\bm{\eta}}}(x))]$$

(1)

where $\ell:[k]\times\Delta_{k}\to{\mathbb{R}}$ is the loss function for using the class predictor $\hat{{\bm{\eta}}}(x)$ when the true class is $y$, and $L:\Delta_{k}\times\Delta_{k}\to{\mathbb{R}}$ is the expected loss of $\hat{{\bm{\eta}}}(x)$ under the true class probability ${\bm{\eta}}(x)$.

In statistical decision theory (Gneiting & Raftery, 2007), the negative proper loss is also called proper scoring rule (i.e., $S(y,\hat{{\bm{\eta}}}(x))=-\ell(y,\hat{{\bm{\eta}}}(x))$), which assesses the utility of the prediction. Properness of a loss is desirable in multi-class classification because it encourages the class probability estimator $\hat{{\bm{\eta}}}$ to match the true class probability function ${\bm{\eta}}$. An important property of proper loss is summarized in the following theorem:

Given the Bregman divergence representation of the point-wise regret in Theorem 1 and the $\ell$-risk in Equation (1), the excess risk of a class probability estimator $\hat{{\bm{\eta}}}$ over the Bayes optimal ${\bm{\eta}}$ is:

	$\displaystyle\mathrm{reg}(\hat{{\bm{\eta}}};M,{\bm{\eta}},\ell)\mathrel{\mathop{:}}=$	$\displaystyle{\mathcal{L}}_{\text{CPE}}(\hat{{\bm{\eta}}};D)-{\mathcal{L}}_{\text{CPE}}({\bm{\eta}};D)=\mathbb{E}_{M(x)}[L({\bm{\eta}}(x),\hat{{\bm{\eta}}}(x))-L({\bm{\eta}}(x),{\bm{\eta}}(x))]$		(4)
	$\displaystyle=$	$\displaystyle\mathbb{E}_{M(x)}[\mathbf{B}_{f}({\bm{\eta}}(x),\hat{{\bm{\eta}}}(x))]$

Csiszaŕ’s $f$-divergence is a popular way to measure the discrepancy between two probability distributions. Specifically, given two distributions $P,Q$ and a convex function $f:\mathbb{R}_{+}\to\mathbb{R}\cup\{\pm\infty\}$ satisfying $f(1)=0$, the $f$-divergence between $P$ and $Q$ is defined as $\mathbf{D}_{f}(P||Q)=\mathbb{E}_{Q}[f(\mathrm{d}P/\mathrm{d}Q)]$. In the following, we will introduce the multi-distribution extension of $f$-divergence (Garcia-Garcia & Williamson, 2012).

Inspired by the definition in Eq. (5), we consider the following canonical density ratio vector (more discussion about this choice can be found in Section 3.2): ${\bm{r}}(x)=(r_{1}(x),\ldots,r_{k}(x))$ where $r_{i}(x)\mathrel{\mathop{:}}=p_{i}(x)/p_{k}(x)$ and $r_{k}(x)=1$. Then we can connect a density ratio vector ${\bm{r}}(x)\in\mathbb{R}_{+}^{k-1}\times\{1\}$ and a class probability vector ${\bm{\eta}}(x)\in\Delta_{k}$ via an invertible link function.

According to Bayes’ theorem, we have:

$$\frac{{\mathbb{P}}(X=x,Y=i)}{{\mathbb{P}}(X=x,Y=k)}=\frac{\pi_{i}p_{i}(x)}{\pi_{k}p_{k}(x)}=\frac{M(x)\eta_{i}(x)}{M(x)\eta_{k}(x)}\Leftrightarrow r_{i}(x)=\frac{p_{i}(x)}{p_{k}(x)}=\frac{\pi_{k}}{\pi_{i}}\cdot\frac{\eta_{i}(x)}{\eta_{k}(x)}.$$

(6)

Thus we define the following multi-distribution link function $\Psi_{\mathrm{dr}}:\Delta^{k}\rightarrow\mathbb{R}_{+}^{k-1}\times\{1\}$ as a natural generalization of the binary DRE link function (Menon & Ong, 2016; Vernet et al., 2011):

$$[\Psi_{\mathrm{dr}}({\bm{\eta}}(x))]_{i}\mathrel{\mathop{:}}=\frac{\pi_{k}}{\pi_{i}}\cdot\frac{\eta_{i}(x)}{\eta_{k}(x)}=r_{i}(x),~{}\text{for all}~{}i\in[k].$$

(7)

Given Eq. (7) and the normalization constraint $\sum_{i\in[k]}\eta_{i}=1$, we obtain the inverse link function:

$$[\Psi^{-1}_{\mathrm{dr}}({\bm{r}}(x))]_{i}\mathrel{\mathop{:}}=\frac{\pi_{i}r_{i}(x)}{\sum_{j\in[k]}\pi_{j}r_{j}(x)}=\eta_{i}(x),~{}\text{for all}~{}i\in[k].$$

(8)

Thus given knowledge of the prior distribution ${\bm{\pi}}$ (which can also be easily estimated from empirical samples), one can transform a class probability estimator into a density ratio estimator via $\hat{{\bm{r}}}(x)=\Psi_{\mathrm{dr}}(\hat{{\bm{\eta}}}(x))$ and vice versa via $\hat{{\bm{\eta}}}(x)=\Psi^{-1}_{\mathrm{dr}}(\hat{{\bm{r}}}(x))$.

Following the basic formulation of multi-class experiments in Section 2.1, we now introduce the problem setup of multi-distribution density ratio estimation (DRE). Recall that ${\mathcal{X}}$ is the common data domain and $P_{1},\ldots,P_{k}$ are $k$ different distributions defined on ${\mathcal{X}}$ with densities $p_{1},\ldots,p_{k}$. Suppose we are given $n_{i}$ i.i.d. samples $\{x_{j}^{(i)}\}_{j=1}^{n_{i}}$ from each distribution $P_{i}$. The goal of multi-distribution DRE is to estimate the density ratios between all pairs of distributions $\{r_{ij}\mathrel{\mathop{:}}=p_{i}/p_{j}\}_{i,j\in[k]}$ from the i.i.d. datasets $\{\{x_{j}^{(i)}\}_{j=1}^{n_{i}}\}_{i=1}^{k}$. In this paper, we assume that the density ratios are always well-defined on domain ${\mathcal{X}}$ (e.g., when the distributions have strictly positive densities), which is also a common assumption in binary DRE problem (Kanamori et al., 2009; Kato & Teshima, 2021).

A naive approach towards this problem is to separately estimate each density $p_{i}$ from $\{x_{j}^{(i)}\}_{j=1}^{n_{i}}$ and then plug in $p_{i}$ and $p_{j}$ to get $r_{ij}$. However, as previous theoretical works (Kpotufe, 2017; Nguyen et al., 2007; Kanamori et al., 2012; Que & Belkin, 2013) suggest, directly estimating density ratios has many advantages in practical settings. Specifically, we know that (1) optimal convergence rates depend only on the smoothness of the density ratio and not on the densities; (2) optimal rates depend only on the intrinsic dimension of data, thus escaping the curse of dimension in density estimation. Inspired by these observations in binary DRE, this paper aims to develop a general framework for directly estimating multi-distribution density ratios. Moreover, we also theoretically prove that various interesting facts (Menon & Ong, 2016; Sugiyama et al., 2012), which hold in the binary case, extend to our multi-distribution case in Section 4.

While most previous works focus on DRE in binary cases, multi-distribution DRE has many important downstream applications. For example, given any integrable function $\phi:{\mathcal{X}}\to{\mathbb{R}}$, suppose we want to use importance sampling to estimate the expectation of $\phi$ with respect to a target distribution $Q$ with density $q$ w.r.t. the base measure:

$$\mathbb{E}_{q(x)}[\phi(x)]=\int_{\mathcal{X}}q(x)\phi(x)\mathrm{d}x=\int_{\mathcal{X}}p(x)\frac{q(x)}{p(x)}\phi(x)\mathrm{d}x=\mathbb{E}_{p(x)}\left[r(x)\cdot\phi(x)\right]$$

(9)

where we use the density ratio $r=p/q$ to correct the bias caused by using samples from the proposal distribution $p$ rather than the target distribution $q$. However, in practice, finding a good proposal is critical yet challenging (Owen & Zhou, 2000). An alternative and more robust strategy is to use a population of different proposals (sampling schemes) and use a set of density ratios to correct the bias, which is also known as multiple importance sampling (MIS) (Cappé et al., 2004; Elvira et al., 2015). Given $k$ different proposals $p_{1},\ldots,p_{k}$, the MIS estimation of the expectation is given by:

$$\mathbb{E}_{q(x)}[\phi(x)]=\sum_{i=1}^{k}\omega_{i}\mathbb{E}_{p_{i}(x)}\left[\frac{q(x)}{p_{i}(x)}\phi(x)\right]$$

(10)

where $\omega_{i}$ is the weight for each proposal $p_{i}$ and satisfies $\sum_{i}\omega_{i}=1$. Thus a more efficient and accurate multi-distribution DRE method will lead to better MIS. In the context of multi-source off-policy policy evaluation (Kallus et al., 2021), the proposals correspond to a set of demonstration policies and the target distribution is the query policy whose performance we want to evaluate from the offline multi-souce demonstrations; in the context of multi-domain transfer learning setting (covariate shift adaptation) (Bickel et al., 2008; Dinh et al., 2013), the proposals correspond to a set of data generating distributions (e.g. multiple source domains or various data augmentation strategies) and the target is the test distribution we care about. Estimating multi-distribution density ratios also allows us to compute important information quantities among multiple random variables such as the multi-distribution $f$-divergence in Equation (5), which can be used to analyze various kinds of discrepancy and correlations between multiple random variables and further has the potential of inspiring new generative models for multiple marginal matching problem (Cao et al., 2019).

Inspired by the success of Bregman divergence minimization for unifying various DRE methods in the binary case (Sugiyama et al., 2012), in this section, we propose a general framework for solving the multi-distribution density ratio estimation problem. First, we discuss our modeling choices. Although our goal is to estimate ${k\choose 2}$ density ratios (between all possible pairs), the solution set $\{r_{ij}\mathrel{\mathop{:}}=p_{i}/p_{j}\}_{i,j\in[k]}$ actually has $k-1$ degrees of freedom (e.g., $r_{ik}=r_{ij}\cdot r_{jk}$). Thus without loss of generality, we parametrize the following $k-1$ density ratio models $\hat{{\bm{r}}}_{\bm{\theta}}=(\hat{r}_{\theta_{1}},\ldots,\hat{r}_{\theta_{k-1}})$ to approximate the true canonical density ratios ${\bm{r}}=(r_{1},\ldots,r_{k-1})$, where $r_{i}\mathrel{\mathop{:}}=p_{i}/p_{k}$ for $i\in[k-1]$. For the simplicity of notation, we will omit the dependence on the parameters ${\bm{\theta}}$ and write our density ratio models as $\hat{{\bm{r}}}=(\hat{r}_{1},\ldots,\hat{r}_{k-1})$. An advantage of such modeling choice is that any density ratio can be recovered within one step of computation $\frac{p_{i}}{p_{j}}=\frac{p_{i}/p_{k}}{p_{j}/p_{k}}=\frac{r_{i}}{r_{j}}$, thus avoiding large compounding error while naturally ensuring consistency within the solution set (i.e., if we parametrize $\hat{r}_{ij}$, $\hat{r}_{jk}$ and $\hat{r}_{ik}$ respectively, we have to make sure they satisfy $\hat{r}_{ik}=\hat{r}_{ij}\cdot\hat{r}_{jk}$).

Since our goal is to optimize $\hat{{\bm{r}}}$ to approximate the true density ratios ${\bm{r}}$, we consider to use Bregman divergence (Def. 2) to measure the discrepancy between ${\bm{r}}$ and $\hat{{\bm{r}}}$. Specifically, for any strictly convex function $f:\mathbb{R}^{k-1}_{+}\to\mathbb{R}$ and $\forall x\in{\mathcal{X}}$, we have the following point-wise optimization problem:

$$\min_{\hat{{\bm{r}}}(x)\in{\mathbb{R}}^{k-1}_{+}}\mathbf{B}_{f}({\bm{r}}(x),\hat{{\bm{r}}}(x))=f({\bm{r}}(x))-f(\hat{{\bm{r}}}(x))-\langle\nabla f(\hat{{\bm{r}}}(x)),{\bm{r}}(x)-\hat{{\bm{r}}}(x)\rangle$$

(11)

which corresponds to the difference between the value of $f$ at ${\bm{r}}$, and the value of the first-order Taylor expansion of $f$ around point $\hat{{\bm{r}}}$ evaluated at point ${\bm{r}}$. Although the current formulation can be understood as a regression problem from $\hat{{\bm{r}}}(x)$ to the true density ratios ${\bm{r}}(x)$, we actually only have i.i.d. samples $x\sim p_{1},\ldots,p_{k}$ instead of the true targets ${\bm{r}}(x)$. In this case, we consider to use the following expected Bregman divergence to measure the overall discrepancy from the true density ratios ${\bm{r}}$ to the density ratio models $\hat{{\bm{r}}}$:

	$\displaystyle{\mathcal{L}}_{\text{DRE}}(\hat{{\bm{r}}};D)=$	$\displaystyle\int_{\mathcal{X}}p_{k}(x)\Big{(}f({\bm{r}}(x))-f(\hat{{\bm{r}}}(x))-\langle\nabla f(\hat{{\bm{r}}}(x)),{\bm{r}}(x)-\hat{{\bm{r}}}(x)\rangle\Big{)}\mathrm{d}x$		(12)
	$\displaystyle=$	$\displaystyle\mathbb{E}_{p_{k}(x)}\left[\langle\nabla f(\hat{{\bm{r}}}(x)),\hat{{\bm{r}}}(x)\rangle-f(\hat{{\bm{r}}}(x))\right]-\sum_{i\in[k-1]}\mathbb{E}_{p_{i}(x)}[\partial_{i}f(\hat{{\bm{r}}}(x))]+C$		(13)

where $C\mathrel{\mathop{:}}=\int_{\mathcal{X}}p_{k}(x)f({\bm{r}}(x))\mathrm{d}x=\mathbf{D}_{f}(P_{1},\ldots,P_{k-1}\|P_{k})$ is a constant with respect to $\hat{{\bm{r}}}$ and the equality comes from the fact that $p_{k}\cdot(r_{1},\ldots,r_{k-1})=(p_{1},\ldots,p_{k-1})$ according to the definition of ${\bm{r}}$. The rationale behind the above choice is that it allows us to get an unbiased estimation of the discrepancy between ${\bm{r}}$ and $\hat{{\bm{r}}}$ only using i.i.d. samples from $p_{1},\ldots,p_{k}$. Specifically, since $C$ is a constant, we have the following optimization problem over $\hat{{\bm{r}}}$ to approximate the true density ratios (where each expectation $\mathbb{E}_{p_{i}}$ can be empirically estimated using samples from $p_{i}$):

$$\min_{\hat{{\bm{r}}}:{\mathcal{X}}\to{\mathbb{R}}^{k-1}_{+}}\mathbb{E}_{p_{k}(x)}\left[\left\langle\nabla f(\hat{{\bm{r}}}(x)),\hat{{\bm{r}}}(x)\right\rangle-f(\hat{{\bm{r}}}(x))\right]-\sum_{i\in[k-1]}\mathbb{E}_{p_{i}(x)}\left[\partial_{i}f(\hat{{\bm{r}}}(x))\right]$$

(14)

Interestingly, the above multi-distribution DRE formulation, which is based on Bregman divergence minimization, can be alternatively derived from the perspective of variational estimation of multi-distribution $f$-divergence. In the following, We briefly discuss such an interpretation of Eq. (14).

Based on Fenchel duality, we can represent any strictly convex function $f:\mathbb{R}_{+}^{k-1}\to\mathbb{R}\cup\{+\infty\}$ through its conjugate function $f^{*}({\bm{s}})\mathrel{\mathop{:}}=\max_{{\bm{r}}\in{\mathbb{R}}^{k-1}_{+}}\langle{\bm{s}},{\bm{r}}\rangle-f({\bm{r}})$ as:

$$f({\bm{r}}(x))=\max_{{\bm{s}}:{\mathcal{X}}\to{\mathbb{R}}^{k-1}}\langle{\bm{r}}(x),{\bm{s}}(x)\rangle-f^{*}({\bm{s}}(x)),~{}~{}\text{for any}~{}x\in{\mathcal{X}}.$$

(15)

In order to estimate the multi-distribution $f$-divergence defined in Eq. (5) only using samples from $P_{1},\ldots,P_{k}$ (instead of their density information), we consider the following variational representation of multi-distribution $f$-divergence by substituting Eq. (15) into Eq. (5):

$\displaystyle\mathbf{D}_{f}(P_{1},\ldots,P_{k-1}||P_{k})=-\min_{{\bm{s}}:{\mathcal{X}}\rightarrow\mathbb{R}^{k-1}}\left[-\sum_{i\in[k-1]}\mathbb{E}_{p_{i}(x)}[{\bm{s}}(x)]_{i}+\mathbb{E}_{p_{k}(x)}f^{*}({\bm{s}}(x))\right]$

(16)

We then have the following lemma revealing the equivalence between the optimization problem in Eq. (14) and Eq. (16).

Multi-class Logistic Regression.  From Section 2.3, we know that there is a one-to-one correspondence between a class probability estimator and a density ratio estimator: $\hat{{\bm{r}}}=\Psi_{\mathrm{dr}}\circ\hat{{\bm{\eta}}}$ and $\hat{{\bm{\eta}}}=\Psi_{\mathrm{dr}}^{-1}\circ\hat{{\bm{r}}}$. For the clarity of presentation, here we assume the class prior distribution ${\bm{\pi}}$ is uniform such that $\hat{r}_{i}(x)=\hat{\eta}_{i}(x)/\hat{\eta}_{k}(x)$ and $\hat{\eta}_{i}(x)=\hat{r}_{i}(x)/\sum_{j=1}^{k}\hat{r}_{j}(x)$. To recover the loss of multi-class logistic regression, we choose the following convex function to be $f(\hat{r}_{1},\ldots,\hat{r}_{k-1})=\frac{1}{k}\sum_{i=1}^{k}\hat{r}_{i}\log\left(\hat{r}_{i}/\sum_{j=1}^{k}\hat{r}_{j}\right)$. In this case, the loss in Eq. (14) reduces to:

$$\frac{1}{k}\mathbb{E}_{p_{k}(x)}\left[\log\left(\sum_{j=1}^{k}\hat{r}_{j}(x)\right)\right]-\frac{1}{k}\sum_{i=1}^{k-1}\mathbb{E}_{p_{i}(x)}\left[\log\left(\frac{\hat{r}_{i}(x)}{\sum_{j=1}^{k}\hat{r}_{j}(x)}\right)\right]=-\left(\frac{1}{k}\sum_{i=1}^{k}\mathbb{E}_{p_{i}(x)}[\log\hat{\eta}_{i}(x)]\right)$$

(21)

We provide discussions for the general case (non-uniform prior ${\bm{\pi}}$) in Appendix A.4.1. Interestingly, we noticed that the above convex function also gives rise to the multi-distribution Jensen-Shannon divergence (Garcia-Garcia & Williamson, 2012) (also known as the information radius (Sibson, 1969), $\mathbf{D}_{f}(P_{1},\ldots,P_{k})=\frac{1}{k}\sum_{i=1}^{k}D_{\mathrm{KL}}(P_{i}\|\frac{1}{k}\sum_{j=1}^{k}P_{j})$) up to a constant of $\log k$.

Multi-LSIF.  When the convex function associated with the Bregman divergence is chosen to be $f(\hat{r}_{1},\ldots,\hat{r}_{k-1})=\frac{1}{2}\sum_{i=1}^{k-1}(\hat{r}_{i}-1)^{2}=\frac{1}{2}\|\hat{{\bm{r}}}-\bm{1}\|^{2}$, the loss in Eq. (14) reduces to:

$$\frac{1}{2}\sum_{i=1}^{k-1}\mathbb{E}_{p_{k}(x)}\left[\hat{r}_{i}^{2}(x)-1\right]-\sum_{i=1}^{k-1}\mathbb{E}_{p_{i}(x)}\left[\hat{r}_{i}(x)-1\right]=\frac{1}{2}\sum_{i=1}^{k-1}\mathbb{E}_{p_{k}(x)}\left[(\hat{r}_{i}(x)-r_{i}(x))^{2}\right]-C$$

(22)

where $C=\mathbb{E}_{p_{k}(x)}\left[\|{\bm{r}}(x)-1\|^{2}\right]$ is a constant w.r.t. $\hat{{\bm{r}}}$ and the minimizer to the above loss function matches the true density ratios, which strictly generalizes the Least-Squares Importance Fitting (LSIF) (Kanamori et al., 2009) method to the multi-distribution case.

Besides, we also consider the following simple convex functions that either strictly generalize their binary DRE counterparts as above, or lead to completely new methods for multi-distribution DRE:

• •

  Multi-KLIEP.  $f(\hat{r}_{1},\ldots,\hat{r}_{k-1})=\sum_{i=1}^{k-1}(\hat{r}_{i}\log\hat{r}_{i}-\hat{r}_{i})=\langle\hat{{\bm{r}}},\log(\hat{{\bm{r}}})\rangle-\|\hat{{\bm{r}}}\|_{1}$. This strictly generalizes the Kullback–Leibler Importance Estimation Procedure (KLIEP) (Sugiyama et al., 2008) to the multi-distribution case. See Appendix A.4.3 for more details.

• •

  Power.  $f(\hat{r}_{1},\ldots,\hat{r}_{k-1})=\sum_{i=1}^{k-1}\hat{r}_{i}^{\alpha}=\|\hat{{\bm{r}}}\|_{\alpha}^{\alpha}$,   for $\alpha>1$. This strictly generalizes the robust DRE method in (Sugiyama et al., 2012), which recovers Multi-KLIEP when $\alpha\to 1$ and Multi-LSIF when $\alpha=2$. See Appendix A.4.4 for more details.

• •

  Quadratic.  $f(\hat{r}_{1},\ldots,\hat{r}_{k-1})=\hat{{\bm{r}}}^{\top}{\bm{H}}\hat{{\bm{r}}}+{\bm{q}}^{\top}\hat{{\bm{r}}}$, for any positive definite matrix ${\bm{H}}\succ 0$. When ${\bm{H}}$ is the identity matrix and ${\bm{q}}=(-2,\ldots,-2)$, this is equivalent to Multi-LSIF.

• •

  LogSumExp.  $f(\hat{r}_{1},\ldots,\hat{r}_{k-1})=\alpha\log\left(\sum_{i=1}^{k-1}\exp(\hat{r}_{i}/\alpha)\right)$ for $\alpha>0$.

In principle, we can use any desired strictly convex function $f:{\mathbb{R}}^{k-1}_{+}\to{\mathbb{R}}$ within the optimization problem in Eq. (14), implying the great potential of our unified framework for discovering novel objectives for multi-distribution DRE. In terms of modeling flexibility, the curvature of different convex functions encode different inductive biases that may favor various downstream applications and we leave the design of more suitable convex functions for DRE as exciting future avenues.

From Section 4, we know that any strictly proper loss $\ell:[k]\times\Delta_{k}\to{\mathbb{R}}$ (or strictly proper scoring rule $S(i,\hat{{\bm{\eta}}})=-\ell(i,\hat{{\bm{\eta}}})$) in conjunction with the link function $\Psi_{\mathrm{dr}}$ also induces valid losses for multi-distribution DRE:

$$\min_{\hat{{\bm{r}}}:{\mathcal{X}}\to{\mathbb{R}}^{k-1}_{+}}\mathbb{E}_{D(x,y)}[\ell(y,\hat{{\bm{\eta}}}(x))]=\mathbb{E}_{x\sim M,y\sim{\bm{\eta}}(x)}[\ell(y,\Psi_{\mathrm{dr}}^{-1}(\hat{{\bm{r}}}(x)))]$$

(23)

In this work, we consider using the following classic proper scoring rules (Gneiting & Raftery, 2007), where $\hat{{\bm{\eta}}}$ is parametrized as $\Psi_{\mathrm{dr}}^{-1}(\hat{{\bm{r}}})$ (i.e. $\hat{\eta}_{i}=\pi_{i}\hat{r}_{i}/\sum_{j=1}^{k}\pi_{j}\hat{r}_{j}$):

• •

  Logarithm score. (Good, 1952) The loss is specified as $\ell(i,\hat{{\bm{\eta}}})=-\log(\hat{\eta}_{i})$, which also recovers the loss of multi-class logistic regression in Section 5.1.

• •

  Brier score. (Brier et al., 1950) The loss is specified as $\ell(i,\hat{{\bm{\eta}}})=-2\hat{\eta}_{i}+\sum_{j=1}^{k}\hat{\eta}_{j}^{2}+1$.

• •

  Logarithm pseudo-spherical score. (Good, 1971; Fujisawa & Eguchi, 2008) The loss is specified as $\ell(i,\hat{{\bm{\eta}}})=-\log\left(\frac{\hat{\eta}_{i}^{\alpha-1}}{(\sum_{j=1}^{k}\hat{\eta}_{j}^{\alpha})^{(\alpha-1)/\alpha}}\right)$ for $\alpha>1$.