\section{Introduction}
\label{sec:intro}

Sound event localization and detection (SELD) is an audio processing task that aims to jointly detect temporally target classes of sound events and localize them in space when active. In that sense it differs from the classic sensor array task of sound source localization (SSL) which utilizes only spatial information to detect, localize, and track sources independently from their signal content \cite{evers2020locata}. It also differs from the popular sound event detection (SED) task which is focused on the temporal detection and classification part, omitting the spatial information of the scene. The spatiotemporal characterization of the scene produced by SELD makes it suitable for a range of applications such as robot audition and machine listening in general \cite{valin2003robust, He2018}, acoustic monitoring \cite{valenzise2007scream, do2018rish}, smart home environments \cite{do2018rish,brdiczka2008learning}, improved human-machine interaction \cite{cech2013active}, speech recognition \cite{he2018deep}, and sonic information visualization \cite{matsinos2008spatio}, among others.

Research interest in SELD grew quickly during the last couple of years, with deep learning methods handling the task jointly \cite{adavanne2018sound}, or fusing information from solving individual subtasks of SED and SSL \cite{Cao2019, nguyen2020sequence}. This interest culminated in the task becoming part of the DCASE Challenge
in 2019, with participants bringing novel approaches to the problem, summarized in \cite{politis2020overview}. The dataset used in the challenge \cite{Politis2020} included sound scenes from two different array formats with sound events spatialized in both azimuth and elevation using spatial room impulse responses (SRIRs) of real rooms. Additionally, spatial ambient noise captured in situ was added to the recordings. For the next iteration of the task in the DCASE Challenge 2020,
 a new dataset was generated based on SRIRs from additional rooms with more realistic and challenging conditions beyond the limitations of the first one \cite{Politis2020}. More specifically, the discrete grid of potential directions-of-arrival (DOAs) of the older dataset was replaced with continuous DOA trajectories and, apart from static events, moving sources using interpolated SRIRs were emulated at different speeds. Furthermore, the newer SRIRs were captured in rooms of more diverse acoustical properties and from a wider range of distances, resulting in longer reverberation times and more challenging direct-to-reverberant ratios (DRRs).

The second iteration of the SELD task in DCASE2020 brought additional innovations, with participants experimenting with homogeneous joint loss functions \cite{shimada2021accdoa, Phan2020}, self-attention layers \cite{Phan2020, wang2021four}, advanced spatial augmentation strategies \cite{shimada2021accdoa, wang2021four}, combinations of model-based localization with learning-based SED \cite{Nguyen2020c, Perez-Lopez2020}, data-based fusion of individual SSL and SED systems \cite{Nguyen2020c, nguyen2021general}, and event- or track-based prediction modeling, instead of class-based prediction \cite{Cao2020, Perez-Lopez2020}. The latter development specifically tried to address the case of same-class events occurring simultaneously \cite{nguyen2020sequence, Cao2020, cao2021improved}, a case that distinguishes the SELD task from SED and becomes possible mainly due to spatial information. Research following the DCASE2020 challenge investigated fusion of pre-trained SED and SSL models \cite{nguyen2021general}, or parameter sharing between joint, semi-joint SELD models, and models fusing SSL and SED subsystems \cite{cao2021improved}.

