# Doubly Calibrated Estimator for Recommendation
on Data Missing Not At Random

Wonbin Kweon  Pohang University of Science and TechnologyPohangRepublic of
Korea<EMAIL_ADDRESS>and  Hwanjo Yu  Pohang University of Science and
TechnologyPohangRepublic of Korea<EMAIL_ADDRESS>
(2024)

###### Abstract.

Recommender systems often suffer from selection bias as users tend to rate
their preferred items. The datasets collected under such conditions exhibit
entries missing not at random and thus are not randomized-controlled trials
representing the target population. To address this challenge, a doubly robust
estimator and its enhanced variants have been proposed as they ensure
unbiasedness when accurate imputed errors or predicted propensities are
provided. However, we argue that existing estimators rely on miscalibrated
imputed errors and propensity scores as they depend on rudimentary models for
estimation. We provide theoretical insights into how miscalibrated imputation
and propensity models may limit the effectiveness of doubly robust estimators
and validate our theorems using real-world datasets. On this basis, we propose
a Doubly Calibrated Estimator that involves the calibration of both the
imputation and propensity models. To achieve this, we introduce calibration
experts that consider different logit distributions across users. Moreover, we
devise a tri-level joint learning framework, allowing the simultaneous
optimization of calibration experts alongside prediction and imputation
models. Through extensive experiments on real-world datasets, we demonstrate
the superiority of the Doubly Calibrated Estimator in the context of debiased
recommendation tasks.

Recommender System; Collaborative Filtering; Personalization; Recommendation
Size; User Utility

††journalyear: 2024††copyright: rightsretained††conference: Proceedings of the
ACM Web Conference 2024; May 13–17, 2024; Singapore, Singapore††booktitle:
Proceedings of the ACM Web Conference 2024 (WWW ’24), May 13–17, 2024,
Singapore, Singapore††doi: 10.1145/3589334.3645417††isbn:
979-8-4007-0171-9/24/05††ccs: Information systems Collaborative
filtering††ccs: Information systems Personalization††ccs: Information systems
Recommender systems

##  1\. Introduction

Real-world recommender systems utilize feedback data collected from user-item
interactions for learning user behaviors. However, there has been recently a
growing concern regarding the issue of selection bias in user feedback
datasets. As users tend to rate their preferred items, the datasets collected
under such conditions exhibit entries missing not at random (MNAR) and thus
are not randomized-controlled trials representing the target population. To
tackle this issue, numerous debiasing methods have been proposed for unbiased
prediction of various user behaviors including explicit ratings (Schnabel et
al., 2016; Wang et al., 2019; Guo et al., 2021), implicit feedback (Saito et
al., 2020; Saito, 2020; Zhu et al., 2020; Chen et al., 2023), and post-click
conversion rate (Dai et al., 2022a; Li et al., 2023f, b; Zhang et al., 2020).
In the early literature, the error-imputation-based (EIB) estimator (Steck,
2010) and the inverse propensity scoring (IPS) estimator (Schnabel et al.,
2016) are two major approaches for designing unbiased estimators, capable of
accurately approximating the ideal loss over the target population with only
biased MNAR data.

Lately, the doubly robust (DR) estimator (Wang et al., 2019) and its enhanced
variants (Guo et al., 2021; Dai et al., 2022b; Li et al., 2023a, f, b) are
proposed to merge the EIB and the IPS estimator for double robustness. DR
estimators ensure the unbiased estimation of the target population when
accurate imputed errors or predicted propensities are provided. However,
rather than focusing on achieving precise estimation, they prioritize the
development of estimators that are robust to the inaccurate estimation of
imputed error and propensity score. They enhance the robustness of DR
estimators for better bias-variance trade-off by manipulating either the
imputed errors (Li et al., 2023b; Dai et al., 2022a; Guo et al., 2021) or the
propensity scores (Li et al., 2023f). While DR estimators exhibit state-of-
the-art debiasing performance based on their double robustness, we argue that
their effectiveness may be limited since they still depend on rudimentary
models (Koren et al., 2009; He et al., 2017) for estimating imputed errors and
propensity scores.

Recent work in the field of machine learning has highlighted that both
logistic regression and neural networks, commonly adopted for imputation and
propensity models, have a tendency to generate overly confident predictions
(Guo et al., 2017; Kull et al., 2017; Kweon et al., 2022; Bai et al., 2021).
The overconfidence problem can be exacerbated especially for DR estimators, as
user-item pairs for training the prediction model and the propensity model
overlap with each other. DR estimators adopt inverse propensity scoring for
the observed pairs that already gave positive signals to the propensity model.
Indeed, in our analysis with real-world datasets, we demonstrate that the
current imputed error and propensity score estimation yields miscalibrated
estimates that fail to accurately reflect the ground-truth likelihood.
Miscalibrated imputed errors can increase the bias and the variance of DR
estimators and overconfident propensity scores can yield a variance
amplification problem. Nonetheless, literature has yet to delve into the
direct methodology for training accurate imputation and propensity models,
highlighting the necessity for such exploration.

This paper claims that DR estimators can be further enhanced by leveraging
model calibration approaches for the imputation and the propensity models. We
first provide theoretical insights into how miscalibrated imputation and
propensity models may limit the effectiveness of doubly robust estimators: the
bias and the variance have an upper bound that is proportional to the
calibration errors. On this basis, we propose Doubly Calibrated Estimator that
involves the calibration of both the imputation and propensity models. To
achieve this, we introduce calibration experts, each assigned to a specific
group of users via the assignment network. By doing so, each expert can learn
specialized knowledge about its group for the calibration of imputation and
propensity models and also get enough training signals from users within the
group. Additionally, we devise a tri-level joint learning framework, allowing
the simultaneous optimization of calibration experts alongside prediction and
imputation models.

The proposed method offers several merits as follows: (1) Our approach is
orthogonal to existing DR estimators and can be seamlessly combined with them,
(2) It enables the simultaneous reduction of both bias and variance of DR
estimators, (3) It does not require any additional unbiased data. The main
contributions of this paper are summarized as follows:

  * •

We provide a theoretical analysis on the correlation between the performance
of DR estimators and the calibration of the imputation and the propensity
models. Then, we demonstrate existing DR estimators may exhibit limited
effectiveness, as they rely on miscalibrated imputed errors and propensity
scores.

  * •

We propose Doubly Calibrated Estimator that involves the calibration of both
the imputation and propensity models. Imputed errors and propensity scores are
calibrated with calibration experts by assigning users to each expert through
the assignment network.

  * •

We validate the superiority of the proposed method by extensive experiments on
real-world datasets. We also provide in-depth quantitative analyses to verify
the effectiveness of each proposed component.

##  2\. Preliminaries

###  2.1. Problem Formulation

Let $\mathcal{U}=\\{u\\}$ and $\mathcal{I}=\\{i\\}$ denote a set of users and
a set of items, respectively. For a pair of $u\in\mathcal{U}$ and
$i\in\mathcal{I}$, an observation indicator $o_{u,i}$ is given as 1 if the
item $i$ is exposed to the user $u$ and 0 otherwise. The observation indicator
$o_{u,i}$ is often referred to as treatment since the rating $r_{u,i}$ is
observed only when an item is exposed to a user ($o_{u,i}=1$) (Li et al.,
2023f, b; Dai et al., 2022a). In this paper, we adopt the binary rating scheme
($r_{u,i}\in\\{0,1\\}$) for illustrating purposes, following recent literature
in the post-click conversion rate scenario (Dai et al., 2022a; Li et al.,
2023f, b; Zhang et al., 2020) or the implicit feedback scenario (Saito et al.,
2020; Saito, 2020; Zhu et al., 2020). To distinguish the entire population and
observed distribution, let $\mathcal{D}=\mathcal{U}\times\mathcal{I}$ denote
the set of all user-item pairs and
$\mathcal{O}=\\{(u,i)|(u,i)\in\mathcal{D},o_{u,i}=1\\}$ denote the set of
observed pairs.

If ratings for all user-item pairs in $\mathcal{D}$ are observed, the
prediction model $\hat{r}_{u,i}=f_{\theta}(x_{u,i})$ with input feature
$x_{u,i}$ can be trained by the ideal loss function on full observation:

(1)
$\mathcal{L}_{\textup{ideal}}(\theta)=\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}e_{u,i},$

where $e_{u,i}=e(\hat{r}_{u,i},r_{u,i})$ represents the prediction error
between the prediction $\hat{r}_{u,i}$ and the target $r_{u,i}$, and can be
Mean Squared Error (MSE) (Koren et al., 2009; Chai and Draxler, 2014) or
Binary Cross Entropy (BCE) (He et al., 2017). Instead, the prediction model is
often trained by a naive estimator on the only observed ratings, since most
ratings are missing due to the users’ selection mechanism:

(2)
$\mathcal{E}_{\textup{naive}}(\theta)=\frac{1}{|\mathcal{O}|}\sum_{u,i\in\mathcal{O}}e_{u,i}.$

For this naive estimator to be an unbiased estimator of the ideal loss, we
need the expected value of its estimation across all the possible observations
$\mathcal{O}$ to precisely match the ideal loss (Wang et al., 2019), i.e.,
$\mathbb{E}_{\mathcal{O}}[\mathcal{E}_{\textup{naive}}]=\mathcal{L}_{\textup{ideal}}$.
However, as addressed in (Schnabel et al., 2016), the users in the recommender
systems tend to rate their preferred items and this process is missing not at
random (MNAR). Therefore, the observed set $\mathcal{O}$ is not a result of
randomized-controlled trials and the naive estimator may induce a large bias
(Schnabel et al., 2016):
$|\mathbb{E}_{\mathcal{O}}[\mathcal{E}_{\textup{naive}}]-\mathcal{L}_{\textup{ideal}}|>0$.

###  2.2. Unbiased Estimators

There have been proposed two major approaches for designing the unbiased
estimators. First, Error-imputation-based (EIB) estimators (Steck, 2010; Wang
et al., 2020) adopt imputation models to directly predict the individual loss
$e_{u,i}$ for each pair. Then, EIB estimators use the imputed error as a proxy
of errors for missing ratings:

(3)
$\mathcal{E}_{\textup{EIB}}(\theta)=\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}(o_{u,i}e_{u,i}+(1-o_{u,i})\hat{e}_{u,i}),$

where $e_{u,i}=e(\hat{r}_{u,i},r_{u,i})$ is the true error for observed
ratings and $\hat{e}_{u,i}=e(\hat{r}_{u,i},\tilde{r}_{u,i})$ is the imputed
error predicted by the imputation model $g_{\phi}(x_{u,i})=\tilde{r}_{u,i}$.
Obviously, the EIB estimator is an unbiased estimator of the ideal loss when
the imputed errors are accurate, i.e, $(\hat{e}_{u,i}=e_{u,i}\ \
\forall(u,i)\in\mathcal{D})\Rightarrow(\mathbb{E}_{\mathcal{O}}[\mathcal{E}_{\textup{EIB}}]=\mathcal{L}_{\textup{ideal}})$.

On the other hand, Inverse Propensity Scoring (IPS) estimators (Schnabel et
al., 2016; Saito et al., 2020) adopt propensity models to predict the
probability of observing the true rating. Then, IPS estimators use the
propensity score to inversely weight the prediction error for observed
ratings:

(4)
$\mathcal{E}_{\textup{IPS}}(\theta)=\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}\frac{o_{u,i}e_{u,i}}{\hat{p}_{u,i}},$

where $\hat{p}_{u,i}=\hat{P}(o_{u,i}=1)$ is the propensity score predicted by
the propensity model $h_{\psi}(x_{u,i})=\hat{p}_{u,i}$. The IPS estimator is
an unbiased estimator of the ideal loss when the predicted propensity scores
are accurate, i.e., $(\hat{p}_{u,i}=p_{u,i}\ \
\forall(u,i)\in\mathcal{O})\Rightarrow(\mathbb{E}_{\mathcal{O}}[\mathcal{E}_{\textup{IPS}}]=\mathcal{L}_{\textup{ideal}})$.

###  2.3. Doubly Robust Estimators

Recently, Doubly Robust (DR) estimator (Wang et al., 2019) and its enhanced
variants (Li et al., 2023a, f, b; Dai et al., 2022b) are proposed to merge the
EIB and the IPS estimator for double robustness. Given imputed error
$\hat{e}=e(\hat{r},\tilde{r})$ and learned propensity score
$\hat{p}=h_{\psi}(x_{u,i})$, DR estimator can be formulated as follows:

(5)
$\mathcal{E}_{\textup{DR}}(\theta)=\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}\Big{(}\hat{e}_{u,i}+\frac{o_{u,i}(e_{u,i}-\hat{e}_{u,i})}{\hat{p}_{u,i}}\Big{)}.$

By utilizing both the imputed error and the propensity score, DR estimators
have double robustness: DR estimator is an unbiased estimator when either the
imputed errors or the propensity scores are accurate. The following lemma
presents the bias and the variance of DR estimator induced by the inaccurate
estimation of imputed errors and propensity scores. Please refer (Wang et al.,
2019; Dai et al., 2022a) for the proofs.

###### Lemma 1 (Bias and Variance of DR estimator).

The bias and variance of DR estimator are computed as:

(6)
$\displaystyle\textup{Bias}[\mathcal{E}_{\textup{DR}}]=\frac{1}{|\mathcal{D}|}\Bigl{|}\sum_{u,i\in\mathcal{D}}\Big{(}\frac{\hat{p}_{u,i}-p_{u,i}}{\hat{p}_{u,i}}\Big{)}(e_{u,i}-\hat{e}_{u,i})\Bigr{|},$
$\displaystyle\textup{Var}[\mathcal{E}_{\textup{DR}}]=\frac{1}{|\mathcal{D}|^{2}}\sum_{u,i\in\mathcal{D}}\frac{p_{u,i}(1-p_{u,i})}{\hat{p}^{2}_{u,i}}(\hat{e}_{u,i}-e_{u,i})^{2}.$