This report introduces the new \textbf{TAU-NIGENS Spatial Sound Events 2021}\footnote{\url{https://doi.org/10.5281/zenodo.4844825}} dataset and the baseline\footnote{\url{https://github.com/sharathadavanne/seld-dcase2021}} of the SELD challenge task in DCASE2021\footnote{\url{http://dcase.community/challenge2021/task-sound-event-localization-and-detection}}. The major difference of this dataset with the previous one is the introduction of localized interfering events outside of the target classes. This condition, naturally encountered in a real environment, introduces new challenges to the task. Apart from the dataset and baseline description, we present an extensive evaluation of the baseline on different versions of the dataset with and without the presence of ambient noise, directional interferers, and reverberation.

\section{Dataset}
\label{sec:dataset}

Similarly to the dataset of the previous iteration, the current one consists of 800 one-minute spatial recordings, of which 600 constitute the development set of the dataset, and the other 200 the evaluation set. The recordings are sampled at 24kHz, and they are offered in two 4-channel spatial audio formats, the raw signals of a tetrahedral microphone array and first-order Ambisonics, abbreviated as MIC and FOA for the rest of the paper. Detailed descriptions of the formats in terms of their directional encoding properties can be found in the previous challenge dataset reports \cite{Adavanne2019a, Politis2020}.

\subsection{Sound events}

The sound event samples are sourced from the \emph{NIGENS general sound events database} \cite{trowitzsch2019nigens}, which consists of 14 classes of specific sound types, and an additional general one with disparate sounds not belonging to any of the other classes. We use the sounds in the 12 classes \emph{alarm, crying baby, crash, barking dog, female scream, female speech, footsteps, knocking on door, male scream, male speech, ringing, phone, piano} as target events, and the sounds in the classes \emph{running engine, burning fire} and the general class as directional interferers. This division results in about 500 distinct sound samples distributed across the target events of the dataset, and about 400 across the interfering events.

\subsection{Dataset synthesis}

The synthesis of the spatial sound reccordings are based on a collection of SRIRs acquired continuously along measurement trajectories inside 13 enclosures of Tampere University. The RIR collection and synthesis process is described in more detail in \cite{Politis2020}. We summarize briefly the acoustical properties of the dataset. SRIRs are extracted along the measurement trajectories with an approximate resolution of 1 degree, resulting on about 1184 to 6480 possible RIRs/DOAs per room, depending on the type (circular/linear) and number of measurement trajectories. Events added in a single recording can be static or moving. The source position for a static event is drawn randomly from the pool of SRIRs of a single room used in that recording, while moving events are synthesized for one of the measured trajectories in the room. Moving events are synthesized to have an approximate speed of 10$^\circ$/sec, 20$^\circ$/sec, or 40$^\circ$/sec, drawn randomly.  The dataset is split into 8 folds with distinct rooms and samples in each of them. Distinct rooms result in different reverberation conditions, and even though similar ranges of DOAs may occur between rooms, the source distance, DRR, and reverberation conditions are distinct between folds for a certain DOA.

The events are laid out in layers in each recording, with the total number of layers determining the maximum polyphony possible. The parameter determining the density of events per layer and, hence, the average per-frame polyphony is the total gap time distributed between events in each layer. A larger gap time results in fewer events per layer and a lower average polyphony, while a smaller gap time results in higher event density and average polyphony. The last event per layer is truncated to fit the total 1 minute duration. For the present dataset there are three layers of target events and an additional layer of interfering events, resulting in a total maximum polyphony of 4. In addition to the spatialized reverberant events, multichannel ambient noise that was collected in each room with the same recording setup as the SRIRs is truncated to 1 minute segments and added to the event mixtures. The noise is scaled to result in signal-to-noise ratios (SNRs) ranging from noiseless (30dB) to noisy (6dB) conditions, with respect to the total energy of the target events in the recording excluding silences. A depiction of the layering of events in one recording is shown in Fig.~\ref{fig:recording}.

\begin{figure}[t]
  \centering
  \centerline{\includegraphics[width=0.4\textwidth ,keepaspectratio]{task3_seld_recording_sketch.png}}
  \caption{A graphic depiction of an emulated recording in the dataset, with colored objects indicating target classes, gray objects indicating interferers and ambient noise, and arrows indicating moving events.}
  \label{fig:recording}
\end{figure}

\subsection{Differences with DCASE2020 task 3 dataset}

Even though the acoustical and synthesis characteristics of the new dataset are similar to the dataset of the previous DCASE2020 Challenge, the following differences make it more challenging:
\begin{enumerate}
    \item Directional interferers, out of the target classes of a detection system, are common in real conditions and they add to the challenge by forcing a strong joint modeling and training strategy that can learn to ignore them.
    \item The overall maximum polyphony is increased from 2 to 3 target events.
    \item The recordings are not anymore divided into recordings with no overlap (polyphony 1), and recordings with two simultaneous events (polyphony 2). Instead all recordings have the maximum level of polyphony, with all intermediate levels (from silence to 3 simultaneous target events + interference) varying during the duration of the recording. This choice reflects more natural recording conditions in a real dataset.
    \item Even though the dataset of DCASE2020 had instances of the same class occurring at the same time, such occurrences were fairly rare. In the present dataset, these occurences have been increased in order to give a clear advantage to systems that can resolve this difficult but realistic case.
\end{enumerate}

\section{Baseline}
\label{sec:baseline}

\begin{figure}[t]
  \centering
  \centerline{\includegraphics[width=0.49\textwidth ,keepaspectratio]{CRNN_SELD_DCASE2021.png}}
  \caption{Convolutional recurrent neural network with ACCDOA loss for SELD.}
  \label{fig:crnn}
\end{figure}

Similar to the previous iterations of the challenge, we adopt a modified version of SELDnet \cite{adavanne2018sound} as the baseline method, due to its conceptual simplicity. Its architecture remains a convolutional recurrent neural network (CRNN) receiving multichannel log-mel spectrograms as inputs, together with acoustic intensity vectors \cite{pulkki2017first} for the FOA dataset, and generalized cross-correlation (GCC-PHAT) sequences for the MIC dataset, added as extra channels. The baseline implementation extracts log-mel spectrograms in  64 mel-bands from 1024-point FFTs, using a 40 ms window and 20 ms hop length at 24kHz. The intensity vectors are similarly extracted for every FFT bin and aggregated in the same number of mel-bands as the spectrograms, while the GCC sequences are also truncated to the same number of lag values as the mel-bands, adopted from \cite{Cao2019}. More details on the architecture and features can be found in \cite{Politis2020}.

The only difference of the current SELDnet baseline with respect to the previous DCASE challenge iteration is the output format and the respective loss function. The original SELDnet architecture employs separate output branches for detection and localization, with as many classification outputs and as many localization regressors as the number of classes. In the current baseline, we adopt the \emph{activity-coupled cartesian direction of arrival} representation (ACCDOA) introduced in DCASE2020 Challenge by Shimada et al. \cite{shimada2021accdoa}, which unifies the SED and SSL losses into a single homogeneous regression loss, simplifying the overall architecture by removing the detection branch, while simultaneously improving its performance. Using the ACCDOA representation, the network receives a sequence of $T$ STFT frames of multichannel features and outputs $T/5\times 3$ Cartesian vector coordinates for each of the target classes, with the direction of each vector indicating DOA and the vector length indicating class activity probability. The reduction in temporal resolution is intended to match the 100 ms resolution of annotations in the challenge. A block diagram of the current baseline architecture is shown in Fig.~\ref{fig:crnn}.

\section{Evaluation}
\label{sec:evaluation}
The dataset and baseline are delivered to the challenge participants at the commencement of the challenge, along with the development set of the dataset consisting of the 6 first folds, while the last two folds are made available during the evaluation phase of the challenge. Participants are required to report results on the test set of the development set using the predefined split of Table~\ref{tab:splits}, so that conclusions can be drawn among the submissions on the same configuration. However, during the evaluation phase, participants have to report results on the evaluation dataset only, using the development dataset for training and validation in any way they see fit.

\begin{table}[h]
\centering
\renewcommand\thetable{1}
\caption{Evaluation setup}
\begin{tabular}{l|ccc}
 & \multicolumn{3}{c}{\textbf{Splits}} \\ \cline{2-4}
\textbf{Dataset} & \textbf{Training} & \textbf{Validation} & \textbf{Testing} \\ \hline
\textbf{Development}  & 1,2,3,4 & 5 & 6 \\
\end{tabular}
\label{tab:splits}
\end{table}

The submissions are evaluated using the same combination of joint detection/localization metrics studied in \cite{mesaros2019joint, politis2020overview} and introduced in DCASE2020. Closer to SED evaluation, the localization-dependent error rate ($ER_X$) and F1-score ($F_X$) express detection performance but they penalize correct detections that occur further from the reference than some threshold distance $X$. On the other hand, the class-dependent localization error ($LE_\mathrm{CD}$) and localization recall ($LR_\mathrm{CD}$) are inspired by classical localization metrics, but are computed for each class individually before being averaged. The $LE_\mathrm{CD}$ is a mean angular localization error after pairing the predicted DOAs to their closest reference DOAs, while $LR_\mathrm{CD}$ is a simple recall metric on the detected localized events without any spatial threshold. Since in the SELD case there can be multiple simultaneous references of the same class, the detection metrics are modified to consider multiple instances of the same class and penalize cases where, e.g., only one of the predictions belong to that class. For the exact formulation of the metrics the reader is referred to \cite{politis2020overview}. The submissions are first ranked for each of the four metrics individually, and the final rank of each system is determined by the sum of the four individual ranks.

\input{table_baselines}

\section{Results}
\label{sec:Results}

In order to evaluate the performance of the new baseline utilizing the ACCDOA loss, we compare it against the previous SELDnet baseline of DCASE2020,
on the development set of DCASE2020 and the current one. Table~\ref{tab:baseline-results} shows a clear improvement of the ACCDOA version in all metrics. Especially in the more challenging new dataset, the ACCDOA loss brings large gains in detection and improves localization accuracy by about 25\%. A significant decrease of performance for both methods is also observed from the DCASE2020 dataset to the DCASE2021 dataset. This suggests that the new dataset is more challenging, as intended.

To get a more detailed picture on the effect of the various components in the scene, namely reverberation, ambient noise, and directional interferers, we generate various versions of the dataset including those components in various combinations. More specifically, the \emph{targets}, \emph{targets+ambience}, \emph{targets+interferers}, \emph{targets+ambience+interferers} develop from the presence of targets only, to the inclusion of ambient noise or interferers separately, to the full dataset combining all components. Excluding the effect of reverberation is less straightforward due to the use of real SRIRs for the synthesis. In order to generate reverberation-free versions of the dataset, the sound events for each recording in the original dataset are spatialized with anechoic IRs of the same Eigenmike spherical microphone array used to capture the SRIRs. The anechoic array IRs are computed for the same measurement trajectories and DOAs as the measured SRIRs in each room, and stored in a similar data structure. Additionally, each IR is delayed and scaled according to the source distance of the respective measured SRIR, following an inverse distance law and a speed of sound of $c=343 m/sec$. Delaying and scaling ensures that the events between the reverberant and non-reverberant versions are approximately time-aligned and with comparable distance-dependent attenuation. The Eigenmike responses were measured in an equirectangular grid of $5^{\circ}$ azimuth and $5^{\circ}$ elevation in the large anechoic chamber of Aalto University, as described in \cite{tervo2015direction}. Since the DOAs in the reverberant dataset do not necessarily coincide with the measurement grid of the array, array response interpolation is performed to recover anechoic IRs at the DOAs of the measured SRIRs, based on a spherical harmonic expansion of the array steering vectors, as in \cite{politis2017comparing}.

The results are presented in Table~\ref{tab:baseline-results}. As expected, reverberation affects negatively all combinations, increasing error rates and decreasing F-scores and localization recall in a consistent manner between the same scenarios. Additionally, it decreases localization accuracy by $2^{\circ}$--$4^{\circ}$. Inclusion of the ambient noise has a small but noticeable effect when added to the targets, without interferers. The small effect may be due to the large range of possible positive SNRs (6--30dB) distributed uniformly across the recordings, with most of them having adequate SNR to be unaffected by the noise presence. Interestingly, together with directional interference, inclusion of ambient noise seems to even improve certain results slightly.  This may be due to potential regularization effects of noise and is worth further investigation.

The most detrimental effects happen with the inclusion of the directional interferers, proving that this challenging case will need to be taken into account for future SELD systems. Error rate $ER$ increases up to about 40\% in the non reverberant case for the FOA recordings, and up to about 33\% for the MIC recordings. Similarly, in reverberant scenarios, the $ER$ increases up to about 28\% for both FOA and MIC formats. F-scores decrease by up to 40\% on the FOA dataset and up to 50\% on the MIC dataset, for both anechoic and reverberant conditions. The localization recall ($LR$) also drops by about 30\% for both formats and both anechoic and reverberant conditions. Finally, localization errors increase by up to about $7^{\circ}$ in the case of FOA recordings and up to $10^{\circ}$ in the case of MIC recordings, for both anechoic and reverberant conditions. In general, the MIC dataset exhibits a worse performance than FOA in all cases. This fact may be attributed to the input features employed in the baseline for each format. GCC sequences for the MIC format may become very noisy in complex scenes with multiple simultaneous events, while the intensity vectors of the FOA format can potentially retain robustness due to their narrowband nature and sparsity of the event signals in the time-frequency domain.

\input{table_results}

\section{Conclusions}
\label{sec:conclusions}

In this report we describe the new dataset and baseline for the SELD task of the DCASE2021 challenge. The differences with the dataset of the previous iteration of the challenge are highlighted; namely, inclusion of directional interferers, higher polyphony, and higher number of multiple simultaneous same-class event occurences. The evaluation task setup is also described, including a predefined fixed split on the development data for straightforward comparison of the submissions. The new baseline adopts the recent ACCDOA SELD representation introduced in the previous challenge to improve its performance, and is evaluated in the testing split of the development dataset. The new dataset is shown to be significantly more challenging for both the previous and the new baseline according to all considered metrics. A detailed analysis of the new baseline on different versions of the dataset shows that between reverberation, ambient noise, and directional interferers, the latter has the most detrimental effect of the three by far, in all evaluation metrics.

\newpage
\bibliographystyle{IEEEtran}


\section{Introduction}
\label{sec:introduction}
\input{sections/10-intro}

\section{Related Work}
\label{sec:related}
\input{sections/20-related}

\section{Estimating Risk under Exposure Bias}
\label{sec:estimating-risk}
\input{sections/30-estimating-risk}

\section{Learning Propensities and Link Probabilities}
\label{sec:learning-estimators}
\input{sections/40-learning-estimators}

\section{Feedback Loops}
\label{sec:feedback-loops}
\input{sections/50-feedback-loops}

\section{Experiments}
\label{sec:experiments}
\input{sections/60-experiments}

\section{Conclusion}
\label{sec:conclusion}
\input{sections/70-conclusion}

\subsection{Semi-Synthetic Dataset} \label{sec:experiments-semi-synthetic}

\begin{table}
\caption{Evaluation metrics
on the test set of the semi-synthetic data computed against
known ground truth citation links.
}
\label{table:semi-synthetic-eval-metrics}
\begin{center}
\begin{small}
\begin{sc}
\begin{tabular}{lcccccr}
\toprule
Model & Prec. & Rec. & AUC & MAP \\
\midrule
\makecell[l]{No Prop.} & 67.24 & 54.81 & 84.45 & 41.87 \\
\makecell[l]{MLE} & 81.04 & 60.19 & 93.12 & 56.77 \\
\makecell[l]{$\widehat{R}_w$} & \textbf{83.28} & 63.73 & \textbf{96.42} & 56.96 \\
\makecell[l]{$\widehat{R}_\text{PU}$} & 82.16 & 63.07 & 94.28 & 58.01 \\
\makecell[l]{$\widehat{R}_\text{AP}$} & 83.01 & \textbf{65.54} & 95.38 & \textbf{59.90} \\
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{center}
\end{table}

Since we do not have ground truth exposure values in the MAG dataset,
we cannot know whether a paper was not cited
due to a lack of exposure or due to irrelevancy.
As a result, we construct a semi-synthetic dataset
with simulated propensity scores and link probabilities.
We use a subset of $41{,}600$ papers.
We generate train-test-validation splits
by taking a topological ordering of the nodes
and use the subgraph created
from the first $70\%$ for training,
next $10\%$ for validation,
and the remaining $20\%$ for testing.
We use the real text and FOS for each paper.
The simulated propensity matrix $\pi$
is a $19 \times 19$ matrix
with its diagonal and off-diagonal entries initialized
from $U(0.7, 1)$ and $U(0.1, 0.3)$, respectively,
where $U(.)$ is the uniform distribution.
The link probability is simulated using
$y_{ij} = \sigma(w^\top (h_i \odot h_j) + b)$,
where $\sigma$ is the sigmoid function,
$\odot$ is element-wise product,
and $w, b$ are fixed \emph{known} vectors.

\begin{table}
\caption{RMSE of the estimated risk
with respect to the true risk
computed using our proposed estimators.
The first column shows the risk
used in the loss function in Eq.~\ref{eq:combined-loss-function}
to learn $\widehat{\pi}$ and $\widehat{y}$.
}
\label{table:semi-synthetic-rmse-metrics}
\begin{center}
\begin{small}
\begin{sc}
\begin{tabular}{l|ccccr}
\toprule
\multirow{2}{*}{\makecell{Trained \\ Using}} & \multicolumn{4}{c}{Estimator Used} \\
& $\widehat{R}_\text{naive}$ & $\widehat{R}_w$ & $\widehat{R}_\text{PU}$ & $\widehat{R}_\text{AP}$ \\
\midrule
No Prop. & 1.50 & - & - & - \\
\makecell[l]{MLE} & 0.67 & \textbf{0.23} & 0.24 & 0.32 \\
\makecell[l]{$\widehat{R}_w$} & 0.43 & \textbf{0.04} & 0.10 & 0.11 \\
\makecell[l]{$\widehat{R}_\text{PU}$} & 0.38 & 0.05 & 0.11 & \textbf{0.04} \\
\makecell[l]{$\widehat{R}_\text{AP}$} & 0.41 & 0.06 & 0.08 & \textbf{0.03} \\
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{center}
\end{table}

We show
the evaluation metrics for five models on the test set
computed against the simulated true citations
(not the observed citations) (Table \ref{table:semi-synthetic-eval-metrics}).
\textit{No Prop} is the model trained naively on the observed data using only
the output model in Eq. \ref{eq:experiments-output-model}.
\textit{MLE} is the model trained using the loss function
in Eq. \ref{eq:combined-loss-function} with $\lambda_R = 0$.
The remaining three are models trained using
$\widehat{R}_w, \widehat{R}_{\text{PU}}$,
and $\widehat{R}_{\text{AP}}$ with $\lambda_R = 10$ and $\lambda_L = 1$.
We see that all other estimators
significantly outperform \textit{No Prop}.
Additionally, our proposed estimators
lead to improved performance over the MLE.
We emphasize that these metrics
are computed against true citations
and thus are a measure of true risk
which is the appropriate metric
for evaluating an RS's performance.
This shows the utility of accounting for exposure bias
and learning using our proposed loss function.

\begin{figure}[t]
\begin{center}
\begin{subfigure}{0.47\columnwidth}
\includegraphics[width=0.9\columnwidth]{figures/ground_truth_propensities.pdf}
\caption{Ground Truth}
\end{subfigure}
\begin{subfigure}{0.47\columnwidth}
\includegraphics[width=0.9\columnwidth]{figures/estimated_propensities.pdf}
\caption{Estimated}
\end{subfigure}
\caption{The estimated propensities propensities are close to the true simulated values when learned using $\widehat{R}_w$.}
\label{fig:ground-truth-and-estimated-propensities}
\end{center}
\end{figure}

In this work, we tackle two separate (but related) challenges.
The first challenge is \textit{learning} link probabilities in such a way that they are not underestimated due to exposure bias.
The second challenge is \textit{evaluating} a RS given learned link probabilities and propensity scores, i.e., computing a good estimate of the true risk.
We demonstrate the efficacy of our methods for the second challenge and show
that our proposed weighting schemes lead to good estimates of the true risk (Table~\ref{table:semi-synthetic-rmse-metrics}).
We show
the RMSE of the risk estimated
using the proposed estimators
with respect to the true risk.
The first column denotes risk function
used to train the model
(as described in Section~\ref{sec:learning-estimators}).
The rest of the columns denote
the estimators used to estimate true risk
using the learned propensities and link probabilities
from the trained model in the first column.
We trained each model $10$ times to compute the RMSE.
The RMSE estimated using $\widehat{R}_{\text{naive}}$
is always greater than that of the other estimators,
which shows that leveraging the learned propensities
leads to substantially better estimates of the true risk
(and thus more accurately evaluates the RS).
The RMSE when trained using the \textit{MLE}
is higher than when trained using our proposed estimators,
showing the benefit of our proposed estimators over the \textit{MLE}.
This also qualitatively validates the generalization bound
proved in Section~\ref{sec:learning-estimators}
by showing that minimizing
$\widehat{R} \in \{ \widehat{R}_w, \widehat{R}_{\text{PU}}, \widehat{R}_{\text{AP}} \}$
also leads to small values of the true risk.

The heatmap of simulated propensities
and estimated propensities when using $\widehat{R}_w$
shows that the estimated propensities
are close to the true propensities (Figure \ref{fig:ground-truth-and-estimated-propensities}).
The mean relative error between the true
and estimated propensities is $19.47\%$,
demonstrating that the training procedure recovers the propensities.
Together, these results show that our methods
successfully mitigate exposure bias in this dataset.

\subsection{Real-World Datasets}\label{sec:real-world-datasets}

\begin{table}
\caption{Evaluation metrics for various models
computed on the test sets
of the two real-world citation datasets.}
\label{table:real-data-standard-metrics}
\begin{center}
\begin{small}
\begin{sc}
\begin{tabular}{lccccccr}
\toprule
Model & Prec. & Rec. & F1 & AUC & \makecell[l]{MAP}\\
\midrule
\multicolumn{6}{l}{\textbf{Dataset 1}} \\
\midrule
\makecell[l]{No Prop.} & 29.45 & 78.30 & 42.81 & 84.44 & 24.10 \\
\makecell[l]{MLE} & 30.24 & 77.84 & 43.56 & 84.41 & 24.60 \\
\makecell[l]{$\widehat{R}_w$} & 31.46 & 78.02 & 44.84 & 84.74 & 25.60  \\
\makecell[l]{$\widehat{R}_\text{PU}$} & 30.98 & \textbf{78.94} & 44.49 & \textbf{85.24} & 25.11 \\
\makecell[l]{$\widehat{R}_\text{AP}$} & \textbf{36.07} & 76.08 & \textbf{48.94} & 84.67 & \textbf{28.58} \\
\midrule
\multicolumn{6}{l}{\textbf{Dataset 2}} \\
\midrule
\makecell[l]{No Prop.} & 44.86 & 70.85 & 54.94 & 83.22 & 33.19 \\
\makecell[l]{MLE} & 44.43 & 74.66 & 55.71 & 84.97 & 34.39 \\
\makecell[l]{$\widehat{R}_w$} & \textbf{48.70} & 71.62 & \textbf{57.98} & 83.90 & \textbf{36.25} \\
\makecell[l]{$\widehat{R}_\text{PU}$} & 42.17 & \textbf{76.15} & 54.28 & \textbf{85.43} & 33.26 \\
\makecell[l]{$\widehat{R}_\text{AP}$} & 47.22 & 71.84 & 56.98 & 83.89 & 35.27 \\
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{center}
\end{table}
We now evaluate our proposed method
on a real-world citation network.
We construct two datasets
by using disjoint subgraphs of the MAG.
The first generated dataset has
$2{,}442{,}008$ papers and $7{,}577{,}886$ edges.
The second dataset has $1{,}328{,}664$ papers and $1{,}469{,}899$ edges.
Thus the second graph is sparser than the first one.
The FOS distribution is also different in both datasets
(see details in Appendix \ref{appendix:experiments}).
We use 70-10-20\% train-validation-test splits
generated similarly to the semi-synthetic dataset.
We do not have access to true exposure values
and thus we evaluate our methods against the observed citation links.

We show
the evaluation metrics for the proposed estimators (Table \ref{table:real-data-standard-metrics}).
Since we do not have access
to the true citation links in the real dataset,
we compute these metrics over the observed links.
In other words, this is a measure of the naive risk.
We see that our proposed estimators
achieve comparable metrics
to \textit{No Prop}.
For both datasets, the best numbers for each metric
are achieved by estimators other than \textit{No Prop}.
Moreover, the models using the weighted estimators
outperform the \textit{MLE} estimator in both datasets.
These results show that our proposed estimators
achieve comparable performance
even when evaluated on the observed citation data.
Similarly, Table \ref{table:real-data-link-pred-metrics} shows
link prediction metrics
for various models computed against \textit{observed} citations.
\textit{Recall@100} refers to the recall in the top 100 recommendations
averaged across all papers in the test set.
\textit{Mean Rank} is the mean rank
of the cited papers
averaged across all the papers.
\textit{Entropy@100 of True Positives} is the entropy in the FOS
of the true positives in the top 100 recommendations for each paper;
we use it to measure the diversity in the FOS of the recommendations.
Our proposed estimators achieve comparable \textit{Recall@100} and \textit{Mean Rank} to \textit{No Prop}
for both datasets.
As expected, propensity based estimators have higher FOS entropy scores than \textit{No Prop}, with $\widehat{R}_w$ achieving the highest FOS entropy in both datasets. Thus our proposed estimators \emph{recommend more relevant papers from different FOS and still maintain comparable performance} to \textit{No Prop}.

At first blush, the comparable performance to \textit{No Prop} may not seem compelling. However, this is a strong result.
Our goal is to correct exposure bias and minimize \textit{true risk}, not observed (or naive) risk.
Since Tables \ref{table:real-data-standard-metrics} and \ref{table:real-data-link-pred-metrics} are computed against observed citations, our proposed methods should not be expected to outperform \textit{No Prop} as they are not trying to optimize metrics against the observed links.
In Section \ref{sec:experiments-semi-synthetic}, we showed that our methods correct exposure bias and achieve lower true risk. Coupled with those results, our goal in this section was to show that our methods do not lower performance even if evaluated against the standard evaluation metrics.
\textcolor{black}{We suspect that the negative log-likelihood term $\mathcal{L}(o|\widehat{\pi}, \widehat{y})$ in Eq.~\ref{eq:combined-loss-function} is likely responsible for the comparable performance against the observed risk. This is because, as seen from the results, the MLE also performs well as compared to \textit{No Prop}.}.

\begin{table}
\caption{Link prediction metrics for various models when evaluated on the test set of the real datasets.}
\label{table:real-data-link-pred-metrics}
\begin{center}
\begin{small}
\begin{sc}
\begin{tabular}{lcccccr}
\toprule
Model & \makecell{Recall\\@100} & \makecell{Mean \\ Rank} & \makecell{Entropy@100 \\ True positives}\\
\midrule
\multicolumn{4}{l}{\textbf{Dataset 1}} \\
\midrule
\makecell[l]{No Prop.} & 24.39 & \textbf{2247.27} & 1.65 \\
\makecell[l]{MLE} & 24.70 & 2891.40 & 1.73 \\
\makecell[l]{$\widehat{R}_w$} & 25.03 & 2836.73 & \textbf{1.74} \\
\makecell[l]{$\widehat{R}_\text{PU}$} & 24.61 & 2875.13 & 1.73 \\
\makecell[l]{$\widehat{R}_\text{AP}$} & \textbf{26.66} & 2425.51 & 1.71 \\
\midrule
\multicolumn{4}{l}{\textbf{Dataset 2}} \\
\midrule
\makecell[l]{No Prop.} & \textbf{6.32} & \textbf{10170.26} & 1.06 \\
\makecell[l]{MLE} & 6.07 & 10731.88 & 1.08 \\
\makecell[l]{$\widehat{R}_w$} & 6.30 & 10873.19 & \textbf{1.12} \\
\makecell[l]{$\widehat{R}_\text{PU}$} & 5.92 & 10717.06 & 1.08 \\
\makecell[l]{$\widehat{R}_\text{AP}$} & 5.99 & 10801.11 & 1.10 \\
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{center}
\end{table}

\subsection{Feedback Loops}

\begin{figure}[t]
\begin{center}
\begin{subfigure}{0.49\columnwidth}
\includegraphics[width=1\columnwidth]{figures/feedback_loops_no_propensity.pdf}
\caption{No Propensity}
\label{fig:feedback-loops-fos-no-propensity}
\end{subfigure}
\begin{subfigure}{0.49\columnwidth}
\includegraphics[width=1\columnwidth]{figures/feedback_loops_with_propensity.pdf}
\caption{With Propensity ($\widehat{R}_w$)}
\label{fig:feedback-loops-fos-with-propensity}
\end{subfigure}
\caption{The fraction of recommended papers from the same FOS over time.}
\label{fig:feedback-loops-fos}
\end{center}
\end{figure}

We run simulations to examine what happens
when a citation recommender is trained repeatedly
on data collected from users
interacting with its recommendations.
We use the iterative training procedure
described in Section~\ref{sec:feedback-loops}.
We construct a training set of $410$ papers from the MAG,
with their corresponding real FOS and text embeddings.
For ease of exposition, we use an arbitrary
but fixed mapping to map the 19 FOS to two FOS.
The synthetic propensities and link probabilities
are simulated similarly to Sec~\ref{sec:experiments-semi-synthetic}.
In the first iteration, the models are trained
using the observed citation network.
For subsequent training iterations,
the training data is generated as follows.
For each paper, we recommend 20 papers.
The probability of recommending a paper is proportional
to its estimated citation probability.
We then simulate the user's interaction with the recommendations
according to the known simulated propensities and citation probabilities.
This generates the training set for the subsequent iteration.
We then repeat this process.

We show how the fraction of recommended papers
from the same FOS changes over multiple training iterations
for models trained without propensity, i.e., \textit{No Prop}
and the model trained using $\widehat{R}_w$
(Figure \ref{fig:feedback-loops-fos}).
We plot this time series for both FOS in our dataset.
For \textit{No Prop}, the fraction of papers recommended
from the same FOS increases over time for both FOS (Figure~\ref{fig:feedback-loops-fos-no-propensity}).
This demonstrates the existence of a feedback loop
that worsens exposure bias and reduces
the number of papers recommended
from a different FOS over time.
On the other hand, when we train our models using $\widehat{R}_w$,
the feedback loop no longer exists
and the fraction of papers recommended
from a different FOS
remains stable over time
(Figure~\ref{fig:feedback-loops-fos-with-propensity}).
This shows that our proposed estimator continues
to recommend relevant papers from a different FOS
and corrects the feedback loop.

\section{Bias and Variance} \label{appendix:bias-variance}

\begin{lemma}\label{lemma:variance-bernoulli}
Let $X \sim \text{Bernoulli}(\theta)$ and $Y = aX + b(1-X)$, where $a$ and $b$ are some constants. Then
\begin{align*}
    \text{Var}(Y) = \theta (1 - \theta) (a - b)^2.
\end{align*}
\end{lemma}

\textbf{Proof of Lemma \ref{lemma:bias-var-naive}}

\begin{lemma*}
The bias and variance of $\widehat{R}_{\text{naive}}(\widehat{o})$ are
\begin{align*}
    \text{B}(\widehat{R}_{\text{naive}}) &= \left|{\bf E}[\widehat{R}_{\text{naive}}] - R(\widehat{o})\right| \\
        &= \frac{\Delta}{|\mathcal{U}|} \left| \sum_{(i,j) \in \mathcal{U}} y_{ij} (1 - \pi_{ij}) (1 - 2 \widehat{o}_{ij}) \right|,\\
    \text{Var}(\widehat{R}_{\text{naive}}) &= \frac{\Delta^2}{|\mathcal{U}|^2} \sum_{(i,j) \in \mathcal{U}} y_{ij}\pi_{ij} (1 - y_{ij}\pi_{ij}).
\end{align*}
\end{lemma*}
\begin{proof}
We have
\begin{align*}
    \widehat{R}_{\text{naive}}(o, \widehat{y}) &= \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \delta(o_{ij}, \widehat{o}_{ij}) \\
        &= \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \left[ o_{ij} \delta(1, \widehat{o}_{ij}) + (1-o_{ij}) \delta(0, \widehat{o}_{ij}) \right] \\
    \therefore\,\, {\bf E}_o[\widehat{R}_{\text{naive}}(o, \widehat{y})] &= \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \left[ y_{ij} \pi_{ij} \delta(1, \widehat{o}_{ij}) + (1-y_{ij} \pi_{ij}) \delta(0, \widehat{o}_{ij}) \right] \\
    &= \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \left[ y_{ij} \pi_{ij} (1 - \widehat{o}_{ij})\delta(1, 0) + (1-y_{ij} \pi_{ij}) \widehat{o}_{ij} \delta(0, 1) \right] \\
    &= \frac{\Delta}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \left[ y_{ij} \pi_{ij} (1 - \widehat{o}_{ij}) + (1-y_{ij} \pi_{ij}) \widehat{o}_{ij} \right].
\end{align*}
The true risk is
\begin{align*}
    R(\widehat{y}) &= \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \left[y_{ij} \delta(1, \widehat{o}_{ij}) + (1 - y_{ij}) \delta(0, \widehat{o}_{ij}) \right] \\
    &= \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \left[y_{ij} (1-\widehat{o}_{ij}) \delta(1, 0) + (1 - y_{ij}) \widehat{o}_{ij} \delta(0, 1) \right] \\
    &= \frac{\Delta}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \left[y_{ij} (1-\widehat{o}_{ij}) + (1 - y_{ij}) \widehat{o}_{ij} \right].
\end{align*}
Thus the bias is
\begin{align*}
    \text{B}(\widehat{R}_{\text{naive}}) &= \left|{\bf E}[\widehat{R}_{\text{naive}}] - R(\widehat{o})\right| \\
    &= \left| \frac{\Delta}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \left[ y_{ij} \pi_{ij} (1 - \widehat{o}_{ij}) + (1-y_{ij} \pi_{ij}) \widehat{o}_{ij} - y_{ij} (1-\widehat{o}_{ij}) - (1 - y_{ij}) \widehat{o}_{ij} \right]  \right| \\
    &= \frac{\Delta}{|\mathcal{U}|} \left| \sum_{(i,j) \in \mathcal{U}} y_{ij} (1 - \pi_{ij}) (1 - 2 \widehat{o}_{ij}) \right|.
\end{align*}
The variance is
\begin{align*}
    \text{Var}(\widehat{R}_{\text{naive}}) &= \text{Var}\left( \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \left[ o_{ij} \delta(1, \widehat{o}_{ij}) + (1-o_{ij}) \delta(0, \widehat{o}_{ij}) \right] \right) \\
    &= \frac{1}{|\mathcal{U}|^2} \sum_{(i,j) \in \mathcal{U}} \text{Var}\left(  o_{ij} \delta(1, \widehat{o}_{ij}) + (1-o_{ij}) \delta(0, \widehat{o}_{ij}) \right) \\
    &= \frac{1}{|\mathcal{U}|^2} \sum_{(i,j) \in \mathcal{U}} y_{ij} \pi_{ij} (1 - y_{ij} \pi_{ij}) \left(\delta(1, \widehat{o}_{ij}) - \delta(0, \widehat{o}_{ij}) \right)^2 \,\,\, \text{(using Lemma \ref{lemma:variance-bernoulli})} \\
    &= \frac{\Delta^2}{|\mathcal{U}|^2} \sum_{(i,j) \in \mathcal{U}} y_{ij} \pi_{ij} (1 - y_{ij} \pi_{ij}).
\end{align*}
\end{proof}
Lemmas \ref{lemma:rw_bias_var}, \ref{lemma:bias-var-pu}, and \ref{lemma:bias-var-ap} can be proved similarly.
\newline

\textbf{Proof of Theorem \ref{thm:comp_var}}

\begin{theorem*}[\textbf{Comparison of Variances}]
For all values of $\widehat{\pi}, \widehat{y}$, we have $\text{Var}(\widehat{R}_{\text{AP}}) < \text{Var}(\widehat{R}_{\text{naive}}), \,\, \text{and} \\
\text{Var}(\widehat{R}_{\text{AP}}) < \text{Var}(\widehat{R}_w) < \text{Var}(\widehat{R}_{\text{PU}})$
\end{theorem*}
\begin{proof}
First we show that $\text{Var}(\widehat{R}_{\text{AP}}) < \text{Var}(\widehat{R}_{\text{naive}})$. We have
\begin{align*}
    \text{Var}(\widehat{R}_{\text{AP}}) &= \frac{\Delta^2}{|\mathcal{U}|^2} \sum_{(i,j) \in \mathcal{U}} y_{ij}\pi_{ij} (1 - y_{ij}\pi_{ij}) \psi^2_{ij}, \\
    \text{where} \,\, & \psi_{ij} = \frac{1-\widehat{y}_{ij}}{1-\widehat{\pi}_{ij}\widehat{y}_{ij}} < 1, \\
    \text{Var}(\widehat{R}_{\text{naive}}) &= \frac{\Delta^2}{|\mathcal{U}|^2} \sum_{(i,j) \in \mathcal{U}} y_{ij}\pi_{ij} (1 - y_{ij}\pi_{ij}).
\end{align*}
Using the fact that $\psi_{ij}^2 < 1 \,\, \forall (i,j) \in \mathcal{U}$, we get $\text{Var}(\widehat{R}_{\text{AP}}) < \text{Var}(\widehat{R}_{\text{naive}})$.

Next, we show that $\text{Var}(\widehat{R}_w) < \text{Var}(\widehat{R}_{\text{PU}})$:
\begin{align*}
    \text{Var}(\widehat{R}_w) &= \frac{\Delta^2}{|\mathcal{U}|^2} \sum_{(i,j) \in \mathcal{U}} y_{ij}\pi_{ij} (1 - y_{ij}\pi_{ij}) \left( \frac{1 - \widehat{o}_{ij}}{\widehat{\pi}^2_{ij}} + \widehat{o}_{ij} \psi^2_{ij} \right) \\
    \text{Var}(\widehat{R}_{\text{PU}}) &= \frac{\Delta^2}{|\mathcal{U}|^2} \sum_{(i,j) \in \mathcal{U}} \frac{y_{ij}\pi_{ij} (1 - y_{ij}\pi_{ij})}{\widehat{\pi}^2_{ij}} \\
        &= \frac{\Delta^2}{|\mathcal{U}|^2} \sum_{(i,j) \in \mathcal{U}} y_{ij}\pi_{ij} (1 - y_{ij}\pi_{ij}) \left( \frac{1 - \widehat{o}_{ij}}{\widehat{\pi}^2_{ij}} + \frac{\widehat{o}_{ij}}{\widehat{\pi}^2_{ij}}  \right) \\
    \therefore\,\, \text{Var}(\widehat{R}_w) - \text{Var}(\widehat{R}_{\text{PU}}) &= \frac{\Delta^2}{|\mathcal{U}|^2} \sum_{(i,j) \in \mathcal{U}} y_{ij}\pi_{ij} (1 - y_{ij}\pi_{ij}) \widehat{o}_{ij} \left( \psi^2_{ij} - \frac{1}{\widehat{\pi}^2_{ij}} \right) \\
    \therefore\,\, \text{Var}(\widehat{R}_w) - \text{Var}(\widehat{R}_{\text{PU}}) &< 0 \,\,\, \left(\text{because $\psi_{ij} < 1$ and $\frac{1}{\widehat{\pi}_{ij}} > 1$}\right) \\
    \therefore\,\, \text{Var}(\widehat{R}_w) &< \text{Var}(\widehat{R}_{\text{PU}}).
\end{align*}
Next, we show that $\text{Var}(\widehat{R}_{\text{AP}}) < \text{Var}(\widehat{R}_w)$:
\begin{align*}
    \text{Var}(\widehat{R}_{\text{AP}}) - \text{Var}(\widehat{R}_w) &= \frac{\Delta^2}{|\mathcal{U}|^2} \sum_{(i,j) \in \mathcal{U}} y_{ij}\pi_{ij} (1 - y_{ij}\pi_{ij}) (1 - \widehat{o}_{ij}) \left( \psi^2_{ij} - \frac{1}{\widehat{\pi}^2_{ij}} \right) \\
    \therefore\,\, \text{Var}(\widehat{R}_{\text{AP}}) - \text{Var}(\widehat{R}_w) &< 0 \,\,\, \left(\text{because $\psi_{ij} < 1$ and $\frac{1}{\widehat{\pi}_{ij}} > 1$}\right) \\
    \therefore\,\, \text{Var}(\widehat{R}_{\text{AP}}) &< \text{Var}(\widehat{R}_w).
\end{align*}
\end{proof}

\textbf{Proof of Theorem \ref{thm:bias-comparisons}}

\begin{theorem*}[\textbf{Comparison of Biases}]
Under the bias approximations, a sufficient condition
for $\text{B}(\widehat{R}_w) = \text{B}(\widehat{R}_{\text{PU}}) < \text{B}(\widehat{R}_{\text{naive}})$ is
\begin{align*}
    \frac{\pi_{ij}}{2 - \pi_{ij}} < \widehat{\pi}_{ij} < 1, \,\, \forall (i, j) \in \mathcal{U},
\end{align*}
and for $\text{B}(\widehat{R}_{\text{AP}}) < \text{B}(\widehat{R}_{\text{naive}})$ is
\begin{align*}
    & \frac{\pi_{ij}}{2 - \pi_{ij}} < \widehat{\pi}_{ij} < 1 \,\, \text{and} \,\, 0 < \widehat{y}_{ij} < c y_{ij}, \,\, \forall (i,j) \in \mathcal{U} \\
    & \text{where}\,\, c = \frac{2 (1 - \pi_{ij})}{1 - \widehat{\pi}_{ij} - \pi_{ij} y_{ij} + (2 - \pi_{ij}) \widehat{\pi}_{ij} y_{ij}} \geq 1.
\end{align*}
\end{theorem*}
\begin{proof}
We first derive the sufficient condition for $\text{B}(\widehat{R}_w) = \text{B}(\widehat{R}_{\text{PU}}) < \text{B}(\widehat{R}_{\text{naive}})$. We have
\begin{align*}
    \text{B}(\widehat{R}_{\text{naive}}) &\approx \frac{\Delta}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}'} y_{ij} (1 - \pi_{ij}), \\
    \text{B}(\widehat{R}_w) \approx \text{B}(\widehat{R}_{\text{PU}}) &\approx \begin{aligned}[t]
        & \frac{\Delta}{|\mathcal{U}|} \left| \sum_{(i,j) \in \mathcal{U}'} y_{ij}\left(1 - \frac{\pi_{ij}}{\widehat{\pi}_{ij}} \right) \right|.
        \end{aligned}
\end{align*}
If $1 > \widehat{\pi}_{ij} > \pi_{ij} \forall (i,j) \in \mathcal{U}$, we have
\begin{align*}
    & \left(1 - \frac{\pi_{ij}}{\widehat{\pi}_{ij}} \right) > 0 \,\,\, \forall (i,j) \in \mathcal{U} \\
    \therefore\,\, & \text{B}(\widehat{R}_w) \approx \text{B}(\widehat{R}_{\text{PU}}) \approx \begin{aligned}[t]
        & \frac{\Delta}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}'} y_{ij}\left(1 - \frac{\pi_{ij}}{\widehat{\pi}_{ij}} \right).
        \end{aligned} \\
    \therefore\,\, & \begin{aligned}[t]
    \text{B}(\widehat{R}_w) - \text{B}(\widehat{R}_{\text{naive}}) &= \frac{\Delta}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}'} y_{ij} \left( \pi_{ij} - \frac{\pi_{ij}}{\widehat{\pi}_{ij}} \right) \\
        &= \frac{\Delta}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}'} y_{ij} \pi_{ij} \left( 1 - \frac{1}{\widehat{\pi}_{ij}} \right) \\
        &< 0 \,\,\, \text{(because $\widehat{\pi}_{ij} < 1$)}
    \end{aligned} \\
    \therefore\,\, & \text{B}(\widehat{R}_w) < \text{B}(\widehat{R}_{\text{naive}}).