###  2.4. Limitations of Existing Estimators

According to Lemma 1, the bias and the variance of DR estimator are positively
correlated with the inaccuracy of the imputation model and the propensity
model. Nevertheless, both models have been learned in a surprisingly
simplistic manner. In the early methods (Steck, 2010; Saito et al., 2020),
heuristic techniques are adopted for the estimation of imputed errors and
propensity scores, e.g., $\hat{e}_{u,i}=\omega|\hat{r}_{u,i}-\gamma|$, where
$\omega$ and $\gamma$ are hyper-parameters. Recent studies have also embraced
straightforward model-based approaches, such as employing logistic regression
on held-out unbiased data (Schnabel et al., 2016; Li et al., 2023f) or
training binary classifiers for the observation indicator (Wang et al., 2019;
Li et al., 2023b; Dai et al., 2022a). These model-based approaches exhibit a
slight improvement over heuristic techniques, resulting in a moderately
accurate prediction model. Nevertheless, as indicated by recent literature
(Guo et al., 2017; Kweon et al., 2022), both logistic regression and neural
networks have a tendency to generate overly confident predictions. An over-
confident propensity model can lead to propensity scores that are either too
low or too high, and these poorly calibrated estimates may further hinder the
effectiveness of DR estimator.

On the other hand, several recent works (Guo et al., 2021; Li et al., 2023a;
Dai et al., 2022a) focus on developing estimators that are robust to the
inaccurate estimation of imputed errors and propensity scores. They enhance
the robustness of DR estimator by reducing the bias and the variance while
having the same level of inaccuracy in imputed errors and propensity scores.
Nevertheless, we argue that their effectiveness may be limited since they
still depend on rudimentary models for estimating imputed errors and
propensity scores. The literature has yet to delve into the direct methodology
for training accurate imputation and propensity models without any unbiased
data, highlighting the necessity for such exploration.

###  2.5. Model Calibration

In this paper, we adopt the concept of model calibration (Guo et al., 2017) to
quantitatively measure the inaccuracy of the imputation and the propensity
models. We say a model is calibrated if its output reflects the ground-truth
likelihood of correctness (Kull et al., 2017). For the propensity model
$h_{\psi}$ and the observation indicator $o$, a formal definition can be
formulated as follows:

(7) $\mathbb{E}[o|h_{\psi}(x)=\hat{p}]=\hat{p}\quad\forall\hat{p}\in[0,1].$

For example, if we have 100 pairs with propensity scores $\hat{p}_{u,i}=0.2$,
we expect exactly of these pairs to be observed ($o_{u,i}=1$). Using the above
definition, the miscalibration of a propensity model can be measured by
Expected Calibration Error (ECE) and Maximum Calibration Error (MCE) (Naeini
et al., 2015):

(8)
$\displaystyle\text{ECE}(h_{\psi})=\mathbb{E}_{\hat{p}}\big{[}\lvert\mathbb{E}[o|h_{\psi}(x)=\hat{p}]-\hat{p}\rvert\big{]},$
$\displaystyle\text{MCE}(h_{\psi})=\max_{\hat{p}}\lvert\mathbb{E}[o|h_{\psi}(x)=\hat{p}]-\hat{p}\rvert.$

ECE and MCE quantify the average and worst-case discrepancy between actual
observation proportion and average predicted propensity across $M$ bins.
Likewise, we can compute ECE and MCE of the imputation model $g_{\phi}$ by
substituting $(h_{\psi},o,\hat{p})$ with $(g_{\phi},r,\tilde{r})$ in Eq.8.

##  3\. Calibration and DR Estimators

In this section, we describe the motivation for proposing our doubly
calibrated estimator. We first theoretically analyze how DR estimator stands
to benefit from the calibration of the imputation and the propensity models.
Then, we empirically demonstrate that existing imputation and propensity
models are indeed miscalibrated and limit the effectiveness of DR estimator.

###  3.1. Theoretical Analysis

We analyze how the miscalibration of the imputation and the propensity models
amplify the bias and the variance of DR estimator and further hinder the
effectiveness of DR estimators.

####  3.1.1. Bias of DR Estimator

We first present a theorem concerning the interplay between the bias of DR
estimator and the calibration of the propensity model.

###### Theorem 1.

The bias of DR estimator exhibits an upper bound proportional to the
calibration error of the propensity model.

$\displaystyle\textup{Bias}[\mathcal{E}_{\textup{DR}}]\leq\rho_{\textup{max}}\cdot\textup{ECE}(h_{\psi}),$
$\displaystyle\textup{Bias}[\mathcal{E}_{\textup{DR}}]\leq\rho_{\textup{max}}\cdot\textup{MCE}(h_{\psi}),$

where
$\rho_{\textup{max}}=\textup{max}_{(u,i)\in\mathcal{D}}\lvert(e_{u,i}-\hat{e}_{u,i})/\hat{p}_{u,i}\rvert$.

Please refer to Appendix A for the proof. The above theorem implies that DR
estimator may become unreliable when the propensity models are over-confident
and yield large calibration errors. Similarly, we can derive the following
corollary of Theorem 1 for the imputation model.

###### Corollary 2.

The bias of DR estimator exhibits an upper bound proportional to the
calibration error of the imputation model.

$\displaystyle\textup{Bias}[\mathcal{E}_{\textup{DR}}]\leq\pi_{\textup{max}}\cdot\textup{ECE}(g_{\phi}),$
$\displaystyle\textup{Bias}[\mathcal{E}_{\textup{DR}}]\leq\pi_{\textup{max}}\cdot\textup{MCE}(g_{\phi}),$

where
$\pi_{\textup{max}}=\textup{max}_{(u,i)\in\mathcal{D}}\lvert(\hat{p}_{u,i}-p_{u,i})(e^{(1)}_{u,i}-e^{(0)}_{u,i})/\hat{p}_{u,i}\rvert$.

$e^{(r)}$ denotes the loss when the target label is $r$, e.g.,
$e^{(r)}=-r\text{log}\hat{r}-(1-r)\text{log}(1-\hat{r})$ for BCE. Please refer
to Appendix A for the proof. Likewise, the above corollary implies that
calibrated imputed errors can reduce the upper bound on the bias of the DR
estimator, leading toward the unbiased estimation of the ideal loss.

####  3.1.2. Variance of DR Estimator

Moreover, we argue that the miscalibration of the imputation model also yields
negative effects for the variance of the DR estimator by the following
theorem.

###### Theorem 3.

The variance of DR estimator exhibits an upper bound proportional to the
square of the calibration error of the imputation model.

$\displaystyle\textup{Var}[\mathcal{E}_{\textup{DR}}]\leq\omega_{\textup{max}}\cdot\big{(}\textup{ECE}(g_{\phi})\big{)}^{2},$
$\displaystyle\textup{Var}[\mathcal{E}_{\textup{DR}}]\leq\frac{\omega_{\textup{max}}}{|\mathcal{D}|}\cdot\big{(}\textup{MCE}(g_{\phi})\big{)}^{2},$

where $\omega_{\textup{max}}=\textup{max}_{(u,i)\in\mathcal{D}}\lvert
p_{u,i}(1-p_{u,i})(e^{(1)}_{u,i}-e^{(0)}_{u,i})^{2}/\hat{p}^{2}_{u,i}\rvert$.

The above theorem implies that the calibration of the imputation model is also
able to reduce the upper bound on the variance of the DR estimator. We expect
this theorem to give us a new remedy for bias-variance trade-off in the DR
estimator as calibrating the pseudo label of the imputation model can reduce
both the bias and the variance of the DR estimator simultaneously.

On the other hand, as shown in Lemma 1, the accuracy of the propensity model
does not have any relation with the variance of the DR estimator.
Nevertheless, the variance is inversely proportional to $\hat{p}^{2}_{u,i}$,
which can be problematic when the estimated propensity scores are
exceptionally low. Indeed, current propensity models often produce too low
propensity scores as they can easily be over-confident (Guo et al., 2017;
Kweon et al., 2022), i.e., over-confidence for non-observation. While existing
methods (Saito et al., 2020; Li et al., 2023f, b) incorporate a propensity
clipping technique (Su et al., 2019) to curtail extremely low propensities,
this approach does not address the underlying issue of overconfidence.
Instead, we argue that model calibration can effectively mitigate the variance
amplification problem, as a calibrated propensity model is less likely to
produce extremely low propensities compared to a miscalibrated one.

###  3.2. Empirical Evidence

In this section, we demonstrate that the imputation and propensity models
deployed in recent methods are miscalibrated and thus limit the effectiveness
of DR estimators, which can be explained by our theoretical analysis. We adopt
two DR estimators including DR-JL (Wang et al., 2019) and TDR (Li et al.,
2023b), and two kinds of imputation and propensity models, including simple
heuristics and model-based approaches. (1) For the heuristic approach, we use
$\hat{p}_{u,i}=(\sum_{u\in\mathcal{U}}Y_{u,i}/\text{max}_{i\in\mathcal{I}}\sum_{u\in\mathcal{U}}Y_{u,i})^{0.5}$
(Saito et al., 2020) for the propensity scores and
$\hat{e}_{u,i}=\omega|\hat{r}_{u,i}-\gamma|$ (Steck, 2010) for the imputed
errors ($\omega$ and $\gamma$ are hyper-parameters). (2) For the model-based
approach, we adopt the neural collaborative filtering (He et al., 2017) for
the imputation and the propensity models. We train the imputation model
$g_{\phi}(x_{u,i})=\tilde{r}_{u,i}$ by using the imputation loss adopted in
(Wang et al., 2019; Li et al., 2023f):

(9)
$\mathcal{L}_{\textup{imp}}(\phi)=\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}\frac{o_{u,i}(e(\hat{r},\tilde{r})-e(\hat{r},r))^{2}}{\hat{p}_{u,i}},$

and train the propensity model $h_{\psi}(x_{u,i})=\hat{p}_{u,i}$ through
binary classification between $\mathcal{O}$ and
$\mathcal{D}\setminus\mathcal{O}$, as done in (Li et al., 2023b).

We investigate the miscalibration of the adopted imputation and propensity
models and the performance of DR estimators when they deploy each of the above
imputation and propensity models. For the quantitative measurement of the
miscalibration, we adopt ECE defined in Eq.8. However, since we cannot observe
the true propensity
$\mathbb{E}[o|h_{\psi}(x)=\hat{p}]=P(o=1|h_{\psi}(x)=\hat{p})$, we cannot
directly compute Eq.8. Instead, we partition the [0,1] range of $\hat{p}$ into
$M$ bins and aggregate the value of pairs in each bin as done in (Naeini et
al., 2015; Guo et al., 2017):

(10) $\text{ECE}_{M}(h_{\psi})=\sum_{m=1}^{M}\frac{\lvert
B_{m}\rvert}{N}\left|\frac{\sum_{(u,i)\in B_{m}}o_{u,i}}{\lvert
B_{m}\rvert}-\frac{\sum_{(u,i)\in B_{m}}\hat{p}_{u,i}}{\lvert
B_{m}\rvert}\right|$

where $B_{m}$ is $m$-th bin and $N$ is the number of samples. The first term
in the absolute value symbols denotes the ground-truth proportion of
observations in $B_{m}$ and the second term denotes the average predicted
propensity score of $B_{m}$. We use $M=15$ as done in (Guo et al., 2017; Kweon
et al., 2022).

Figure 1. Performance of DR estimators with various imputation/propensity
models on Yahoo!R3 dataset.

Figure 1 shows the performance of DR estimators with various imputation and
propensity models on Yahoo!R3
dataset111http://research.yahoo.com/Academic_Relations. TDR_Cal denotes TDR
with our calibrated imputation model and calibrated propensity model (we show
this point as a reference). Please note that figures for $\text{MSE vs
ECE}(h_{\psi})$ and $\text{AUC vs ECE}(g_{\phi})$ can be readily reconstructed
by Figure 1. We have the following findings: (1) Existing approaches for
estimating imputed errors and propensity scores produce poorly calibrated
probabilities. This finding is consistent with the previous work (Guo et al.,
2017; Kweon et al., 2022) concerning that the output of machine learning
models, from logistic regression to neural networks, does not necessarily
indicate the accurate correctness likelihood. (2) The performance of DR
estimators is positively correlated with the calibration of the imputation and
the propensity models. We observe that even a naive DR estimator with the
model-based approach outperforms TDR with the heuristic approach. This
experimental result serves as evidence to support our theoretical analysis in
Section 3.1.

##  4\. Doubly Calibrated Estimator

We propose a Doubly Calibrated Estimator on the basis of our comprehensive
analysis on the correlation between the performance of DR estimators and the
calibration of the imputation and the propensity models. A doubly calibrated
estimator entails the calibration of both the imputation and propensity
models, achieved through the calibration experts (Section 4.2). These experts
are simultaneously optimized alongside prediction and imputation models within
our tri-level joint learning framework (Section 4.3).

###  4.1. Motivation

As elaborated upon in Section 3, current imputation and propensity models are
miscalibrated, which, as substantiated by our theoretical findings, could
potentially constrain the effectiveness of DR estimators. Therefore, DR
estimators can benefit from the calibration of imputation and propensity
models. The straightforward approach for the calibration is to adopt a post-
processing calibration function $c_{\omega}:[0,1]\rightarrow[0,1]$ that maps
the miscalibrated $\tilde{r}_{u,i}=g_{\phi}(x_{u,i})$ or
$\hat{p}_{u,i}=f_{\psi}(x_{u,i})$ to the well-calibrated probability. Among
various forms of calibration functions, we adopt Platt scaling
$c_{\omega}(p)=\sigma(a\cdot\sigma^{-1}(p)+b)$ (Platt et al., 1999), a general
form of the temperature scaling (Guo et al., 2017). Here, $\sigma^{-1}$ is the
inverse of the sigmoid function and $\sigma^{-1}(p)$ represents the logit of
the probability $p$. Platt scaling has found widespread application across
diverse domains, including computer vision (Minderer et al., 2021; Ding et
al., 2021), natural language processing (Desai and Durrett, 2020), and
recommender system (Kweon et al., 2022, 2024).

However, employing a single global calibration function for all users falls
short of achieving satisfactory performance. The learning parameters
$\omega=\\{a,b\\}$ represent the characteristics of the logit distribution
(Kull et al., 2017), e.g., $a=\mu_{1}/\sigma_{1}^{2}-\mu_{0}/\sigma_{0}^{2}$,
where $\mu$ and $\sigma$ is the mean and the variance of the logit
distribution. As a result, a global calibration function would blend
information from users with varying preferences, unable to fully capture the
distinct logit distribution of individual users. A naive solution to this
challenge is to create a dedicated calibration function for each user:
$c_{\omega_{u}}(p)=\sigma(a_{u}\cdot\sigma^{-1}(p)+b_{u}),\ \forall
u\in\mathcal{U}$. Here, the user-specific parameters
$\omega_{u}=\\{a_{u},b_{u}\\}$ account for the unique logit distributions of
user $u$. Nevertheless, this approach also has a limitation. Each calibration
function requires training using the respective user’s interactions and
therefore cold-start users with minimal interactions may not possess
sufficient training signals for their calibration function. In the following
subsection, we present our solution to strike a balance between the single
global calibration and the user-specific calibration.

###  4.2. Calibration Experts

We introduce calibration experts that consider distinct logit distributions
across users while alleviating the cold-start problem of user-specific
calibration. Inspired by Mixture-of-Experts (Jacobs et al., 1991; Kang et al.,
2020), we deploy $K$ calibration experts, denoted as
$\\{c_{\omega_{k}}(\cdot)\\}_{k\in[K]}$,222$[K]$ denotes $\\{1,...,K\\}$ with
each expert assigned to a specific group of users. By doing so, each expert
can learn specialized knowledge about its group for the calibration of
imputation and propensity models and also get enough training signals from the
users in the group.

Specifically, we devise an assignment network to assign each user to an
expert, by mapping the user embedding to the assignment probability:

(11) $A(E_{u})=\alpha_{u}\in\mathbb{R}^{K},$

where $A:\mathbb{R}^{d}\rightarrow\mathbb{R}^{K}$ is the assignment network
with a softmax output layer. $E_{u}\in\mathbb{R}^{d}$ is the user embedding
and $\alpha_{u}\in\mathbb{R}^{K}$ is the assignment probability vector of $u$.
Each element of the assignment probability vector $\alpha_{u,k}$ represents
the probability of the user $u$ being assigned to expert
$c_{\omega_{k}}(\cdot)$. Then, we sample an assignment vector from the
assignment probability vector.

(12) $\beta_{u}\sim\text{Categorical}(\alpha_{u}),$

where $\beta_{u}$ is a $K$-dimensional one-hot assignment vector sampled from
the categorical distribution with the probability
$\\{\alpha_{u,k}\\}_{k\in[K]}$. However, the sampling process is non-
differentiable and blocks the gradient flow. To tackle this challenge, we
adopt Gumbel-Softmax (Jang et al., 2016), a continuous relaxation with the
reparameterization trick. With Gumbel-Softmax, the assignment vector
$\beta_{u}$ can be pseudo-sampled as follows:

(13)
$\beta_{u,k}=\frac{\textrm{exp}((\textrm{log}\alpha_{u,k}+g_{k})/\tau)}{\Sigma_{l=1}^{K}\textrm{exp}((\textrm{log}\alpha_{u,l}+g_{l})/\tau)}\quad\forall
k\in[K],$

where $g$ is i.i.d drawn from Gumbel distribution with ($\text{location}=0$,
$\text{scale}=1$), and $\tau$ is the temperature of the softmax. We use a
simple exponential annealing $\tau=T_{0}(\frac{T_{Q}}{T_{0}})^{\frac{q}{Q}}$,
where $q$ is the current epoch, $Q$ is the total epochs, $T_{0}$ is initial
temperature and $T_{Q}$ is the terminal temperature $(T_{0}>>T_{Q})$. With a
large $\tau$ the assignment network explores the combinations of the
calibration experts, and with a small $\tau$ the assignment network can select
a specific calibration expert for a user.

It is noted that we distinguish the set of calibration experts for imputation
and propensity models. We represent
$\\{c^{\textup{imp}}_{\omega_{k}}(\cdot)\\}_{k\in[K]}$ as calibration experts
for the imputation model and
$\\{c^{\textup{prop}}_{\omega_{k}}(\cdot)\\}_{k\in[K]}$ as calibration experts
for the propensity model. Accordingly, the assignment network for the
imputation model $A^{\textup{imp}}$ takes the embedding of the imputation
model $g_{\phi}$ as an input and the assignment network for the propensity
model $A^{\textup{prop}}$ takes the embedding of the propensity model
$h_{\psi}$ as an input. Calibration experts and assignment networks are
trained in an end-to-end manner as follows: (1) For the imputation model,
$\\{c^{\textup{imp}}_{\omega_{k}}(\cdot)\\}_{k\in[K]}$ and $A^{\textup{imp}}$
are trained by the loss

(14)
$\displaystyle\mathcal{L}_{\textup{impcal}}(c^{\textup{imp}}_{\omega_{k}},A^{\textup{imp}})$
$\displaystyle=\sum_{(u,i)\in\mathcal{O}_{\textup{val}}}-\frac{r_{u,i}}{\bar{p}_{u,i}}\text{log}(\sum_{k=1}^{K}\beta^{\textup{imp}}_{u,k}\cdot
c^{\textup{imp}}_{\omega_{k}}(\tilde{r}))$
$\displaystyle-(1-\frac{r_{u,i}}{\bar{p}_{u,i}})\text{log}(1-\sum_{k=1}^{K}\beta^{\textup{imp}}_{u,k}\cdot
c^{\textup{imp}}_{\omega_{k}}(\tilde{r})),$

where $\tilde{r}_{u,i}=g_{\phi}(x_{u,i})$ is the miscalibrated pseudo label of
the imputation model and
$\bar{p}_{u,i}=\sum_{k=1}^{K}\beta^{\textup{prop}}_{u,k}\cdot
c^{\textup{prop}}_{\omega_{k}}(\hat{p})$ is the calibrated propensity score.
$\mathcal{O}_{\textup{val}}$ is the set of observed pairs in the validation
set. It is a common practice to adopt validation set as the calibration set
(Guo et al., 2017; Kull et al., 2017; Ding et al., 2021; Kweon et al., 2022).
We use the binary cross-entropy loss, a widely adopted loss function for the
calibration of binary classifiers (Guo et al., 2017; Kweon et al., 2022; Kull
et al., 2017). (2) For the propensity model, similarly,
$\\{c^{\textup{prop}}_{\omega_{k}}(\cdot)\\}_{k\in[K]}$ and
$A^{\textup{prop}}$ are trained by the loss

(15)
$\displaystyle\mathcal{L}_{\textup{propcal}}(c^{\textup{prop}}_{\omega_{k}},A^{\textup{prop}})$
$\displaystyle=\sum_{(u,i,o)\in\mathcal{D}_{\textup{val}}}-o_{u,i}\text{log}(\sum_{k=1}^{K}\beta^{\textup{prop}}_{u,k}\cdot
c^{\textup{prop}}_{\omega_{k}}(\hat{p}))$
$\displaystyle-(1-o_{u,i})\text{log}(1-\sum_{k=1}^{K}\beta^{\textup{prop}}_{u,k}\cdot
c^{\textup{prop}}_{\omega_{k}}(\hat{p})),$

where $\hat{p}_{u,i}=h_{\psi}(x_{u,i})$ is the miscalibrated propensity score
and
$\mathcal{D}_{\textup{val}}=\\{(u,i,o_{u,i}=1)|(u,i)\in\mathcal{O}_{\textup{val}}\\}\cup\\{(u,i,o_{u,i}=0)|(u,i)\in\mathcal{D}\setminus(\mathcal{O}_{\textup{val}}\cup\mathcal{O})\\}$.
During the training of the calibration experts, the imputation model
$g_{\phi}(x_{u,i})=\tilde{r}_{u,i}$ and the propensity model
$h_{\psi}(x_{u,i})=\hat{p}_{u,i}$ are frozen.

###  4.3. Tri-Level Joint Learning Framework

We introduce a tri-level joint learning framework where the proposed
calibration experts can be simultaneously optimized along with the existing DR
estimators. We have four core components including prediction model
$f_{\theta}(x_{u,i})=\hat{r}_{u,i}$, imputation model
$g_{\phi}(x_{u,i})=\tilde{r}_{u,i}$, propensity model
$h_{\psi}(x_{u,i})=\hat{p}_{u,i}$, and calibration experts
$\\{c_{\omega_{k}}^{\textup{imp}}(\cdot)\\}_{k\in[K]}$,
$\\{c_{\omega_{k}}^{\textup{prop}}(\cdot)\\}_{k\in[K]}$. Among them, the
propensity model $h_{\psi}(\cdot)$ can be trained by itself as the other three
components are not related to its optimization process. In this paper, we
train the propensity model through binary classification between $\mathcal{O}$
and $\mathcal{D}\setminus\mathcal{O}$, as done in Section 3.2. Then, on top of
the frozen propensity model, we train the calibration experts
$\\{c^{\textup{prop}}_{\omega_{k}}(\cdot)\\}_{k\in[K]}$ and the assignment
network $A^{\textup{prop}}$ with $\mathcal{L}_{\textup{propcal}}$ in Eq.15.

Input : $\mathcal{D}$, $\mathcal{O}$, $\mathcal{O}_{\textup{val}}$, pre-
trained propensity model $f_{\psi}$, pre-trained calibration experts
$\\{c^{\textup{prop}}_{\omega_{k}}(p)\\}_{k\in[K]}$, pre-trained assignment
network $A^{\textup{prop}}$.

1 while  _epochs remain or stopping criteria is not satisfied_ do

2 for  _mini-batch $\mathcal{O}^{s}$ of $\mathcal{O}$_ do

3 Update imputation model $g_{\phi}$ with $\mathcal{L}^{\textup{cal}}_{imp}$
(Eq.16) Sample a batch of user-item pairs from $\mathcal{D}$

4 Update prediction model $f_{\theta}$ with $\mathcal{E}^{\textup{cal}}_{DR}$
(Eq.17)

5 for  _mini-batch $\mathcal{O}^{s}_{\textup{val}}$ of
$\mathcal{O}_{\textup{val}}$_ do

6 Update calibration experts
$\\{c^{\textup{imp}}_{\omega_{k}}(p)\\}_{k\in[K]}$ and assignment network
$A^{\textup{imp}}$ with $\mathcal{L}_{\textup{impcal}}$ (Eq.14)

7

Algorithm 1 Tri-level Joint Learning Framework

The other components need to be trained together since their loss functions
incorporate predictions from others. In our tri-level joint learning
framework, each component is trained as follows:

(1) The imputation model $g_{\phi}(x_{u,i})=\tilde{r}_{u,i}$ is optimized by
using the imputation loss function with the calibrated propensity scores. For
illustrative purposes, we show the imputation loss adopted in (Wang et al.,
2019; Li et al., 2023f) with the calibrated propensity scores:

(16)
$\mathcal{L}^{\textup{cal}}_{\textup{imp}}(\phi)=\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}\frac{o_{u,i}(e(\hat{r},\tilde{r})-e(\hat{r},r))^{2}}{\bar{p}_{u,i}},$

where $\bar{p}_{u,i}=\sum_{k=1}^{K}\beta^{\textup{prop}}_{u,k}\cdot
c^{\textup{prop}}_{\omega_{k}}(\hat{p})$ is the calibrated propensity.

(2) The calibration experts $\\{c^{\textup{imp}}_{\omega_{k}}(p)\\}_{k\in[K]}$
and the assignment network $A^{\textup{imp}}$ are trained with
$\mathcal{L}_{\textup{impcal}}$ in Eq.14 on top of the frozen imputation
model.

(3) Lastly, the prediction model $f_{\theta}$ is trained with the calibrated
imputed errors and the calibrated propensity scores:

(17)
$\mathcal{E}^{\textup{cal}}_{\textup{DR}}(\theta)=\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}\Big{(}\bar{e}_{u,i}+\frac{o_{u,i}(e_{u,i}-\bar{e}_{u,i})}{\bar{p}_{u,i}}\Big{)},$

where $\bar{p}_{u,i}$ is the calibrated propensity score and
$\bar{e}_{u,i}=e(\hat{r},\bar{r})$ is the calibrated imputed error computed
with calibrated pseudo label
$\bar{r}_{u,i}=\sum_{k=1}^{K}\beta^{\textup{imp}}_{u,k}\cdot
c^{\textup{imp}}_{\omega_{k}}(\tilde{r})$. Algorithm 1 presents the entire
procedure of the tri-level joint learning framework. It is noted that the
proposed method is orthogonal to existing DR estimators and can be used in
tandem with them. We can adopt any other existing imputation loss (Guo et al.,
2021; Li et al., 2023b; Dai et al., 2022a) instead of Eq.16, by substituting
the miscalibrated propensity score $\hat{p}$ with our calibrated propensity
score $\bar{p}$.

##  5\. Experiments

###  5.1. Experiment Setup

####  5.1.1. Datasets