\end{align*}
If $0 < \widehat{\pi}_{ij} \leq \pi_{ij} \forall (i,j) \in \mathcal{U}$, we have
\begin{align*}
    & \left(1 - \frac{\pi_{ij}}{\widehat{\pi}_{ij}} \right) \leq 0 \,\,\, \forall (i,j) \in \mathcal{U} \\
    \therefore\,\, & \text{B}(\widehat{R}_w) \approx \text{B}(\widehat{R}_{\text{PU}}) \approx \begin{aligned}[t]
        & \frac{\Delta}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}'} y_{ij}\left( \frac{\pi_{ij}}{\widehat{\pi}_{ij}} - 1 \right).
        \end{aligned}.
\end{align*}
Then, a sufficient condition for $\text{B}(\widehat{R}_w) = \text{B}(\widehat{R}_{\text{PU}}) < \text{B}(\widehat{R}_{\text{naive}})$ is
\begin{align*}
    y_{ij}\left( \frac{\pi_{ij}}{\widehat{\pi}_{ij}} - 1 \right) &< y_{ij} (1 - \pi_{ij}) \,\,\, \forall (i,j) \in \mathcal{U} \\
    \therefore\,\, y_{ij}\left( \frac{\pi_{ij}}{\widehat{\pi}_{ij}} - 1 \right) &< y_{ij} (1 - \pi_{ij}) \,\,\, \forall (i,j) \in \mathcal{U} \\
    \therefore\,\, \widehat{\pi}_{ij} &> \frac{2}{2 - \pi_{ij}} \,\,\, \forall (i,j) \in \mathcal{U}.