To evaluate the performance of unbiased recommendations, we need an unbiased
test set along with an MNAR training set. Following recent studies (Schnabel
et al., 2016; Wang et al., 2019; Li et al., 2023b; Dai et al., 2022a), we
adopt two real-world datasets, including Coat (Schnabel et al., 2016) and
Yahoo!R3333http://research.yahoo.com/Academic_Relations. Additionally, we
adopt KuaiRec (Gao et al., 2022), a recently published dataset that has a
separate unbiased test set where the users are asked to rate all test items.
For Coat and Yahoo!R3, the observed rating $r_{u,i}$ is set to 1 if the
explicit rating is greater than 3 and is set to 0 otherwise, as done in
compared studies (Saito, 2020; Dai et al., 2022a; Li et al., 2023b, f). For
KuaiRec, the rating $r_{u,i}$ is set to 1 if the watching ratio is greater
than 1 and is set to 0 otherwise. Data statistics are presented in Table 3
(Appendix B). We hold out 10% of the training set as the validation set.

####  5.1.2. Compared methods

We validate the superiority of Doubly Calibrated Estimator with various
debiasing methods as follows: (1) IPS estimators: IPS (Schnabel et al., 2016),
SNIPS (Swaminathan and Joachims, 2015), AT (Saito, 2020), (2) EIB estimator:
CVIB (Wang et al., 2020), (3) Multi-task approaches: Multi-DR (Zhang et al.,
2020), ESCM2-DR (Wang et al., 2022), (4) DR estimators: DR-JL (Wang et al.,
2019), MRDR-JL (Guo et al., 2021), BRD-DR (Dai et al., 2022b), MR (Li et al.,
2023a), DR-MSE (Dai et al., 2022a), StableDR (Li et al., 2023f), TDR-CL (Li et
al., 2023b). Please note that we exclude methods utilizing a fraction of the
unbiased test set in the training phase (Bonner and Vasile, 2018; Chen et al.,
2021; Wang et al., 2021; Ding et al., 2022; Li et al., 2023d) for a fair
comparison.

####  5.1.3. Implementation details

Compared methods. For all methods compared, we use their author codes and
strictly follow the parameters and configurations as documented in their
papers and codes. For methods without any public source code, we use PyTorch
(Paszke et al., 2019) for the implementation. Additionally, in cases where the
best configurations are not provided for each dataset, we perform a grid
search for hyperparameters on the validation set. Following previous studies
(Schnabel et al., 2016; Wang et al., 2019; Saito et al., 2020; Saito, 2020; Li
et al., 2023f, b), the base architectures of the prediction, imputation, and
propensity models are chosen from either matrix factorization (Koren et al.,
2009) or neural collaborative filtering (He et al., 2017), based on the
validation performance. It is noted that we do not use the unbiased test set
for the propensity estimation for a fair comparison.

Table 1. Recommendation performance. The numbers in bold indicate the best results, while those underlined represent the best competitor. * denotes the statistical significance of the paired t-test at a 0.05 level when compared to the best competitor. Method | Coat | Yahoo!R3 | KuaiRec
---|---|---|---
MSE | AUC | N@5 | N@10 | MSE | AUC | N@5 | N@10 | MSE | AUC | N@50 | N@100
Naive | 0.2547 | 0.6828 | 0.6087 | 0.6747 | 0.2603 | 0.6528 | 0.6303 | 0.7584 | 0.2738 | 0.7510 | 0.7350 | 0.7571
SNIPS | 0.2438 | 0.7061 | 0.6145 | 0.6875 | 0.2493 | 0.6815 | 0.6451 | 0.7701 | 0.2673 | 0.7608 | 0.7507 | 0.7844
IPS | 0.2423 | 0.6935 | 0.6130 | 0.6824 | 0.2496 | 0.6757 | 0.6348 | 0.7641 | 0.2699 | 0.7601 | 0.7442 | 0.7768
CVIB | 0.2311 | 0.7006 | 0.6221 | 0.6991 | 0.2438 | 0.6823 | 0.6487 | 0.7709 | 0.2702 | 0.7540 | 0.7381 | 0.7612
AT | 0.2332 | 0.7020 | 0.6302 | 0.6984 | 0.2480 | 0.6816 | 0.6419 | 0.7667 | 0.2631 | 0.7581 | 0.7471 | 0.7803
DR-JL | 0.2302 | 0.7122 | 0.6382 | 0.7042 | 0.2459 | 0.6837 | 0.6515 | 0.7731 | 0.2580 | 0.7683 | 0.7689 | 0.7932
|  $\pm$ 0.0013 |  $\pm$ 0.0053 |  $\pm$ 0.0088 |  $\pm$ 0.0054 |  $\pm$ 0.0016 |  $\pm$ 0.0025 |  $\pm$ 0.0036 |  $\pm$ 0.0020 |  $\pm$ 0.0018 |  $\pm$ 0.0037 |  $\pm$ 0.0034 |  $\pm$ 0.0035
Multi-DR | 0.2288 | 0.7201 | 0.6410 | 0.7081 | 0.2407 | 0.6863 | 0.6587 | 0.7712 | 0.2557 | 0.7714 | 0.7771 | 0.7985
MRDR-JL | 0.2193 | 0.7192 | 0.6360 | 0.7016 | 0.2394 | 0.6842 | 0.6602 | 0.7727 | 0.2542 | 0.7705 | 0.7825 | 0.8023
ESCM2-DR | 0.2201 | 0.7256 | 0.6353 | 0.7022 | 0.2401 | 0.6932 | 0.6683 | 0.7693 | 0.2520 | 0.7811 | 0.7921 | 0.8049
BRD-DR | 0.2196 | 0.7286 |  0.6441 | 0.7094 | 0.2382 | 0.6850 | 0.6592 | 0.7708 | 0.2505 | 0.7831 |  0.7945 | 0.8051
MR | 0.2181 | 0.7243 | 0.6436 | 0.7118 |  0.2312 | 0.6923 | 0.6691 | 0.7747 | 0.2496 | 0.7851 | 0.7912 | 0.8017
DR-MSE |  0.2152 | 0.7214 | 0.6417 | 0.7089 | 0.2377 | 0.6872 | 0.6660 | 0.7724 | 0.2510 | 0.7792 | 0.7823 |  0.8067
|  $\pm$ 0.0018 |  $\pm$ 0.0042 |  $\pm$ 0.0066 |  $\pm$ 0.0064 |  $\pm$ 0.0012 |  $\pm$ 0.0018 |  $\pm$ 0.0022 |  $\pm$ 0.0011 |  $\pm$ 0.0011 |  $\pm$ 0.0018 |  $\pm$ 0.0027 |  $\pm$ 0.0032
StableDR | 0.2238 | 0.7167 | 0.6335 | 0.7073 | 0.2425 | 0.6891 | 0.6613 | 0.7699 | 0.2543 | 0.7714 | 0.7863 | 0.7976
|  $\pm$ 0.0028 |  $\pm$ 0.0061 |  $\pm$ 0.0087 |  $\pm$ 0.0086 |  $\pm$ 0.0033 |  $\pm$ 0.0080 |  $\pm$ 0.0048 |  $\pm$ 0.0029 |  $\pm$ 0.0020 |  $\pm$ 0.0061 |  $\pm$ 0.0052 |  $\pm$ 0.0048
TDR-CL | 0.2173 |  0.7302 | 0.6425 |  0.7121 | 0.2363 |  0.6951 |  0.6708 |  0.7763 |  0.2477 |  0.7885 | 0.7907 | 0.8037
|  $\pm$ 0.0005 |  $\pm$ 0.0045 |  $\pm$ 0.0098 |  $\pm$ 0.0050 |  $\pm$ 0.0010 |  $\pm$ 0.0015 |  $\pm$ 0.0010 |  $\pm$ 0.0013 |  $\pm$ 0.0006 |  $\pm$ 0.0008 |  $\pm$ 0.0011 |  $\pm$ 0.0012
DCE-DR | 0.2109 | 0.7384 | 0.6581 | 0.7232 | 0.2281 | 0.7041 | 0.6826 | 0.7923 | 0.2430 | 0.8058 | 0.8079 | 0.8122
|  $\pm$ 0.0004 |  $\pm$ 0.0019 |  $\pm$ 0.0030 |  $\pm$ 0.0043 |  $\pm$ 0.0027 |  $\pm$ 0.0013 |  $\pm$ 0.0018 |  $\pm$ 0.0014 |  $\pm$ 0.0005 |  $\pm$ 0.0006 |  $\pm$ 0.0022 |  $\pm$ 0.0016
DCE-TDR | 0.2054* | 0.7412* | 0.6623* | 0.7258* | 0.2247* | 0.7087* | 0.6901* | 0.8014* | 0.2394* | 0.8097* | 0.8162* | 0.8249*
|  $\pm$ 0.0006 |  $\pm$ 0.0032 |  $\pm$ 0.0045 |  $\pm$ 0.0061 |  $\pm$ 0.0014 |  $\pm$ 0.0015 |  $\pm$ 0.0016 |  $\pm$ 0.0012 |  $\pm$ 0.0010 |  $\pm$ 0.0012 |  $\pm$ 0.0021 |  $\pm$ 0.0026

Proposed method. For our Doubly Calibrated Estimator (DCE), we devise two
variants: DCE-DR utilizing $\mathcal{L}^{\textup{cal}}_{\textup{imp}}$ in
Eq.16 and DCE-TDR utilizing the imputation loss of TDR (Li et al., 2023b)
instead of Eq.16 with substituting the miscalibrated propensity score
$\hat{p}$ with our calibrated propensity score $\bar{p}$. For the assignment
network, we use a single-layer MLP with the softmax output layer. For the
Gumbel-Softmax, the initial temperature is set to 1 and the terminal
temperature is set to $10^{-3}$. We set the batch size as 128 for Coat, 4096
for Yahoo!R3 and KuaiRec. The other hyperparameters are tuned by using grid
searches on the validation set. We use Adam optimizer (Kingma and Ba, 2015)
where the learning rate and weight decay are chosen from the range
$[10^{-4},10^{-1}]$. The number of experts $K$ is chosen from $\\{5,10,20\\}$
and the embedding size is chosen from $\\{16,32,64\\}$. We provide the GitHub
repository of the proposed
method.444https://github.com/WonbinKweon/DCE_WWW2024

###  5.2. Performance Comparison

Table 1 shows the recommendation performance of the proposed and the compared
debiasing methods. We use MSE (Chai and Draxler, 2014), AUC (Huang and Ling,
2005), NDCG@$K$ (N@$K$) (Järvelin and Kekäläinen, 2002) for the evaluation
metrics, following the conventions in the recent studies (Li et al., 2023b, f;
Dai et al., 2022a). Due to the lack of space, we report standard deviation of
five independent runs only for DR-JL and its recent variants. We first observe
that both the IPS estimators and the EIB estimator outperform the Naive
estimator, indicating the importance of addressing the MNAR problem in
recommendation datasets. Also, DR-JL and its variants perform better than the
IPS and the EIB estimator based on their double robustness. The state-of-the-
art DR estimators manipulate the imputed errors (Li et al., 2023b; Dai et al.,
2022a) or the propensity scores (Li et al., 2023f) for the better bias-
variance trade-off and show enhanced unbiased recommendation performance.
However, their effectiveness has been limited as they still utilize the
miscalibrated imputed errors and propensity scores. The proposed methods (DCE-
DR and DCE-TDR) employ calibration experts for both the imputation and the
propensity models, effectively generating calibrated estimates which are
critical requirements for ensuring the unbiasedness of DR estimators. As a
result, the proposed methods successfully increase the recommendation
performance from the existing DR estimators by up to 8.37% (DCE-DR vs DR-JL in
MSE on Coat).

###  5.3. Ablation Study

Figure 2. Ablation study on Yahoo!R3.

Figure 2 shows the ablation study of the proposed doubly calibrated estimator.
We introduce four ablated methods crafted to showcase the superiority of our
design choices: (1) ’Global’ indicates DCE-DR with global calibration function
$c_{\omega}$ for all users. (2) ’User-specific’ represents DCE-DR with user-
specific calibration functions $\\{c_{\omega_{u}}\\}_{u\in\mathcal{U}}$. (3)
’Only_ImpCal’ denotes DCE-DR with calibration experts only for the imputation
model. (4) ’Only_PropCal’ denotes DCE-DR with calibration experts only for the
propensity model. (5) ’DR-JL’ (Wang et al., 2019) does not adopt any
calibration.

First, we observe that our calibration experts outperform the global and the
user-specific calibration (i.e., DCE-DR vs Global, User-specific). Calibration
experts can consider distinct logit distributions and learn specialized
information about their group. Also, they alleviate the cold-start problem of
user-specific calibration by gathering enough training signals from the users
in the group. As a result, the calibration experts produce more calibrated
imputed errors and propensity scores (Figure 2 left), and accordingly, DCE-DR
with calibration experts outperforms DCE-DR with global or user-specific
calibration in terms of the unbiased recommendation performance (Figure 2
right).

Second, we observe that employing calibration experts for both the imputation
and the propensity model is important for the effectiveness of DCE-DR (i.e.,
DCE-DR vs Only_ImpCal, Only_PropCal). As theoretically demonstrated in Section
3.1, the bias and the variance of DR estimators are positively correlated with
the miscalibration of both the imputation and the propensity models. Moreover,
the calibrated propensity scores enhance the calibration (Eq.14) and the
optimization (Eq.16) of the imputation model. Therefore, DCE-DR exhibits
superior performance when utilizing both calibrated imputed errors and
calibrated propensity scores, as compared to employing either one
individually.

###  5.4. Time and Space Analysis

Table 2. Time and space analysis. Time denotes the wall time (in s) used for
the training and #params. denotes the number of learnable parameters. Please
note that the base model and the embedding size are chosen with the grid
search on the validation set.