\end{align*}
Next, we derive the sufficient condition for $\widehat{R}_{\text{AP}} < \widehat{R}_{\text{naive}}$. Observe that
\begin{align*}
    & \frac{\pi_{ij}}{2 - \pi_{ij}} < \widehat{\pi}_{ij} < 1 \,\,\, \forall (i, j) \in \mathcal{U} \\
    \therefore\,\, & (1 -\pi_{ij}) y_{ij} - (1-\pi_{ij}y_{ij}) \tau_{ij} \geq 0 \,\,\, \forall (i, j) \in \mathcal{U} \\
    \therefore\,\, & \widehat{R}_{\text{AP}} \approx \frac{\Delta}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}'} \left[ (1-\pi_{ij}) y_{ij} - (1-\pi_{ij}y_{ij}) \tau_{ij} \right], \,\, \text{where} \,\, \tau_{ij} = \left( \frac{\widehat{y}_{ij}(1 - \widehat{\pi}_{ij})}{1 - \widehat{\pi}_{ij} \widehat{y}_{ij}} \right).
\end{align*}
Therefore, when $\frac{\pi_{ij}}{2 - \pi_{ij}} < \widehat{\pi}_{ij} < 1 \,\,\, \forall (i, j) \in \mathcal{U}$, a sufficient condition for $\widehat{R}_{\text{AP}} < \widehat{R}_{\text{naive}}$ is
\begin{align*}
    (1 -\pi_{ij}) y_{ij} - (1-\pi_{ij}y_{ij}) \tau_{ij} < y_{ij} (1 - \pi_{ij})  \,\,\, \forall (i, j) \in \mathcal{U} \\
    \therefore\,\, 0 < \widehat{y}_{ij} < \left( \frac{2 (1 - \pi_{ij})}{1 - \widehat{\pi}_{ij} - \pi_{ij} y_{ij} + (2 - \pi_{ij}) \widehat{\pi}_{ij} y_{ij}} \right) y_{ij} \,\,\, \forall (i, j) \in \mathcal{U}.