Method | Coat | Yahoo!R3 | KuaiRec
---|---|---|---
Time | #params. | Time | #params. | Time | #params.
MF | 8.16 | 9,440 | 98.87 | 524,864 | 211.47 | 1,136,576
IPS | 12.14 | 18,913 | 217.32 | 1,049,793 | 481.25 | 2,273,281
DR-JL | 19.04 | 28,353 | 460.04 | 3,149,313 | 538.39 | 3,409,857
DR-MSE | 22.14 | 37,772 | 308.01 | 2,099,488 | 697.18 | 4,546,432
TDR-CL | 24.58 | 28,354 | 319.86 | 1,574,660 | 702.74 | 3,409,864
DCE-DR | 24.94 | 28,523 | 320.52 | 1,575,317 | 710.91 | 3,412,457
DCE-TDR | 26.21 | 28,524 | 346.43 | 1,575,532 | 752.13 | 3,412,464

Figure 3. Hyper-parameter study on number of experts ($K$).

Table 2 shows the training time and the number of learnable parameters of
compared methods. In this experiment, We use GTX Titan Xp GPU and Intel
Xeon(R) E5-2667 v4 CPU. First, we observe that the proposed doubly calibrated
estimators do not increase the training time significantly (less than 10%
compared to TDR-CL). Each calibration expert is equivalent to the logistic
regression and can be efficiently solved. Also, the calibrated experts are
updated within 10-20 iterations on the small validation set. Second, doubly
calibrated estimators introduce only marginal additional learning parameters
in comparison to existing methods. Each calibration expert has two parameters
and thus the total number of additional parameters is $(4+2d)K$ (4$K$ for
calibration experts and $2dK$ for the assignment networks). Here, $d\leq 64$
and $K\leq 20$ (the hyper-parameter study on $K$ is presented in Figure 3).

##  6\. Related Work

Over the past recent years, debiasing methods have been extremely investigated
for accurately approximating the ideal loss over the entire population with
only biased MNAR data. EIB estimator (Steck, 2010; Wang et al., 2020) and IPS
estimator (Schnabel et al., 2016; Saito et al., 2020; Raghu et al., 2018) are
two major approaches for designing unbiased estimators. Lately, DR estimator
(Wang et al., 2019) and its enhanced variants (Guo et al., 2021; Dai et al.,
2022b; Li et al., 2023f, b; Zhang et al., 2020; Li et al., 2023a; Wang et al.,
2022) are proposed to merge EIB and IPS estimators for double robustness, and
have shown state-of-the-art performance for various user behaviors including
explicit ratings (Schnabel et al., 2016; Wang et al., 2019; Guo et al., 2021),
implicit feedback (Saito et al., 2020; Saito, 2020; Zhu et al., 2020), uplift
modeling (Sato et al., 2019; Saito et al., 2019), and post-click conversion
rate (Dai et al., 2022a; Li et al., 2023f, b; Zhang et al., 2020). The recent
trend is focusing on the development of estimators that are robust to the
inaccurate estimation of imputed error and propensity score, rather than
achieving precise estimation of them. Recent studies enhance the robustness of
DR estimators by manipulating the imputed errors (Li et al., 2023b; Dai et
al., 2022a; Guo et al., 2021) and the propensity scores (Li et al., 2023f;
Raghu et al., 2018), or by adopting multiple imputation and propensity models
(Li et al., 2023a). Nonetheless, their effectiveness may be limited since they
still depend on miscalibrated imputed errors and propensity scores.

On the other hand, some debiasing methods (Bonner and Vasile, 2018; Chen et
al., 2021; Wang et al., 2021; Ding et al., 2022; Li et al., 2023d; Liu et al.,
2020; Li et al., 2023e, c; Liu et al., 2022) utilize a small set of unbiased
data during training to enhance debiasing performance. For instance,
approaches like (Schnabel et al., 2016; Saito et al., 2020; Li et al., 2023f)
employ logistic regression on held-out unbiased data for accurate propensity
estimation. Other studies leverage unbiased data to address unobserved
confounding in selection bias (Li et al., 2023d; Wang et al., 2021) and for
causal inference (Bonner and Vasile, 2018; Chen et al., 2021; Luo et al.,
2023; Wu et al., 2022). These methods exhibit a slight improvement compared to
scenarios without unbiased uniform data, resulting in a moderately accurate
prediction model. Nonetheless, obtaining unbiased data in real-world
applications can be challenging, as employing a random exposure policy may
degrade user satisfaction. Indeed, they utilize a small fraction of an
unbiased test set for utilizing held-out uniform data. In contrast, our doubly
calibrated estimator operates without any unbiased data and effectively
estimates precise imputed errors and propensity scores.

##  7\. Conclusion

We demonstrate that current DR estimators rely on miscalibrated imputed errors
and propensity scores and provide theoretical analysis on limitations when
faced with miscalibration. Building upon these findings, we propose a Doubly
Calibrated Estimator, which involves calibrating both the imputation and
propensity models. To achieve this, we introduce calibration experts that
consider different logit distributions across users. Furthermore, we design a
tri-level joint learning framework, enabling the simultaneous optimization of
calibration experts alongside prediction and imputation models. Through
extensive experiments conducted on real-world datasets, we demonstrate the
enhanced performance of the doubly calibrated estimator in debiased
recommendation tasks. We anticipate that our doubly calibrated estimator can
be seamlessly integrated with various DR estimators in future work.

## Acknowledgements

This work was supported by the IITP grant funded by the MSIT (South Korea,
No.2018-0-00584, No.2019-0-01906), the NRF grant funded by the MSIT (South
Korea, No.2020R1A2B5B03097210, No. RS-2023-00217286), and NRF grant funded by
the MOE (South Korea, 2022R1A6A1A03052954).

## References

  * (1)
  * Bai et al. (2021) Yu Bai, Song Mei, Huan Wang, and Caiming Xiong. 2021\.  Don’t just blame over-parametrization for over-confidence: Theoretical analysis of calibration in binary classification. In _ICML_.
  * Bonner and Vasile (2018) Stephen Bonner and Flavian Vasile. 2018.  Causal embeddings for recommendation. In _RecSys_.
  * Chai and Draxler (2014) Tianfeng Chai and Roland R Draxler. 2014.  Root mean square error (RMSE) or mean absolute error (MAE)?–Arguments against avoiding RMSE in the literature.  _Geoscientific model development_ 7, 3 (2014), 1247–1250.
  * Chen et al. (2021) Jiawei Chen, Hande Dong, Yang Qiu, Xiangnan He, Xin Xin, Liang Chen, Guli Lin, and Keping Yang. 2021\.  AutoDebias: Learning to debias for recommendation. In _SIGIR_.
  * Chen et al. (2023) Jiawei Chen, Hande Dong, Xiang Wang, Fuli Feng, Meng Wang, and Xiangnan He. 2023\.  Bias and debias in recommender system: A survey and future directions.  _ACM Transactions on Information Systems_ 41, 3 (2023), 1–39.
  * Dai et al. (2022a) Quanyu Dai, Haoxuan Li, Peng Wu, Zhenhua Dong, Xiao-Hua Zhou, Rui Zhang, Rui Zhang, and Jie Sun. 2022a.  A generalized doubly robust learning framework for debiasing post-click conversion rate prediction. In _SIGKDD_.
  * Dai et al. (2022b) Quanyu Dai, Haoxuan Li, Peng Wu, Zhenhua Dong, Xiao-Hua Zhou, Rui Zhang, Rui Zhang, and Jie Sun. 2022b.  A generalized doubly robust learning framework for debiasing post-click conversion rate prediction. In _SIGKDD_.
  * Desai and Durrett (2020) Shrey Desai and Greg Durrett. 2020.  Calibration of Pre-trained Transformers. In _EMNLP_.
  * Ding et al. (2022) Sihao Ding, Fuli Feng, Xiangnan He, Jinqiu Jin, Wenjie Wang, Yong Liao, and Yongdong Zhang. 2022.  Interpolative distillation for unifying biased and debiased recommendation. In _SIGIR_.
  * Ding et al. (2021) Zhipeng Ding, Xu Han, Peirong Liu, and Marc Niethammer. 2021\.  Local temperature scaling for probability calibration. In _ICCV_.
  * Gao et al. (2022) Chongming Gao, Shijun Li, Wenqiang Lei, Jiawei Chen, Biao Li, Peng Jiang, Xiangnan He, Jiaxin Mao, and Tat-Seng Chua. 2022.  KuaiRec: A fully-observed dataset and insights for evaluating recommender systems. In _CIKM_.
  * Guo et al. (2017) Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q Weinberger. 2017\.  On calibration of modern neural networks. In _ICML_.
  * Guo et al. (2021) Siyuan Guo, Lixin Zou, Yiding Liu, Wenwen Ye, Suqi Cheng, Shuaiqiang Wang, Hechang Chen, Dawei Yin, and Yi Chang. 2021.  Enhanced doubly robust learning for debiasing post-click conversion rate estimation. In _SIGIR_.
  * He et al. (2017) Xiangnan He, Lizi Liao, Hanwang Zhang, Liqiang Nie, Xia Hu, and Tat-Seng Chua. 2017\.  Neural collaborative filtering. In _WWW_.
  * Huang and Ling (2005) Jin Huang and Charles X Ling. 2005.  Using AUC and accuracy in evaluating learning algorithms.  _IEEE Transactions on knowledge and Data Engineering_ 17, 3 (2005), 299–310.
  * Jacobs et al. (1991) Robert A Jacobs, Michael I Jordan, Steven J Nowlan, and Geoffrey E Hinton. 1991.  Adaptive mixtures of local experts.  _Neural computation_ 3, 1 (1991), 79–87.
  * Jang et al. (2016) Eric Jang, Shixiang Gu, and Ben Poole. 2016.  Categorical Reparameterization with Gumbel-Softmax. In _ICLR_.
  * Järvelin and Kekäläinen (2002) Kalervo Järvelin and Jaana Kekäläinen. 2002.  Cumulated gain-based evaluation of IR techniques.  _ACM Transactions on Information Systems (TOIS)_ 20, 4 (2002), 422–446.
  * Kang et al. (2020) SeongKu Kang, Junyoung Hwang, Wonbin Kweon, and Hwanjo Yu. 2020\.  DE-RRD: A Knowledge Distillation Framework for Recommender System. In _CIKM_.
  * Kingma and Ba (2015) Diederik P Kingma and Jimmy Ba. 2015.  Adam: A method for stochastic optimization. In _ICLR_.
  * Koren et al. (2009) Yehuda Koren, Robert Bell, and Chris Volinsky. 2009\.  Matrix factorization techniques for recommender systems.  _Computer_ 42, 8 (2009), 30–37.
  * Kull et al. (2017) Meelis Kull, Telmo Silva Filho, and Peter Flach. 2017\.  Beta calibration: a well-founded and easily implemented improvement on logistic calibration for binary classifiers. In _AISTATS_.
  * Kumar et al. (2019) Ananya Kumar, Percy S Liang, and Tengyu Ma. 2019.  Verified uncertainty calibration. In _NeurIPS_.
  * Kweon et al. (2024) Wonbin Kweon, SeongKu Kang, Sanghwan Jang, and Hwanjo Yu. 2024\.  Top-Personalized-K Recommendation. In _WWW_.
  * Kweon et al. (2022) Wonbin Kweon, SeongKu Kang, and Hwanjo Yu. 2022.  Obtaining Calibrated Probabilities with Personalized Ranking Models. In _AAAI_.
  * Li et al. (2023a) Haoxuan Li, Quanyu Dai, Yuru Li, Yan Lyu, Zhenhua Dong, Xiao-Hua Zhou, and Peng Wu. 2023a.  Multiple robust learning for recommendation. In _AAAI_.
  * Li et al. (2023b) Haoxuan Li, Yan Lyu, Chunyuan Zheng, and Peng Wu. 2023b.  TDR-CL: Targeted Doubly Robust Collaborative Learning for Debiased Recommendations. In _ICLR_.
  * Li et al. (2023c) Haoxuan Li, Kunhan Wu, Chunyuan Zheng, Yanghao Xiao, Hao Wang, Zhi Geng, Fuli Feng, Xiangnan He, and Peng Wu. 2023c.  Removing hidden confounding in recommendation: a unified multi-task learning approach. In _NeurIPS_.
  * Li et al. (2023d) Haoxuan Li, Yanghao Xiao, Chunyuan Zheng, and Peng Wu. 2023d.  Balancing unobserved confounding with a few unbiased ratings in debiased recommendations. In _WWW_.
  * Li et al. (2023e) Haoxuan Li, Yanghao Xiao, Chunyuan Zheng, Peng Wu, and Peng Cui. 2023e.  Propensity matters: Measuring and enhancing balancing for recommendation. In _ICML_. 20182–20194.
  * Li et al. (2023f) Haoxuan Li, Chunyuan Zheng, and Peng Wu. 2023f.  StableDR: Stabilized Doubly Robust Learning for Recommendation on Data Missing Not at Random. In _ICLR_.
  * Liu et al. (2020) Dugang Liu, Pengxiang Cheng, Zhenhua Dong, Xiuqiang He, Weike Pan, and Zhong Ming. 2020\.  A general knowledge distillation framework for counterfactual recommendation via uniform data. In _SIGIR_.
  * Liu et al. (2022) Haochen Liu, Da Tang, Ji Yang, Xiangyu Zhao, Hui Liu, Jiliang Tang, and Youlong Cheng. 2022.  Rating distribution calibration for selection bias mitigation in recommendations. In _WWW_. 2048–2057.
  * Luo et al. (2023) Huishi Luo, Fuzhen Zhuang, Ruobing Xie, Hengshu Zhu, Deqing Wang, Zhulin An, and Yongjun Xu. 2023.  A survey on causal inference for recommendation.  _The Innovation_ (2023).
  * Minderer et al. (2021) Matthias Minderer, Josip Djolonga, Rob Romijnders, Frances Hubis, Xiaohua Zhai, Neil Houlsby, Dustin Tran, and Mario Lucic. 2021\.  Revisiting the calibration of modern neural networks. In _NeurIPS_.
  * Naeini et al. (2015) Mahdi Pakdaman Naeini, Gregory Cooper, and Milos Hauskrecht. 2015\.  Obtaining well calibrated probabilities using bayesian binning. In _AAAI_.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. 2019\.  PyTorch: An imperative style, high-performance deep learning library. In _NeurIPS_.
  * Platt et al. (1999) John Platt, Alexander Smola, Peter Bartlett, Bernhard Scholkopf, and Dale Schuurmans. 1999\.  Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods.  _Advances in large margin classifiers_ 10, 3 (1999), 61–74.
  * Raghu et al. (2018) Aniruddh Raghu, Omer Gottesman, Yao Liu, Matthieu Komorowski, Aldo Faisal, Finale Doshi-Velez, and Emma Brunskill. 2018.  Behaviour policy estimation in off-policy policy evaluation: Calibration matters.  _arXiv preprint arXiv:1807.01066_ (2018).
  * Saito (2020) Yuta Saito. 2020\.  Asymmetric tri-training for debiasing missing-not-at-random explicit feedback. In _SIGIR_.
  * Saito et al. (2019) Yuta Saito, Hayato Sakata, and Kazuhide Nakata. 2019\.  Doubly robust prediction and evaluation methods improve uplift modeling for observational data. In _SDM_.
  * Saito et al. (2020) Yuta Saito, Suguru Yaginuma, Yuta Nishino, Hayato Sakata, and Kazuhide Nakata. 2020.  Unbiased recommender learning from missing-not-at-random implicit feedback. In _WSDM_.
  * Sato et al. (2019) Masahiro Sato, Janmajay Singh, Sho Takemori, Takashi Sonoda, Qian Zhang, and Tomoko Ohkuma. 2019\.  Uplift-based evaluation and optimization of recommenders. In _RecSys_.
  * Schnabel et al. (2016) Tobias Schnabel, Adith Swaminathan, Ashudeep Singh, Navin Chandak, and Thorsten Joachims. 2016\.  Recommendations as treatments: Debiasing learning and evaluation. In _ICML_.
  * Steck (2010) Harald Steck. 2010\.  Training and testing of recommender systems on data missing not at random. In _SIGKDD_.
  * Su et al. (2019) Yi Su, Lequn Wang, Michele Santacatterina, and Thorsten Joachims. 2019.  Cab: Continuous adaptive blending for policy evaluation and learning. In _ICML_.
  * Swaminathan and Joachims (2015) Adith Swaminathan and Thorsten Joachims. 2015.  The self-normalized estimator for counterfactual learning. In _NeurIPS_.
  * Wang et al. (2022) Hao Wang, Tai-Wei Chang, Tianqiao Liu, Jianmin Huang, Zhichao Chen, Chao Yu, Ruopeng Li, and Wei Chu. 2022\.  Escm2: Entire space counterfactual multi-task model for post-click conversion rate estimation. In _SIGIR_.
  * Wang et al. (2019) Xiaojie Wang, Rui Zhang, Yu Sun, and Jianzhong Qi. 2019\.  Doubly robust joint learning for recommendation on data missing not at random. In _ICML_.
  * Wang et al. (2021) Xiaojie Wang, Rui Zhang, Yu Sun, and Jianzhong Qi. 2021\.  Combating selection biases in recommender systems with a few unbiased ratings. In _WSDM_.
  * Wang et al. (2020) Zifeng Wang, Xi Chen, Rui Wen, Shao-Lun Huang, Ercan Kuruoglu, and Yefeng Zheng. 2020\.  Information theoretic counterfactual learning from missing-not-at-random feedback. In _NeurIPS_.
  * Wu et al. (2022) Peng Wu, Haoxuan Li, Yuhao Deng, Wenjie Hu, Quanyu Dai, Zhenhua Dong, Jie Sun, Rui Zhang, and Xiao-Hua Zhou. 2022.  On the Opportunity of Causal Learning in Recommendation Systems: Foundation, Estimation, Prediction and Challenges. In _IJCAI_.
  * Zhang et al. (2020) Wenhao Zhang, Wentian Bao, Xiao-Yang Liu, Keping Yang, Quan Lin, Hong Wen, and Ramin Ramezani. 2020.  Large-scale causal approaches to debiasing post-click conversion rate estimation with multi-task learning. In _WWW_.
  * Zhu et al. (2020) Ziwei Zhu, Yun He, Yin Zhang, and James Caverlee. 2020\.  Unbiased implicit recommendation and propensity estimation via combinational joint learning. In _RecSys_.