\end{align*}
\end{proof}

\section{Generalization Bound} \label{appendix:generalization-bound}

\textbf{Proof of Theorem \ref{thm:generalization-bound}}

\begin{theorem*}[\textbf{Generalization Bound}]
Let $\mathcal{F}$ be a class of functions
$(\widehat{\pi}, \widehat{y})$.
Let $\delta(o_{ij}, \widehat{y}_{ij}) \leq \eta\,\,\forall (i,j) \in \mathcal{U}$
and $\widehat{\pi}_{ij} \geq \epsilon > 0\,\,\forall (i,j) \in \mathcal{U}$.
Then, for $\widehat{R} \in \left\{ \widehat{R}_w,  \widehat{R}_{\text{PU}}, \widehat{R}_{\text{AP}} \right\}$, with probability at least $1 - \delta$, we have
\begin{align}
    R(\widehat{y}) &\leq \widehat{R}(\widehat{y}, \widehat{\pi}) + B(\widehat{R}) + 2 \mathcal{G}(\mathcal{F}, \widehat{R}) + M \label{eq:rademacher-bound} \\
        &\leq \widehat{R}(\widehat{y}, \widehat{\pi}) + B(\widehat{R}_w) + 2 \widehat{\mathcal{G}}(\mathcal{F}, \widehat{R}_w) + 3 M, \label{eq:rademacher-empirical-bound}
\end{align}
where $M = \sqrt{\frac{4 \eta^2}{\epsilon^2 |\mathcal{U}|} \log(\frac{2}{\delta})}$ and $B(\widehat{R})$ is the bias of $\widehat{R}$ derived in Section~\ref{sec:estimating-risk}.
\end{theorem*}
\begin{proof}
We proceed similarly to the standard Rademacher complexity generalization bound
proof \citep{shalev2014understanding}[Ch. 26]. Observe that
\begin{align}
    R(\widehat{y}) &= R(\widehat{y}) - {\bf E}_o[\widehat{R}(o, \widehat{y}, \widehat{\pi})] + {\bf E}_o[\widehat{R}(o, \widehat{y}, \widehat{\pi})] \nonumber \\
        &\leq B(\widehat{R}) + {\bf E}_o[\widehat{R}(o, \widehat{y}, \widehat{\pi})]. \label{eq:rademacher-proof-1}
\end{align}
Let $\Phi(o) = \sup_{(\widehat{\pi}, \widehat{y}) \in \mathcal{F}} \left[ {\bf E}_o[\widehat{R}(o, \widehat{y}, \widehat{\pi})] - \widehat{R}(o, \widehat{y}, \widehat{\pi}) \right]$. Then
\begin{align}\label{eq:rademacher-proof-2}
     {\bf E}_o[\widehat{R}(o, \widehat{y}, \widehat{\pi})]| \leq \widehat{R}(o, \widehat{y}, \widehat{\pi}) + \Phi(o).
\end{align}
Now we upper bound $\Phi(o)$.
Since $\delta(o_{ij}, \widehat{y}_{ij}) \leq \eta\,\,\forall (i,j)$
and $\widehat{\pi}_{ij} \geq \epsilon > 0,\,\,\forall (i,j)$ and $\forall \, \widehat{R} \in \left\{ \widehat{R}_w,  \widehat{R}_{\text{PU}}, \widehat{R}_{\text{AP}} \right\}$,
we have
\begin{align*}
    |\Phi(o) - \Phi(\Tilde{o})| \leq \frac{2 \eta}{\epsilon},
\end{align*}
if $o$ and $\Tilde{o}$ differ in only one coordinate, i.e.,
$o_{ij} \neq \Tilde{o}_{ij}$ for some $(i,j) \in \mathcal{U}$
and $o_{lm} = \Tilde{o}_{lm} \forall (l,m) \in \mathcal{U}$
s.t. $(i, j) \neq (l,m)$.
Using McDiarmid’s Inequality, with probability at least $1 - \delta$,
we have
\begin{align} \label{eq:rademacher-proof-3}
    \Phi(o) \leq {\bf E}[\Phi(o)] + C.
\end{align}
Next, we upper bound ${\bf E}[\Phi(o)]$.
Let $\bar{o}$ be a ghost sample independently drawn having the same distribution as $o$. We have
\begin{align}
    {\bf E}[\Phi(o)] &= {\bf E}_o\left[ \sup_{(\widehat{\pi}, \widehat{y}) \in \mathcal{F}} \left[ {\bf E}_o[\widehat{R}(o, \widehat{y}, \widehat{\pi})] - \widehat{R}(o, \widehat{y}, \widehat{\pi}) \right] \right] \nonumber \\
    &= {\bf E}_o\left[ \sup_{(\widehat{\pi}, \widehat{y}) \in \mathcal{F}} {\bf E}_{\bar{o}}\left[\widehat{R}(\bar{o}, \widehat{y}, \widehat{\pi}) - \widehat{R}(o, \widehat{y}, \widehat{\pi}) \,\big| \, o \right] \right] \nonumber \\
    &= {\bf E}_o\left[ \sup_{(\widehat{\pi}, \widehat{y}) \in \mathcal{F}} {\bf E}_{\bar{o}}\left[ \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} r(\bar{o}_{ij}, \widehat{\pi}_{ij}, \widehat{y}_{ij}) - \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} r(o_{ij}, \widehat{\pi}_{ij}, \widehat{y}_{ij}) \,\Big| \, o \right] \right] \nonumber \\
    &\leq {\bf E}_{o, \bar{o}} \left[ \sup_{(\widehat{\pi}, \widehat{y}) \in \mathcal{F}} \left[ \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} r(\bar{o}_{ij}, \widehat{\pi}_{ij}, \widehat{y}_{ij}) - \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} r(o_{ij}, \widehat{\pi}_{ij}, \widehat{y}_{ij}) \right] \right] \,\, \text{(Jensen's Inequality)} \nonumber \\
    &= {\bf E}_{o, \bar{o}, \sigma} \left[ \sup_{(\widehat{\pi}, \widehat{y}) \in \mathcal{F}} \left[ \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \sigma_{ij} r(\bar{o}_{ij}, \widehat{\pi}_{ij}, \widehat{y}_{ij}) - \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \sigma_{ij} r(o_{ij}, \widehat{\pi}_{ij}, \widehat{y}_{ij}) \right] \right] \nonumber \\
    &= {\bf E}_{o, \bar{o}, \sigma} \left[ \sup_{(\widehat{\pi}, \widehat{y}) \in \mathcal{F}} \left[ \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \sigma_{ij} r(\bar{o}_{ij}, \widehat{\pi}_{ij}, \widehat{y}_{ij}) + \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \sigma_{ij} r(o_{ij}, \widehat{\pi}_{ij}, \widehat{y}_{ij}) \right] \right] \nonumber \\
    &\leq {\bf E}_{o, \bar{o}, \sigma} \left[ \sup_{(\widehat{\pi}, \widehat{y}) \in \mathcal{F}} \left[ \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \sigma_{ij} r(\bar{o}_{ij}, \widehat{\pi}_{ij}, \widehat{y}_{ij}) \right] + \sup_{(\widehat{\pi}, \widehat{y}) \in \mathcal{F}} \left[ \frac{1}{|\mathcal{U}|} \sum_{(i,j) \in \mathcal{U}} \sigma_{ij} r(o_{ij}, \widehat{\pi}_{ij}, \widehat{y}_{ij}) \right] \right] \nonumber \\
    &= 2 \mathcal{G}(\mathcal{F}, \widehat{R}). \label{eq:rademacher-proof-4}
\end{align}
Combining Eqs. \ref{eq:rademacher-proof-1},
\ref{eq:rademacher-proof-2}, \ref{eq:rademacher-proof-3},
and \ref{eq:rademacher-proof-4}, we get Eq. \ref{eq:rademacher-bound}.
Another application of McDiarmid’s Inequality
allows us to obtain Eq. \ref{eq:rademacher-empirical-bound}
from Eq. \ref{eq:rademacher-bound}.
\end{proof}

\section{Feedback Loops} \label{appendix:feedback-loops}

\begin{lemma}[\textbf{Binomial Tail Bound}]\label{lemma:appendix-binomial-tail}
If the random variable $X_n \sim \frac{1}{n} \text{Binomial}(n, \theta)$, then for $\epsilon > 0$, we have
\begin{align*}
    {\bf P}(|X_n - \theta| > \epsilon) \leq 2 \exp\left( -2n \epsilon^2 \right).
\end{align*}
\begin{proof}
Observe that $X_n \in [0, 1]$. Applying Hoeffding's inequality gives us the desired result.
\end{proof}
\end{lemma}

\begin{lemma}\label{lemma:appendix-thm-proof-lemma-1}
Let $n \in \mathbb{N}$ and $\kappa$ be a fixed $C - 1$ simplex such that $\kappa_v n \in \mathbb{N} \,\,\, \forall v \in [C]$. The random variable $\Tilde{q}_v \sim \frac{1}{\kappa_v n} \text{Binomial}(\kappa_v n, q_v)$, where $q_v \in (0, 1)$. Assume that $q_v > q_w$ if $v > w$. We denote as $\widehat{e}$ the following $C-1$ simplex:
\begin{align*}
    \widehat{e} = \frac{1}{Z} \left[ \kappa_1 \Tilde{q}_1, \kappa_2 \Tilde{q}_2, \hdots, \kappa_C \Tilde{q}_C \right], \,\,\, \text{where} \,\, Z = \sum_{i \in [C]} \kappa_i \Tilde{q}_i.
\end{align*}
Let $\widehat{e}_{vw} = \frac{\widehat{e}_v}{\widehat{e}_v + \widehat{e}_w} = \frac{\kappa_v \Tilde{q}_v}{\kappa_w \Tilde{q}_w + \kappa_w \Tilde{q}_w}$ and $\kappa_{vw} = \frac{\kappa_v}{\kappa_v + \kappa_w}$. Then for a constant $\rho_{vw}$ such that
\begin{align*}
    0 < \rho_{vw} < \frac{\kappa_v \kappa_w (q_v - q_w)}{q_v \kappa^2_v + (q_v+q_w)\kappa_v \kappa_w + q_w \kappa^2_w},
\end{align*}
we have
\begin{align*}
    & |\Tilde{q}_v - q_v| < \epsilon_{vw}, \,\, |\Tilde{q}_w - q_w| < \epsilon_{vw} \implies \widehat{e}_{vw} - \kappa_{vw} > \rho_{vw}, \\
    & \text{for some constant $\epsilon_{vw}$ s.t.} \,\, 0 < \epsilon_{vw} < \frac{\rho_{vw} q_v \kappa^2_v - \kappa_v \kappa_w (q_w - q_v) + q_w \rho_{vw} \kappa_v (\kappa_v - \kappa_w)}{\rho_{vw} (\kappa^2_v - \kappa^2_w ) - 2 \kappa_v \kappa_w }.
\end{align*}
This is saying that, for $(v, w)$ s.t. $v > w$ the simplex $\widehat{e}$ will be more skewed towards $v$ than the simplex $\kappa$ if the sampled $\Tilde{q}_v$ and $\Tilde{q}_w$ are close to their mean values $q_v$ and $q_w$, respectively.
\end{lemma}
\begin{proof}
Observe that if $|\Tilde{q}_v - q_v| < \epsilon_{vw}$ and $|\Tilde{q}_w - q_w| < \epsilon_{vw}$, then the lowest value that $\widehat{e}_{vw}$ can take is
\begin{align*}
    & \widehat{e}^{(\text{min})}_{vw} = \frac{\kappa_v (q_v - \epsilon_{vw})}{\kappa_v (q_v - \epsilon_{vw}) + \kappa_w (q_w + \epsilon_{vw})}, \,\,\, \text{and} \\
    & \widehat{e}^{(\text{min})}_{vw} - \kappa_{vw} > \rho_{vw} \implies \widehat{e}_{vw} - \kappa_{vw} > \rho_{vw}.
\end{align*}
Therefore, we have
\begin{align*}
    & \widehat{e}^{(\text{min})}_{vw} - \kappa_{vw} > \rho_{vw} \,\,\, \text{and} \,\,\, \epsilon_{vw} < q_w \\
    & \impliedby \underbrace{\frac{\kappa_v (q_v - \epsilon_{vw})}{\kappa_v (q_v - \epsilon_{vw}) + \kappa_w (q_w + \epsilon_{vw})} - \frac{\kappa_v}{\kappa_v + \kappa_w} > \rho_{vw}}_{(1)} \,\,\, \text{and} \,\,\, \rho_{vw} < \frac{\kappa_v \kappa_w (q_v - q_w)}{q_v \kappa^2_v + (q_v+q_w)\kappa_v \kappa_w + q_w \kappa^2_w}.
\end{align*}
The inequality (1) above can further be simplified as
\begin{align*}
    & \frac{\kappa_v (q_v - \epsilon_{vw})}{\kappa_v (q_v - \epsilon_{vw}) + \kappa_w (q_w + \epsilon_{vw})} - \frac{\kappa_v}{\kappa_v + \kappa_w} > \rho_{vw} \\
    \impliedby & \epsilon_{vw} < \frac{\rho_{vw} q_v \kappa^2_v - \kappa_v \kappa_w (q_w - q_v) + q_w \rho_{vw} \kappa_v (\kappa_v - \kappa_w)}{\rho_{vw} (\kappa^2_v - \kappa^2_w ) - 2 \kappa_v \kappa_w }.
\end{align*}
This completes the proof.
\end{proof}

\begin{lemma}\label{lemma:appendix-thm-proof-lemma-2}
Let $\alpha$ be a fixed $C-1$ simplex and $\widehat{e}$ be the following $G-1$ simplex, $\widehat{e} = \frac{1}{Z} [\alpha_1 \Tilde{q}_1, \alpha_2 \Tilde{q}_2, \hdots, \alpha_C \Tilde{q}_C]$, where $Z = \sum_{z \in [C]} \alpha_z \Tilde{q}_z$ and the vector $\kappa \sim \frac{1}{n} \text{Multinomial}(n, \widehat{e})$. Let $\widehat{e}_{vw} = \frac{\widehat{e}_v}{\widehat{e}_v + \widehat{e}_w} = \frac{\Tilde{q}_v}{\Tilde{q}_w + \Tilde{q}_w}$ and $\kappa_{vw} = \frac{\kappa_v}{\kappa_v + \kappa_w}$.

Assume that $|\Tilde{q}_z - q_z| < \epsilon \,\,\, \forall z \in [C]$ where $q_z \in (0, 1)$ are fixed. If $|\kappa_v - \widehat{e}_v < \frac{\eta_{nw}}{C}|$ and $|\kappa_w - \widehat{e}_w < \frac{\eta_{nw}}{C}|$, then for some constant $\rho$, we have
\begin{align*}
    \widehat{e}_{vw} - \kappa_{vw} < \rho, \,\,\,
    \text{when} \,\,\, \eta_{vw} < \rho \left( \frac{q_v + q_w}{\max_{z \in [C]}q_z + \epsilon} \right).
\end{align*}
\begin{proof}
If $|\kappa_v - \widehat{e}_v < \frac{\eta_{nw}}{C}|$ and $|\kappa_w - \widehat{e}_w < \frac{\eta_{nw}}{C}|$, then the smallest value that $\kappa_{vw}$ can achieve is
\begin{align*}
    \kappa^{\text{(min)}}_{vw} = \frac{\widehat{e}_v - \frac{\eta_{nw}}{C}}{\widehat{e}_v + \widehat{e}_w}.
\end{align*}
This means that
\begin{align*}
    & \widehat{e}_{vw} - \kappa_{vw} < \rho \\
    \impliedby & \widehat{e}_{vw} - \kappa^{\text{(min)}}_{vw} < \rho \\
    \iff & \frac{\frac{\eta_{nw}}{C}}{\widehat{e}_v + \widehat{e}_w} < \rho \\
    \iff & \frac{\eta_{nw}}{C} < \rho (\widehat{e}_v + \widehat{e}_w).
\end{align*}
Since $|\Tilde{q}_z - q_z| < \epsilon \,\,\, \forall z \in [C]$, we have
\begin{align*}
    \widehat{e}_v &= \frac{\alpha_v \Tilde{q}_v}{\sum_{z \in [C]} \alpha_z \Tilde{q}_z} \\
        &> \frac{\alpha_v (q_v - \epsilon)}{\alpha_v (q_v - \epsilon) + \sum_{z \in [C], z \neq v} \alpha_z (q_z + \epsilon)} \\
        &> \frac{\alpha_v (q_v - \epsilon)}{C \max_{z \in [C]} \alpha_z (q_z + \epsilon)}, \\
    \text{and similarly} \,\,\, \widehat{e}_w &> \frac{\alpha_w (q_w - \epsilon)}{C \max_{z \in [C]} \alpha_w (q_z + \epsilon)} \\
    \therefore\,\, \widehat{e}_v + \widehat{e}_w &> \frac{ \alpha_v (q_v - \epsilon) + \alpha_w (q_w - \epsilon) }{C \max_{z \in [C]} (q_z + \epsilon)}.
\end{align*}
Therefore, we can set $\eta_{vw}$ such that
\begin{align*}
    \frac{\eta_{nw}}{C} &< \rho \left(\frac{\alpha_v (q_v - \epsilon) + \alpha_w (q_w - \epsilon) }{C \max_{z \in [C]} (q_z + \epsilon)} \right) \\
    \therefore\,\, \eta_{nw} &< \rho \left( \frac{\alpha_v (q_v - \epsilon) + \alpha_w (q_w - \epsilon) }{\max_{z \in [C]} q_z + \epsilon} \right).
\end{align*}
\end{proof}
\end{lemma}

\textbf{Proof of Theorem \ref{thm:feedback-loop-finite-sample}}

\begin{theorem*}
Suppose that $q_v > q_w$ if $v > w$.
Let $\kappa^{(t)}_{vw} = \frac{\kappa^{(t)}_v}{\kappa^{(t)}_v + \kappa^{(t)}_w}$.
Let $A^{(t)}_{vw}$ represent the event
that relative fraction of recommendations
from $c_v$ to that from $c_w$ increases at time $t$,
i.e., $\kappa^{(t+1)}_{vw} > \kappa^{(t)}_{vw}$.
Let $A^{(t)}$ be the event
that all relative fractions get skewed
towards $c_v$ from $c_w$ if $q_v > q_w$,
i.e., $A^{(t)} = \bigcap_{(v,w) \in \mathcal{S}} A^{(t)}_{vw}$,
where $\mathcal{S} = \{ (v, w): v \in [C], w \in [C], v > w \}$.
Then, for constants $\epsilon, \eta > 0$
that only depend on $\kappa^{(t)}$ and $q$, we have
\begin{align*}
    & {\bf P}(A^{(t)} |\kappa^{(t)}) \begin{aligned}[t]
        & \geq 1 - 2C \exp\left( -2 n \left[ \epsilon^2 + \frac{\eta^2}{C^2} \right]  \right) \\
        & \geq 1 - 2C \exp\left( - \mathcal{O}\left(\frac{n}{C^2}\right) \right).
    \end{aligned}
\end{align*}
\end{theorem*}
\begin{proof}
We know that the estimated probabilities $\widehat{q}^{(t)}_v$ have distribution $\widehat{q}^{(t)}_v | \kappa^{(t)} \sim \frac{1}{n} \text{Binomial}(n \kappa^{(t)}_v, q_v)$. The simplex with normalized probabilities is $\widehat{e}^{(t+1)} = \frac{1}{Z} [\widehat{q}^{(t)}_1, \widehat{q}^{(t)}_2, \hdots, \widehat{q}^{(t)}_C]$, where $Z = \sum_{z \in [C]} \widehat{q}^{(t)}_z$.

Let $\Tilde{q}^{(t)}_v = \frac{\widehat{q}^{(t)}_v}{\kappa^{(t)}_v}$.
Observe that $\Tilde{q}^{(t)}_v | \kappa^{(t)} \sim \frac{1}{n \kappa^{(t)}_v} \text{Binomial}(n \kappa^{(t)}_v, q_v)$. We denote by $\widehat{e}^{(t+1)}_{vw}$,
\begin{align*}
    \widehat{e}^{(t+1)}_{vw} = \frac{\widehat{e}^{(t+1)}_v}{\widehat{e}^{(t+1)}_v + \widehat{e}^{(t+1)}_w} = \frac{\kappa^{(t)}_v \Tilde{q}^{(t)}_v}{\kappa^{(t)}_v \Tilde{q}^{(t)}_v + \kappa^{(t)}_w \Tilde{q}^{(t)}_w}.
\end{align*}
There are two main parts to the proof.
First, we show that, with high probability,
$\widehat{e}^{(t+1)}_{vw} - \kappa^{(t)}_{vw} > \rho \,\, \forall (v,w) \in \mathcal{S}$ for some constant $\rho$.
Then, we show that, with high probability,
$\widehat{e}^{(t+1)}_{vw} - \kappa^{(t+1)}_{vw} < \rho \,\, \forall (v,w) \in \mathcal{S}$. We combine these two results to show that, with high probability, $\kappa^{(t+1)}_{vw} > \kappa^{(t)}_{vw} \,\, \forall (v,w) \in \mathcal{S}$.

Using Lemma \ref{lemma:appendix-thm-proof-lemma-1}, for some $(v, w) \in \mathcal{S}$, we know that for some constant $\rho_{vw}$ such that
\begin{align*}
    0 < \rho_{vw} < \frac{\kappa^{(t)}_v \kappa^{(t)}_w (q_v - q_w)}{q_v (\kappa^{(t)}_v)^2 + (q_v+q_w)\kappa^{(t)}_v \kappa^{(t)}_w + q_w (\kappa^{(t)}_w)^2},
\end{align*}
we have
\begin{align*}
    & |\Tilde{q}^{(t)}_v - q_v| \leq \epsilon_{vw} \,\, \textit{and} \,\, |\Tilde{q}^{(t)}_w - q_w| \leq \epsilon_{vw} \implies \widehat{e}^{(t+1)}_{vw} - \kappa^{(t)}_{vw} \geq \rho_{vw}, \\
    & \text{for a constant $\epsilon_{vw}$ s.t.} \,\, 0 < \epsilon_{vw} < \frac{\rho_{vw} q_v (\kappa^{(t)}_v)^2 - \kappa^{(t)}_v \kappa^{(t)}_w (q_w - q_v) + q_w \rho_{vw} \kappa^{(t)}_v (\kappa^{(t)}_v - \kappa^{(t)}_w)}{\rho_{vw} ((\kappa^{(t)}_v)^2 - (\kappa^{(t)}_w)^2 ) - 2 \kappa^{(t)}_v \kappa^{(t)}_w } \\
    \implies & {\bf P}\left(\widehat{e}^{(t+1)}_{vw} - \kappa^{(t)}_{vw} \geq \rho_{vw}\right) \geq {\bf P}\left(|\Tilde{q}^{(t)}_v - q_v| \leq \epsilon_{vw}, |\Tilde{q}^{(t)}_w - q_w| \leq \epsilon_{vw}\right).
\end{align*}
Intuitively, this is saying
that $\widehat{e}^{(t+1)}_{vw} - \kappa^{(t)}_{vw} > \rho_{vw}$ if $\Tilde{q}^{(t)}_v$ and $\Tilde{q}^{(t)}_w$ are close to $q_v$ and $q_w$, respectively.
Let $\rho = \min_{(v,w) \in \mathcal{S}} \rho_{vw}$ and $\epsilon = \min_{(v,w) \in \mathcal{S}} \epsilon_{vw}$. Then we have
\begin{align}
    {\bf P} \left( \bigcap_{(v,w) \in \mathcal{S}} \widehat{e}^{(t+1)}_{vw} - \kappa^{(t)}_{vw} \geq \rho \right) &\geq {\bf P} \left( \bigcap_{z \in [C]} |\Tilde{q}^{(t)}_z - q_z| \leq \epsilon \right) \\
    &= 1 - {\bf P} \left( \bigcup_{z \in [C]} |\Tilde{q}^{(t)}_z - q_z| \geq \epsilon \right) \\
    &\geq 1 - \sum_{z=1}^{C} {\bf P} \left( |\Tilde{q}^{(t)}_z - q_z| \geq \epsilon \right) \,\,\, \text{(Union Bound)} \\
    &\geq 1 - \sum_{z=1}^{C} 2 \exp\left( -2n \epsilon^2 \right)
    \,\,\, \text{(using Lemma \ref{lemma:appendix-binomial-tail})} \nonumber \\
    &= 1 - 2C \exp\left( -2n \epsilon^2 \right). \label{eq:feedback-finite-sample-proof-1}
\end{align}
Now, we show that $\widehat{e}^{(t+1)}_{vw}$
is close to $\kappa^{(t+1)}_{vw}$.
We know that $\kappa^{(t+1)} \sim \frac{1}{n} \text{Multinomial}(n, \widehat{e}^{(t+1)})$.
Let the event $Q^{(t)} = \bigcap_{z \in [C]} |\Tilde{q}^{(t)}_z - q_z| \leq \epsilon$.
Using Lemma \ref{lemma:appendix-thm-proof-lemma-2}, we know that, under $Q^{(t)}$, for some constant $\eta_{vw}$, we have
\begin{align*}
    & \left|\widehat{e}^{(t+1)}_v - \kappa^{(t+1)}_v \right| < \frac{\eta_{vw}}{C} \,\,\, \text{and} \,\,\, \left|\widehat{e}^{(t+1)}_w - \kappa^{(t+1)}_w \right| < \frac{\eta_{vw}}{C} \implies \widehat{e}^{(t+1)}_{vw} - \kappa^{(t+1)}_{vw} < \rho, \\
    & \text{where} \,\, 0 < \eta_{vw} < \frac{ \kappa^{(t+1)}_v (q_v - \epsilon) + \kappa^{(t+1)}_w (q_w - \epsilon)}{\max_{z \in [C]} \kappa^{(t+1)}_z (q_z + \epsilon)} \\
    \implies & {\bf P} \left( \widehat{e}^{(t+1)}_{vw} - \kappa^{(t+1)}_{vw} < \rho \, \big| \, Q^{(t)} \right) \geq {\bf P}\left(\left|\widehat{e}^{(t+1)}_v - \kappa^{(t+1)}_v \right| < \frac{\eta_{vw}}{C}, \,\,\, \left|\widehat{e}^{(t+1)}_w - \kappa^{(t+1)}_w \right| < \frac{\eta_{vw}}{C} \right).
\end{align*}
Intuitively, this is saying that $\widehat{e}^{(t+1)}_{vw} - \kappa^{(t+1)}_{vw} < \rho$ if $\kappa^{(t+1)}_v$ and $\kappa^{(t+1)}_w$ are close to $\widehat{e}^{(t+1)}_v$ and $\widehat{e}^{(t+1)}_w$, respectively. Thus, for $\eta = \min_{(v,w) \in \mathcal{S}} \eta_{vw}$, we have
\begin{align}
    {\bf P} \left( \bigcap_{(v,w) \in \mathcal{S}} \widehat{e}^{(t+1)}_{vw} - \kappa^{(t+1)}_{vw} \leq \rho \, \Big| \, Q^{(t)} \right) &\geq {\bf P}\left( \bigcap_{z \in [C]} |\widehat{e}^{(t+1)}_z - \kappa^{(t+1)}_z| \leq \frac{\eta}{C} \right) \nonumber \\
    &= 1 - {\bf P}\left( \bigcup_{z \in [C]} |\widehat{e}^{(t+1)}_z - \kappa^{(t+1)}_z| \geq \frac{\eta}{C} \right) \nonumber \\
    &\geq 1 - \sum_{z=1}^{C} {\bf P}\left(|\widehat{e}^{(t+1)}_z - \kappa^{(t+1)}_z| > \frac{\eta}{C} \right) \,\,\, \text{(Union Bound)} \nonumber \\
    &\geq 1 - 2C \exp\left( -\frac{2n \eta^2}{C^2} \right) \,\,\, \text{(using Lemma \ref{lemma:appendix-binomial-tail})}. \label{eq:feedback-finite-sample-proof-2}
\end{align}
Combining Eq. \ref{eq:feedback-finite-sample-proof-1} and \ref{eq:feedback-finite-sample-proof-2},
we get the desired result as follows:
\begin{align*}
    {\bf P} \left( \bigcap_{(v,w) \in \mathcal{S}} A^{(t)}_{vw} \right) &= {\bf P} \left( \bigcap_{(v,w) \in \mathcal{S}} \kappa^{(t+1)}_{vw} > \kappa^{(t)}_{vw} \right) \\
    &\geq {\bf P} \left( \bigcap_{(v,w) \in \mathcal{S}} \left( \widehat{e}^{(t+1)}_{vw} - \kappa^{(t)}_{vw} \geq \rho, \,\,\, \widehat{e}^{(t+1)}_{vw} - \kappa^{(t+1)}_{vw} \leq \rho  \right) \right) \\
    &\geq {\bf P}\left( \bigcap_{z \in [C]} \left( |\widehat{e}^{(t+1)}_z - \kappa^{(t+1)}_z| \leq \frac{\eta}{C}, \,\,\, |\Tilde{q}^{(t)}_z - q_z| \leq \epsilon \right) \right) \\
    &= {\bf P}\left( \bigcap_{z \in [C]} |\widehat{e}^{(t+1)}_z - \kappa^{(t+1)}_z| \leq \frac{\eta}{C} \, \Big| \, Q^{(t)} \right) {\bf P}\left( \bigcap_{z \in [C]}  |\Tilde{q}^{(t)}_z - q_z| \leq \epsilon \right) \\
    &\geq \left( 1 - 2C \exp\left( -\frac{2n \eta^2}{C^2} \right) \right) \left( 1 - 2C \exp\left( -2n \epsilon^2 \right) \right) \\
    &\geq 1 - 2C \left[ \exp\left( -2n \epsilon^2 \right) + \exp\left( -\frac{2n \eta^2}{C^2} \right) \right] \\
    &\geq 1 - 2C \exp\left( - \mathcal{O}\left(\frac{n}{C^2}\right) \right).
\end{align*}
\end{proof}

\textbf{Proof of Theorem \ref{thm:worsen_rate}}

\begin{lemma}[Convergence in Probability]\label{lemma:appendix-convergence-in-prob}
Let $X_n, Y_n$, and $Z$ be random variables such that $X_n \overset{p}{\to} Y_n$ and $Y_n \overset{p}{\to} Z$, then $X_n \overset{p}{\to} Z$.
\begin{proof}
For any $\epsilon > 0$, we have
\begin{align*}
    {\bf P}(|X_n - Z| \geq \epsilon) &= {\bf P}(|X_n - Y_n + Y_n - Z| \geq \epsilon) \\
    &\leq {\bf P}(|X_n - Y_n| + |Y_n-Z| \geq \epsilon) \\
    &\leq {\bf P}\left(|X_n - Y_n| \geq \frac{\epsilon}{2}\right) + {\bf P}\left(|X_n - Y_n| \geq \frac{\epsilon}{2}\right) \\
    &= 0.
\end{align*}
Therefore, $X_n \overset{p}{\to} Z$.
\end{proof}
\end{lemma}

\begin{theorem*}
Suppose that $q_v > q_w$. As $n \rightarrow \infty$, $\kappa^{(t)}_{vw} \overset{p}{\rightarrow} 1 - \frac{1}{1 + c^t}$, where $c = \frac{q_v}{q_w}$.
\end{theorem*}
\begin{proof}
At time step $t$, the fraction of recommendations
from each group is $\kappa_t$.
From group $g_v$, the user cites papers
according to probability $q_v$.
Therefore, $\widehat{q}^{(t)}_v \overset{p}{\rightarrow} \kappa^{(t)}_v q_v$.
And the normalized estimate is
$\widehat{e}^{(t+1)} = \frac{1}{S} [\kappa^{(t)}_1 q_1, \hdots, \kappa^{(t)}_C q_C]$, where $S = \sum_{z \in [C]} \kappa^{(t)}_z q_z$.
Since $\kappa^{(t+1)} \sim \frac{1}{n} \text{Multinomial}(n, \widehat{e}^{(t+1)})$,
we have
\begin{align}
    \kappa^{(t+1)} &\overset{p}{\rightarrow} \widehat{e}^{(t+1)} \nonumber \\
    \frac{\kappa^{(t+1)}_v}{\kappa^{(t+1)}_w} &\overset{p}{\rightarrow} \frac{q_v \kappa^{(t)}_v}{q_w \kappa^{(t)}_w} \nonumber \\
    &= c \frac{ \kappa^{(t)}_v}{ \kappa^{(t)}_w}. \label{eq:appendix-feedback-asymptotic-proof}
\end{align}
We know that $\frac{ \kappa^{(1)}_v}{ \kappa^{(1)}_w} \overset{p}{\rightarrow} c$.
Combining this with Eq. \ref{eq:appendix-feedback-asymptotic-proof} and using Lemma \ref{lemma:appendix-convergence-in-prob} recursively, we get
\begin{align*}
    \frac{\kappa^{(t)}_v}{\kappa^{(t)}_w} &\overset{p}{\rightarrow} c^t \\
    \therefore\,\, 1 - \frac{1}{1 + \frac{\kappa^{(t)}_v}{\kappa^{(t)}_w}} &\overset{p}{\rightarrow} 1 - \frac{1}{1 + c^t} \,\,\, \text{(Continuous mapping theorem)} \\
    \therefore\,\, \frac{\kappa^{(t)}_v}{\kappa^{(t)}_v + \kappa^{(t)}_w} &\overset{p}{\rightarrow} 1 - \frac{1}{1 + c^t} \\
    \therefore\,\, \kappa^{(t)}_{vw} &\overset{p}{\rightarrow} 1 - \frac{1}{1 + c^t}.
\end{align*}
\end{proof}

\section{Experiments} \label{appendix:experiments}

\begin{table}
\caption{The distribution of the FOS in the two real-world datasets.
}
\label{table:fos-distributions-in-dataset}
\begin{center}
\begin{sc}
\begin{tabular}{lcc}
\toprule
FOS & Dataset 1 & Dataset 2 \\
\midrule
Art & $0.03\%$ & $0.08\%$ \\
Biology & $26.48\%$ & $23.43\%$ \\
Business & $0.38\%$ & $0.10\%$ \\
Chemistry & $10.11\%$ & $15.67\%$ \\
Computer Science & $9.40\%$ & $3.42\%$ \\
Economics & $2.51\%$ & $0.03\%$ \\
Engineering & $6.24\%$ & $17.98\%$ \\
Environmental Science & $0.13\%$ & $0.03\%$ \\
Geography & $0.48\%$ & $0.40\%$ \\
Geology & $1.45\%$ & $0.46\%$ \\
History & $0.04\%$ & $0.03\%$ \\
Materials Science & $3.06\%$ & $19.09\%$ \\
Mathematics & $7.17\%$ & $1.03\%$ \\
Medicine & $21.28\%$ & $13.90\%$ \\
Philosophy & $0.03\%$ & $0.01\%$ \\
Physics & $2.99\%$ & $3.14\%$ \\
Political Science & $0.18\%$ & $0.01\%$ \\
Psychology & $7.49\%$ & $1.14\%$ \\
Sociology & $0.55\%$ & $0.05\%$ \\
\bottomrule
\end{tabular}
\end{sc}
\end{center}
\end{table}

Table \ref{table:fos-distributions-in-dataset} provides the distribution of the various FOS in both the datasets used for the real-world dataset experiments (Section~\ref{sec:real-world-datasets}). We can see that the FOS distributions are different. For example, Dataset 2 has substantially more \textit{Materials Science} and \textit{Engineering} papers.