##  Appendix A Proofs

###### Theorem 2.

The bias of DR estimator exhibits an upper bound proportional to the
calibration error of the propensity model.

$\displaystyle\textup{Bias}[\mathcal{E}_{\textup{DR}}]\leq\rho_{\textup{max}}\cdot\textup{ECE}(h_{\psi}),$
$\displaystyle\textup{Bias}[\mathcal{E}_{\textup{DR}}]\leq\rho_{\textup{max}}\cdot\textup{MCE}(h_{\psi}),$

where
$\rho_{\textup{max}}=\textup{max}_{(u,i)\in\mathcal{D}}\lvert(e_{u,i}-\hat{e}_{u,i})/\hat{p}_{u,i}\rvert$.

###### Proof.

(1)

$\displaystyle\textup{Bias}[\mathcal{E}_{\textup{DR}}]$
$\displaystyle=\frac{1}{|\mathcal{D}|}\Bigl{|}\sum_{u,i\in\mathcal{D}}\Big{(}\frac{\hat{p}_{u,i}-p_{u,i}}{\hat{p}_{u,i}}\Big{)}(e_{u,i}-\hat{e}_{u,i})\Bigr{|}$
$\displaystyle\leq\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}|\hat{p}_{u,i}-p_{u,i}|\Bigl{|}\frac{e_{u,i}-\hat{e}_{u,i}}{\hat{p}_{u,i}}\Bigr{|}$
$\displaystyle\leq\rho_{\textup{max}}\cdot\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}|\hat{p}_{u,i}-p_{u,i}|,$

where
$\rho_{\textup{max}}=\text{max}_{(u,i)\in\mathcal{D}}\lvert(e_{u,i}-\hat{e}_{u,i})/\hat{p}_{u,i}\rvert$.
From Eq.8, we get the empirical ECE over $\mathcal{D}$.555Since the propensity
score is continuous in [0,1], we make the empirical assumption that there are
no pairs with the same propensity score. Please note that empirical ECE
underestimates the true ECE (Kumar et al., 2019), and thus the bound of
Theorem 2 is tighter than without the assumption.

$\displaystyle\text{ECE}(h_{\psi})$
$\displaystyle=\mathbb{E}_{\hat{p}}\big{[}\lvert\mathbb{E}[o|h_{\psi}(x)=\hat{p}]-\hat{p}\rvert\big{]}$
$\displaystyle=\sum_{\hat{p}}\lvert
P(o=1|h_{\psi}(x)=\hat{p})-\hat{p}\rvert\cdot P(\hat{p})$
$\displaystyle=\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}\lvert
P(o_{u,i}=1)-\hat{p}_{u,i}\rvert$
$\displaystyle=\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}\lvert
p_{u,i}-\hat{p}_{u,i}\rvert$

By combining the above two equations, we get

$\textup{Bias}[\mathcal{E}_{\textup{DR}}]\leq\rho_{\textup{max}}\cdot\textup{ECE}(h_{\psi}).$

(2)

Similar to (1), we get

$\displaystyle\textup{Bias}[\mathcal{E}_{\textup{DR}}]$
$\displaystyle\leq\rho_{\textup{max}}\cdot\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}|\hat{p}_{u,i}-p_{u,i}|$
$\displaystyle\leq\rho_{\textup{max}}\cdot\max_{u,i\in\mathcal{D}}|\hat{p}_{u,i}-p_{u,i}|$
$\displaystyle\text{MCE}(h_{\psi})$
$\displaystyle=\max_{\hat{p}\in[0,1]}\lvert\mathbb{E}[o|h_{\psi}(x)=\hat{p}]-\hat{p}\rvert$
$\displaystyle=\max_{\hat{p}\in[0,1]}\lvert
P(o=1|h_{\psi}(x)=\hat{p})-\hat{p}\rvert$
$\displaystyle=\max_{u,i\in\mathcal{D}}\lvert
P(o_{u,i}=1)-\hat{p}_{u,i}\rvert$
$\displaystyle=\max_{u,i\in\mathcal{D}}\lvert p_{u,i}-\hat{p}_{u,i}\rvert.$

Combining the above two equations, we get

$\textup{Bias}[\mathcal{E}_{\textup{DR}}]\leq\rho_{\textup{max}}\cdot\textup{MCE}(h_{\psi}).$

∎

###### Corollary 3.

The bias of DR estimator exhibits an upper bound proportional to the
calibration error of the imputation model.

$\displaystyle\textup{Bias}[\mathcal{E}_{\textup{DR}}]\leq\pi_{\textup{max}}\cdot\textup{ECE}(g_{\phi}),$
$\displaystyle\textup{Bias}[\mathcal{E}_{\textup{DR}}]\leq\pi_{\textup{max}}\cdot\textup{MCE}(g_{\phi}).$

where
$\pi_{\textup{max}}=\textup{max}_{(u,i)\in\mathcal{D}}\lvert(\hat{p}_{u,i}-p_{u,i})(e^{(1)}_{u,i}-e^{(0)}_{u,i})/\hat{p}_{u,i}\rvert$.

###### Proof.

(1)

$\displaystyle\textup{Bias}[\mathcal{E}_{\textup{DR}}]$
$\displaystyle=\frac{1}{|\mathcal{D}|}\Bigl{|}\sum_{u,i\in\mathcal{D}}\Big{(}\frac{\hat{p}_{u,i}-p_{u,i}}{\hat{p}_{u,i}}\Big{)}(e_{u,i}-\hat{e}_{u,i})\Bigr{|}$
$\displaystyle\leq\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}\Bigl{|}\frac{\hat{p}_{u,i}-p_{u,i}}{\hat{p}_{u,i}}\Bigr{|}|e_{u,i}-\hat{e}_{u,i}|$
$\displaystyle=\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}\Bigl{|}\frac{\hat{p}_{u,i}-p_{u,i}}{\hat{p}_{u,i}}\Bigr{|}|P(r_{u,i}=1)-\tilde{r}_{u,i}||e^{(1)}_{u,i}-e^{(0)}_{u,i}|$
$\displaystyle\because e-\hat{e}$
$\displaystyle=P(r=1)e^{(1)}+(1-P(r=1))e^{(0)}-\tilde{r}e^{(1)}-(1-\tilde{r})e^{(0)}$
$\displaystyle=(P(r=1)-\tilde{r})(e^{(1)}-e^{(0)})$

With
$\pi_{\textup{max}}=\text{max}_{(u,i)\in\mathcal{D}}\lvert(\hat{p}_{u,i}-p_{u,i})(e^{(1)}_{u,i}-e^{(0)}_{u,i})/\hat{p}_{u,i}\rvert$,
we get

$\textup{Bias}[\mathcal{E}_{\textup{DR}}]\leq\pi_{\textup{max}}\cdot\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}|P(r_{u,i}=1)-\tilde{r}_{u,i}|.$

Also, similar to the proof of Theorem 2, we get

$\text{ECE}(g_{\phi})=\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}|P(r_{u,i}=1)-\tilde{r}_{u,i}|.$

By combining the above two equations, we get

$\textup{Bias}[\mathcal{E}_{\textup{DR}}]\leq\pi_{\textup{max}}\cdot\textup{ECE}(g_{\phi}).$

(2)

Similar to (1), we get

$\displaystyle\textup{Bias}[\mathcal{E}_{\textup{DR}}]$
$\displaystyle\leq\pi_{\textup{max}}\cdot\frac{1}{|\mathcal{D}|}\sum_{u,i\in\mathcal{D}}|P(r_{u,i}=1)-\tilde{r}_{u,i}|$
$\displaystyle\leq\pi_{\textup{max}}\cdot\max_{u,i\in\mathcal{D}}|P(r_{u,i}=1)-\tilde{r}_{u,i}|$
$\displaystyle\text{MCE}(g_{\phi})$
$\displaystyle=\max_{\tilde{r}\in[0,1]}\lvert\mathbb{E}[r|g_{\phi}(x)=\tilde{r}]-\tilde{r}\rvert$
$\displaystyle=\max_{\tilde{r}\in[0,1]}\lvert
P(r=1|g_{\phi}(x)=\tilde{r})-\tilde{r}\rvert$
$\displaystyle=\max_{u,i\in\mathcal{D}}\lvert
P(r_{u,i}=1)-\tilde{r}_{u,i}\rvert.$

By combining the above two equations, we get

$\textup{Bias}[\mathcal{E}_{\textup{DR}}]\leq\pi_{\textup{max}}\cdot\textup{MCE}(g_{\phi}).$

∎

###### Theorem 4.

The variance of DR estimator exhibits an upper bound proportional to the
square of the calibration error of the imputation model.

$\displaystyle\textup{Var}[\mathcal{E}_{\textup{DR}}]\leq\omega_{\textup{max}}\cdot\big{(}\textup{ECE}(g_{\phi})\big{)}^{2},$
$\displaystyle\textup{Var}[\mathcal{E}_{\textup{DR}}]\leq\frac{\omega_{\textup{max}}}{|\mathcal{D}|}\cdot\big{(}\textup{MCE}(g_{\phi})\big{)}^{2},$

where $\omega_{\textup{max}}=\textup{max}_{(u,i)\in\mathcal{D}}\lvert
p_{u,i}(1-p_{u,i})(e^{(1)}_{u,i}-e^{(0)}_{u,i})^{2}/\hat{p}^{2}_{u,i}\rvert$.

###### Proof.

(1)

$\displaystyle\textup{Var}[\mathcal{E}_{\textup{DR}}]$
$\displaystyle=\frac{1}{|\mathcal{D}|^{2}}\sum_{u,i\in\mathcal{D}}\frac{p_{u,i}(1-p_{u,i})}{\hat{p}^{2}_{u,i}}(\hat{e}_{u,i}-e_{u,i})^{2}$
$\displaystyle=\frac{1}{|\mathcal{D}|^{2}}\sum_{u,i\in\mathcal{D}}\Bigl{|}\frac{p_{u,i}(1-p_{u,i})}{\hat{p}^{2}_{u,i}}\Bigr{|}|P(r_{u,i}=1)-\tilde{r}_{u,i}|^{2}|e^{(1)}_{u,i}-e^{(0)}_{u,i}|^{2}$
$\displaystyle\leq\omega_{\textup{max}}\cdot\frac{1}{|\mathcal{D}|^{2}}\sum_{u,i\in\mathcal{D}}|P(r_{u,i}=1)-\tilde{r}_{u,i}|^{2}$
$\displaystyle\leq\omega_{\textup{max}}\cdot\frac{1}{|\mathcal{D}|^{2}}\biggr{(}\sum_{u,i\in\mathcal{D}}|P(r_{u,i}=1)-\tilde{r}_{u,i}|\biggr{)}^{2}$
$\displaystyle=\omega_{\textup{max}}\cdot\big{(}\textup{ECE}(g_{\phi})\big{)}^{2},$

where $\omega_{\textup{max}}=\text{max}_{(u,i)\in\mathcal{D}}\lvert
p_{u,i}(1-p_{u,i})(e^{(1)}_{u,i}-e^{(0)}_{u,i})^{2}/\hat{p}^{2}_{u,i}\rvert$.
The second equal comes from the proof (1) of the Corollary 3.

(2)

Similar to (1), we get

$\displaystyle\textup{Var}[\mathcal{E}_{\textup{DR}}]$
$\displaystyle\leq\omega_{\textup{max}}\cdot\frac{1}{|\mathcal{D}|^{2}}\sum_{u,i\in\mathcal{D}}|P(r_{u,i}=1)-\tilde{r}_{u,i}|^{2}$
$\displaystyle\leq\omega_{\textup{max}}\cdot\frac{1}{|\mathcal{D}|}\max_{u,i\in\mathcal{D}}|P(r_{u,i}=1)-\tilde{r}_{u,i}|^{2}$
$\displaystyle=\omega_{\textup{max}}\cdot\frac{1}{|\mathcal{D}|}\cdot\Bigl{(}\max_{u,i\in\mathcal{D}}|P(r_{u,i}=1)-\tilde{r}_{u,i}|\Bigr{)}^{2}$
$\displaystyle=\frac{\omega_{\textup{max}}}{|\mathcal{D}|}\cdot\big{(}\textup{MCE}(g_{\phi})\big{)}^{2}$

∎

##  Appendix B Data Statistics

Table 3. Data statistics. Dataset | #Users | #Items | #Training data | #Test data
---|---|---|---|---
Coat | 290 | 300 | 6,960 | 4,640
Yahoo!R3 | 15,401 | 1,001 | 311,704 | 54,000
KuaiRec | 7,163 | 10,596 | 201,171 | 117,113


# Deep Neural Convolutive Matrix Factorization for Articulatory Representation
Decomposition

###### Abstract

Most of the research on data-driven speech representation learning has focused
on raw audios in an end-to-end manner, paying little attention to their
internal phonological or gestural structure. This work, investigating the
speech representations derived from articulatory kinematics signals, uses a
neural implementation of convolutive sparse matrix factorization to decompose
the articulatory data into interpretable gestures and gestural scores. By
applying sparse constraints, the gestural scores leverage the discrete
combinatorial properties of phonological gestures. Phoneme recognition
experiments were additionally performed to show that gestural scores indeed
code phonological information successfully. The proposed work thus makes a
bridge between articulatory phonology and deep neural networks to leverage
informative, intelligible, interpretable,and efficient speech representations.
The code is made publicly available at https://github.com/Berkeley-Speech-
Group/ema_gesture.

Index Terms: Articulatory Phonology, Gesture, Gestural Score

##  1 Introduction

Research on speech representation learning has been dominated by deep learning
recently in the areas such as speaker recognition [1, 2], automatic speech
recognition [3], voice conversion [4, 5] and text-to-speech [6, 7], etc. The
goal of speech representation learning is to optimize both the performance of
the model architectures and the interpretability of the learned
representations. As there is growing demand of real-life applications of
speech interfaces [8], the performance is emphasized to a larger extent,
enabling human-machine interactions highly accurate and robust. Consequently,
in most of these works the interpretability of representations has not been
explored to an equivalent extent, which is one of the most significant
bottlenecks that keeps the speech research from going farther. In general,
speech representations need to be better understood and developed.

People usually represent speech via audio because human perceive speech
through hearing and audio is cheap to record, collect and process. However,
speech processing is quite a lot different from audio processing. It might not
need any evidence to indicate that any information that can be perceived via
human can be perceived anywhere from source to destination. Perceiving the
speech signal from the source and leveraging how it is produced are the most
straightforward way to interpret it. The speech signal is the result of
respiratory, phonatory and articulatory processes that generate the
perceivable acoustic resonances to encode an intended linguistic message [9].
In that sense, perceiving the speech signal from articulatory data is a
preferred way to derive interpretable, natural and robust speech
representations.

Figure 1: Neural Convolutive Matrix Factorization to Interpret Gestures and
Gestural Scores.The panel in the bottom visualizes four activated gestures in
the example utterance for ”Five” . The gestures capture the moving patterns of
articulators. Each moving pattern goes from thinner line to thicker line
capturing about 200ms.

The framework of articulatory phonology [10] has offered a lawful approach to
modeling the relation between phonological representations as a set of
discrete compositional units, or gestures, and the variability in time that
derives from variation in the activation of the gestures in real-time: the
magnitude of their activation, and the temporal intervals of activation as
represented in gestural scores. However, the gestures and gestural scores of
particular utterances have never been estimated in a completely data-driven
manner. [11] utilized the convolutive sparse non-negative matrix factorization
(CSNMF) to decompose the non-negative articulatory data into the gestures and
gestural scores, both of which are pretty much interpretable. The downsides of
such method are that all the training utterances have to be concatenated into
a large matrix, resulting in both memory and training efficiency issues.
Additionally, such a model is not compatible with the modern deep learning
based speech models so that it is challenging to perform end-to-end training
on articulatory data.

To handle the aforementioned problem, [12] proposed an auto-encoder based
model to replace non-negatve matrix factorization for speech separation task.
Inspired by this work, we propose a convolutional auto-encoder as the neural
implementation of convolutive matrix factorization. Such auto-encoder based
matrix factorization method is compatible with modern deep neural network and
the batch-wise optimization improves the convergence rate to the huge extent.
Under such framework, the articulatory signal is decomposed into gestures and
gestural scores which are still interpretable. The gestural scores are the
learned articulatory speech representations and are constrained to be sparse.
In the last stage, the phoneme recognition experiments were performed to show
that the learned gestural scores are also intelligible and consistent in time
domain. All the experiments are performed using MNGU0 EMA (Electromagnetic
midsagittal articulography) [13] corpus. The intention is that the proposed
work could bridge the gap between explainable articulatory phonology and
modern deep neural networks to deliver interpretable, intelligible,
informative, and efficient speech representations.

##  2 Proposed Methods

###  2.1 Neural Convolutive Sparse Matrix Factorization

Denote EMA data as $X\in\mathbb{R}^{C\times t}$, where (C,t) is (number of
channels, segment length). By convolutive matrix factorization [11]:

$X\approx\Sigma_{i=0}^{T-1}W(i)\cdot\overrightarrow{H}^{i}$ (1)

$W\in\mathbb{R}^{T\times C\times D}$ is gestures and $H\in\mathbb{R}^{D\times
t}$ is gestural scores, where $D$ is number of gestures and $T$ is the kernel
size. $\overrightarrow{H}^{i}$ indicates that $i$ columns of $H$ are shifted
to the right.

It is observable that Eq. 1 is actually the 1-d convolution with kernel $W$
and input matrix $H$. By auto-encoder matrix factorization [12], H should be
the hidden representation derived from the encoder which takes the pseudo-
inverse of $W$ as parameters. However, calculating the pseudo-inverse of high
dimensional matrix is challenging. We experimentally justified that the
encoder can be any types of neural networks with any number of layers. The
proposed neural convolutive sparse matrix factorization is formularized as
follows:

$H=max(f(X),0)$ (2) $\hat{X}=W\odot H$ (3)

where $f(.)$ denotes any type of neural network. In the original non-negative
matrix factorization problem, all components ($X,W,H$) have to be non-
negative. However, in such neural implementation, only $H$ is required to be
non-negative so that the gestures are always additive. There is no constraint
for $W$ and $X$.

###  2.2 Loss Objectives

There are a couple of items in the loss function. The first one is the
reconstruction loss, which is L2 loss. The second one is sparseness. According
to [14], the sparseness of a vector in time dimension and channel dimension
are defined as in Eq. (4) and Eq. 5 respectively:

$S_{1}(H_{i})=\frac{\sqrt{t}-\frac{L_{1}(H_{i})}{L_{2}(H_{i})}}{\sqrt{t}-1}$
(4)
$S_{2}(H_{i})=\frac{\sqrt{D}-\frac{L_{1}(H^{T}_{i})}{L_{2}(H^{T}_{i})}}{\sqrt{D}-1}$
(5)

where $H_{i}$and $H^{T}_{i}$ denote the i-th row vector and column vector
respectively. $L_{1}$ and $L_{2}$ denote $L_{1}$ norm and $L_{2}$ norm
respectively. $n$ is the length of the vector. The sparseness of gestural
score matrix in time dimension and channel dimension are shown in Eq. (6)(7)
respectively. The intuition is that at each time step, there are not too many
gestures activated, and each gesture is not activate for a long time.

$S_{1}(H)=\frac{1}{D}\Sigma_{i=1}^{D}S(H_{i})$ (6)
$S_{2}(H)=\frac{1}{D}\Sigma_{i=1}^{t}S(H^{T}_{i})$ (7)

Following [15], we also introduce the entropy of the sparseness, denoted as:

$E(H)=\frac{1}{D}\Sigma_{i=1}^{D}(-\frac{S_{1}(H_{i})}{\Sigma_{i=1}^{D}S_{1}(H_{i})}\log(\frac{S_{1}(H_{i})}{\Sigma_{i=1}^{D}S_{1}(H_{i})}))$
(8)

It should be noticed that the sparseness cannot control the number of gestures
that are activated accurately. For instance, the H matrix with only one
gesture activated for a long time interval might have the same sparsity with
the matrix with multiple gestures activated for shorter time intervals.
Typically we expect that a proper number of gestures should be activated. Fig.
2 gives an intuition of entropy loss. All three H matrices have the same
sparsity. If the entropy is pretty low, only one gesture is activated, which
leaves many other gestures unused. If the entropy is pretty high, some of
activated gestures are redundant, which makes the gestural score less
explainable. We do not consider the entropy loss on the channel dimension
since it empirically does not make significant difference to the gestural
score.

(a) Low Entro

(b) Medium Entro

(c) High Entro

Figure 2: Three types of gestural score matrices with the same sparsity but
different entropy values.

We introduce balanced factors $\lambda_{1},\lambda_{2}$ and $\lambda_{3}$ to
limit both sparsity and entropy to a certain range. L2 loss is used as the
reconstruction loss, as shown below:

$L_{rec}=||X-\hat{X}||_{2}$ (9)

For EMA resynthesis task, the loss function is shown in Eq. 10, where
$\mathbb{E}_{X}$ means the loss is computed by taking the average in the mini-
batch.

$L_{res}=\mathbb{E}_{X}[L_{rec}-\lambda_{1}S_{1}(H)-\lambda_{2}S_{2}(H)+\lambda_{3}E(H)]$
(10)

For phoneme recognition experiments, CTC [16] loss $L_{CTC}$ is used. For
joint resynthesis-phoneme recognition task, the loss function is shown as in
Eq. 11, where $\lambda_{4}$ is a balanced factor.

$L_{joint}=L_{res}+\lambda_{4}L_{CTC}$ (11)

##  3 Experiments

###  3.1 Dataset

MNGU0 EMA (Electromagnetic midsagittal articulography) [13] dataset is used in
this work. There are in total 1263 utterances recorded from one single
speaker. During the recording, six transducer coils were placed in the
midsagittal plane at the upper lip, lower lip, lower incisors, tongue tip,
tongue blade and tongue dorsum to record the coordinates (x and y) of their
positions, and thus each EMA data frame takes 12 coordinates, as shown in Fig.
1. The sampling rate of EMA is 200 Hz. The Mel-Spectrogram is used as acoustic
feature with the framing configuration of 25ms/16ms and feature dimension of
80. The unaligned phonemes extracted from text transcriptions via the CMU
pronouncing dictionary111http://www.speech.cs.cmu.edu/cgi-bin/cmudict, are
used as labels for phoneme recognition task. The train/test split is 8:2,
which is the same for all experiments.

###  3.2 Tasks and Evaluation Methods

We perform two sets of experiments: (i) EMA Resynthesis. By resynthesizing the
EMA data, we extract, visualize and interpret the gestures and gestural
scores. The reconstruction loss $L_{rec}$ shown in Eq. 9 averaged over all
test samples is used to measure the informativeness of gestural scores [17].
The sparsity defined in Eq. 6 and 7 is used to measure the efficiency of
gestural scores. (ii) Phoneme Recognition (PR). PER (Phoneme Error Rate) is
used as metric for this task. PR on EMA is performed to measure the
intelligibility of EMA data. PER on melspectrogram is performed to measure the
intelligibility gap between articulartory and acoustics data. Lastly, the
joint training of EMA resynthesis and phoneme recognition on gestural scores
is performed to measure the intelligibility [18] of learned sparse speech
representations. Considering that EMA is not able to capture the difference
between voiced and voiceless phones, we also relabel the phoneme sequence by
assigning the same label to the phonemes with the same articulatory
representation in EMA222Specifically,these tuples are expected to have the
same articulatory labels: (p,b,m), (t,d,n), (ch,jh), (f,v), (sh,zh), (k,g,ng),
(s,z), (th,dh), and compute PER on new labels. We call the latter metric as
PER-V, which is reported for all PR experiments. The interpretability of
gestures and gestural scores is evaluated by subjective analysis given a set
of utterances. We also subjectively measure the consistency of the gestural
scores via visualizing a set of phonetic units across different utterances.

###  3.3 Model Architectures

The overall model backbone is shown in Fig. 1. The encoder takes EMA data $X$
in and outputs the gestural scores $H$. The decoder takes $H$ in and
resynthesizes EMA data $\hat{X}$. For phoneme recognition or joint
resynthesis-CTC experiments, the phoneme recognizer takes EMA, melspectrogram
or $H$ in and predicts the alignment. Beamsearch algorithm is used for
decoding with beam width of 50 in phoneme recognition task.

Module Name Block name Configurations Encoder Conv1d_1 (15,64)$\times 1$
Conv1d_2 (5,D)$\times 1$ Decoder Conv1d_3 (41,C)$\times 1$ Phoneme Recognizer
Conv1d_4 (5,64)$\times 3$ Bi-LSTM (256)$\times 3$ Linear (128)$\times 2$

Table 1: Model Details. For Conv1d Block, the configuration is (kernel size,
output channels)$\times$# of layers. For Bi-LSTM and Linear layers, the
configuration is (output dimension)$\times$# of layers.the number of gestures
$D$ is a hyperparameter and $C$ is 12. For all convolutional layers, the
stride is 1 and paddings are made so that the output length keeps the same
across layers. Batchnorm1D [19] is applied after each convolutional layer.

###  3.4 Implementation Details

For EMA resynthesis experiments, we randomly extract a segment with fixed
length of 300 frames as the input of model for each iteration. For phoneme
recognition experiments, the full utterance is taken as input. All experiments
were trained on Nvidia Tesla V100 GPU. It takes one GPU hour to run a single
EMA resynthesis experiment and 5GPU hours to run phoneme recognition as well
as resynthesis-CTC experiments. Optimizer is Adam [20] with the initial
learning rate of 1e-3, which is decayed every 5 epoches with a factor of 5.
Weight decay is 1e-4. Batchsize is 8. The weights of decoder (gestures) are
initialized by the centers from a kmeans algorithm: Slide the window of size
41 with a stride of 1 on EMA kinematics data, concatenate all 41 vectors into
a supervector and perform kmeans on all supervectors. For the loss function in
Eq. 10 and Eq. 11, we set $\lambda_{1}=\lambda_{2}=\lambda_{3}=10$ and
$\lambda_{3}=1$. For resynthesis and resynthesis-CTC experiments, we explore
different values of number of gestures: 20, 40, 60 and 80 as ablation studies.
The results of EMA resynthesis, joint resynthesis-CTC and independent PR on
EMA and melspectrogram are recorded in Table. 2 and Table. 3 respectively. To
interpret the gestures and gestural scores, a set of utterances are fed into
the encoder-decoder framework and we visualize the gestural scores as well as
activated gestures. We perform subjective analysis for each utterance and also
observe the consistency of the phonetic units across different utterances. A
tiny example which takes ”Five” as input is shown in Fig. 1.

Table 2: Resynthesis and Resynthesis-CTC

$\\#$gestures 20 40 60 80 Resynthesis Rec Loss $L_{rec}$ $\%$ 27.16 25.17
24.17 22.99 Sparsity $S_{1}(H)$ $\%$ 94.10 94.50 94.17 94.90 Sparsity
$S_{2}(H)$ $\%$ 92.09 91.78 93.77 90.66 Resynthesis-CTC Rec Loss $L_{rec}$
$\%$ 24.70 19.65 18.95 17.72 Sparsity $S_{1}(H)$ $\%$ 92.90 92.54 93.10 92.50
Sparsity $S_{2}(H)$ $\%$ 90.01 90.52 92.99 91.33 PER $\%$ 20.75 14.10 15.44
15.71 PER-V $\%$ 16.55 11.02 11.88 12.09

Table 3: PER on EMA and Melspec

Feature EMA Melspec PER $\%$ 13.27 7.54 PER-V $\%$ 10.24 6.18

###  3.5 Discussion

We discuss the results in terms of five aspects of the learned gestural
scores: Informativeness, Intelligibility, Efficiceny, Interpretability and
Consistency.

#### Informativeness

Lower reconstruction loss shows that the gestural scores are more informative.
By making the comparison between the input EMA and synthesized EMA, we
empirically observe that the reconstruction loss that is below 40$\%$ would
not loss too much information. As shown in Table. 2, the larger the number of
gestures, the more informative the gestures are. As it is not hard to overfit
the EMA data with an auto-encoder architecture, however, we apply hard
sparsity constraint which makes reconstruction much more challenging.

#### Intelligibility

Based on Table. 3, EMA gives higher PER and PER-V than melspectrogram because
EMA data is sparsely collected from articulators. PER-V of EMA is lower than
PER, which is consistent to the fact that EMA is not able to differentiate
voiced and voiceless phones. Based on Table. 2, when number of gestures is 40,
both PER and PER-V are comparable to the results obtained from EMA
representation, which shows that gestures scores are intelligible. Note that
when increasing the number of gestures, the PER is not always decreasing,
indicating that the intelligibility is not always positive correlated to the
informativeness. We believe that better intelligibility will be achieved when
using more fine-grained articulatory representations such as [21].

#### Efficiency

Based on Table. 2, when number of gestures is 40, the sparsity of gestural
scores in both time and channel dimension is more than $90\%$, however, the
intelligibility does not degrade too much in comparison to the original EMA
representation (13.27 versus 14.10 for PER and 10.24 versus 11.02 for PER-V).
As mentioned in Intelligibility, if more fine-grained articulatory
representations such as [21] are applied, it would be possible to still
achieve the gestural scores that as intelligible as melspectrogram but much
more efficient.

#### Interpretability and Consistency

We pick the resynthesis-CTC model with the best phoneme recognition
performance ($\\#$gestures is 40). Three pairs of short utterances(words) are
passed to generate the gestural scores. We visualize both gestural scores and
gestures that are activated during the pronunciation. In each gesture, the
moving pattern of articulators goes from thinner to thicker. The results are
presented in Fig. 3. The left six figures are the gestural scores and four
gestures that are activated. The solid red line denotes the rough estimation
of phoneme duration. The first pair is ”Five” and ”I” which share the same
vowel /ai/. Looking at the word ”I”, only gesture 9 is activated and gesture 9
exactly carries the moving patterns of ”I”: when pronouncing ”I”, we first
lower the lower lip, lower incisor and tongue articulators and then raise
them. When looking at gestural score of ”Five”, we observe consistent pattern
for /ai/, which is still activated at gesture 9, as indicated by the red
block. Other than /ai/, the phoneme /f/ is activated at gesture 10, where all
articulators except for upper lib are moving down, which reflects the real
articulatory patterns of /f/. For the phoneme /v/, gesture 36 is mainly
activated. In gesture 36, all articulators except for upper lip and tongue
dorsum first move down and then move up, which exactly reflects the
articulatory patterns of /v/. When looking at ”sharp” and ”jar”, it is
observed that /ar/ is the consistent pattern that is activated at gesture 22,
where both the lower lip and lower incisor move down and tongue articulators
move up, which exactly reflects the articulatory patterns of /ar/. For
/$\int$/ and /p/, they are mainly activated at gesture 10 where each
articulator is lowered. When looking at ”shame” and ”it”, we also observe
consistent patterns for /ei/ and /i/, which are mainly activated at gesture 9.
The phoneme /$\int$/ is also activated at gesture 9. Note that in gesture 9,
each articulator first moves down and then moves up. Such pattern is able to
model different phoneme units depending on how it is activated. If the latter
part (thicker) of the gesture 9 is activated, it indicates that both the lower
lip and lower incisor move down and the tongue articulators move up, which
corresponds to the articulatory patterns of /$\int$/. If the former part
(thinner) of gesture 9 is activated, it indicates that we just lower all
articulators, which corresponds to /ei/ and /i/. Note the articulatory
patterns for /m/ and /t/ are also consistent. When gesture 36 is activated, we
first lower all articulators except for upper lip and then raise them. In
conclusion, both gestures and gestural scores are interpretable and the
gestural scores also leverage consistent articulatory pattern across
utterances.

Figure 3: Gestural Scores Visualization.

##  4 Conclusion and Limitations

This work proposes a neural convolutive sparse matrix algorithm which
decomposes the EMA data into gestures and gestural scores. The learned
representations a.k.a gestural scores are informative, intelligent,
consistent, efficient and interpretable. This method bridges the gap between
articulatory phonology and deep learning techniques. Hopefully the proposed
work could become a paradigm that benefits the downstream explorations that
are helpful for patients with vocal cord disorders. One limitation is that EMA
data is sparsely sampled from articulators and the learned representation is
still less intelligible than the acoustic features. The future work will focus
on fine-grained articulatory representations such as [21] to deliver more
generalizable and intelligible representations.

## References

  * [1] J. Lian, A. V. Kumar, H. Dhamyal, B. Raj, and R. Singh, “Masked proxy loss for text-independent speaker verification,” _arXiv preprint arXiv:2011.04491_ , 2020.
  * [2] B. Desplanques, J. Thienpondt, and K. Demuynck, “ECAPA-TDNN: Emphasized channel attention, propagation and aggregation in TDNN based speaker verification,” in _Interspeech 2020_. ISCA, oct 2020. [Online]. Available: https://doi.org/10.21437%2Finterspeech.2020-2650
  * [3] A. Gulati, J. Qin, C.-C. Chiu, N. Parmar, Y. Zhang, J. Yu, W. Han, S. Wang, Z. Zhang, Y. Wu, and R. Pang, “Conformer: Convolution-augmented transformer for speech recognition,” 2020. [Online]. Available: https://arxiv.org/abs/2005.08100
  * [4] J. Lian, C. Zhang, and D. Yu, “Robust disentangled variational speech representation learning for zero-shot voice conversion,” in _IEEE ICASSP_. IEEE, 2022.
  * [5] J. Lian, C. Zhang, G. K. Anumanchipalli, and D. Yu, “Towards improved zero-shot voice conversion with conditional dsvae,” _arXiv preprint arXiv:2205.05227_ , 2022.
  * [6] J. Kim, J. Kong, and J. Son, “Conditional variational autoencoder with adversarial learning for end-to-end text-to-speech,” 2021. [Online]. Available: https://arxiv.org/abs/2106.06103
  * [7] J. Lian, C. Zhang, G. K. Anumanchipalli, and D. Yu, “Utts: Unsupervised tts with conditional disentangled sequential variational auto-encoder,” _arXiv preprint arXiv:2206.02512_ , 2022.
  * [8] C. Herff and T. Schultz, “Automatic speech recognition from neural signals: a focused review,” _Frontiers in neuroscience_ , vol. 10, p. 429, 2016.
  * [9] P. F. MacNeilage, _The origin of speech_. Oxford University Press, 2010, no. 10.
  * [10] C. P. Browman and L. Goldstein, “Articulatory phonology: An overview,” _Phonetica_ , vol. 49, no. 3-4, pp. 155–180, 1992.
  * [11] V. Ramanarayanan, L. Goldstein, and S. S. Narayanan, “Spatio-temporal articulatory movement primitives during speech production: Extraction, interpretation, and validation,” _The Journal of the Acoustical Society of America_ , vol. 134, no. 2, pp. 1378–1394, 2013.
  * [12] P. Smaragdis and S. Venkataramani, “A neural network alternative to non-negative audio models,” in _2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2017, pp. 86–90.
  * [13] K. Richmond, P. Hoole, and S. King, “Announcing the electromagnetic articulography (day 1) subset of the mngu0 articulatory corpus,” in _Twelfth Annual Conference of the International Speech Communication Association_ , 2011.
  * [14] P. O. Hoyer, “Non-negative matrix factorization with sparseness constraints.” _Journal of machine learning research_ , vol. 5, no. 9, 2004.
  * [15] A. Baevski, H. Zhou, A. Mohamed, and M. Auli, “wav2vec 2.0: A framework for self-supervised learning of speech representations,” _arXiv preprint arXiv:2006.11477_ , 2020.
  * [16] A. Graves, S. Fernández, F. Gomez, and J. Schmidhuber, “Connectionist temporal classification: labelling unsegmented sequence data with recurrent neural networks,” in _Proceedings of the 23rd international conference on Machine learning_ , 2006, pp. 369–376.
  * [17] A. M. Saxe, Y. Bansal, J. Dapello, M. Advani, A. Kolchinsky, B. D. Tracey, and D. D. Cox, “On the information bottleneck theory of deep learning,” _Journal of Statistical Mechanics: Theory and Experiment_ , vol. 2019, no. 12, p. 124020, 2019.
  * [18] K. Lakhotia, E. Kharitonov, W.-N. Hsu, Y. Adi, A. Polyak, B. Bolte, T.-A. Nguyen, J. Copet, A. Baevski, A. Mohamed _et al._ , “Generative spoken language modeling from raw audio,” _arXiv preprint arXiv:2102.01192_ , 2021\.
  * [19] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” in _International conference on machine learning_. PMLR, 2015, pp. 448–456.
  * [20] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [21] T. Sorensen, Z. I. Skordilis, A. Toutios, Y.-C. Kim, Y. Zhu, J. Kim, A. C. Lammert, V. Ramanarayanan, L. Goldstein, D. Byrd _et al._ , “Database of volumetric and real-time vocal tract mri for speech science.” in _Interspeech_ , 2017, pp. 645–649.