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

Recommender system provides personalized service for users to seek information, playing an increasingly important role in a wide range of online applications, such as e-commerce, news portal, content-sharing platform, and social media.
However, the system faces popularity bias issues, which stand on the opposite of personalization.
On one hand, the user-item interaction data usually exhibits long-tail distribution on item popularity --- a few head items occupy most of the interactions whereas the majority of items receive relatively little attention~\cite{FeedbackLoop, Himan_reg}.
On the other hand, the recommender model trained on such long-tail data not only inherits the bias, but also amplifies the bias, making the popular items dominate the top recommendations~\cite{UMAP_recbias,Himan-thesis}.
Worse still, the feedback loop ecology of recommender system further intensifies such Matthew effect~\cite{Mattheweffect}, causing notorious issues like echo chamber~\cite{EchoChambers} and filter bubble~\cite{Calibrated}.
Therefore, it is essential to consider the popularity bias issue in recommender systems.

Existing work on popularity bias-aware recommendation mainly performs unbiased learning or ranking adjustment, which can be categorized into:
\begin{itemize}[leftmargin=*]

    \item \textit{Inverse Propensity Scoring} (IPS), which adjusts the data distribution to be even by reweighting the interaction examples for model training~\cite{IPS_rec,rec_ips_cap_norm}. Although IPS methods have sound theoretical foundations, they hardly work well in practice due to the difficulties in estimating propensities and high model variance.

    \item \textit{Causal Embedding}, which uses bias-free uniform data to guide the model to learn unbiased embedding~\cite{2018causalE,distillation}, forcing the model to discard item popularity. However, obtaining such uniform data needs to randomly expose items to users, which has the risk of hurting user experience. Thus the data is usually of a small scale, making the learning less stable.

    \item \textit{Ranking Adjustment}, which performs post-hoc re-ranking on the recommendation list~\cite{xquad2019,BPR_PC} or model regularization on the training~\cite{Himan_reg,BPR_PC}.

    Both types of methods are heuristically designed to intentionally increase the scores of less popular items, which however lack theoretical foundations for effectiveness.

\end{itemize}

Instead of eliminating the effect of popularity bias, we argue that the recommender system should \textbf{leverage} the popularity bias. The consideration is that not all popularity biases in the data mean bad effect. For example,
some items demonstrate higher popularity because of better intrinsic quality or representing current trends, which deserve more recommendations.
Blindly eliminating the effect of popular bias will lose some important signals implied in the data, improperly suppressing the high-quality or fashionable items.
Moreover, some platforms have the need of introducing desired bias into the system, \emph{e.g., } promoting the items that have the potential to be popular in the future.
This work aims to fill the research gap of effectively leveraging popularity bias to improve the recommendation accuracy.

To understand how item popularity affects the recommendation process, we resort to the language of \textit{causal graph}~\cite{pearl2009causality} for a qualitative analysis first. Figure \ref{fig:rec_graph}(a) illustrates that traditional methods mainly perform user-item matching to predict the affinity score: $U$ (user node) and $I$ (item node) are the causes, and $C$ is the effect node to denote the interaction probability. An example is the prevalent latent factor model~\cite{mf,lightgcn}, which forms the prediction as the inner product between user embedding and item embedding. Since how a model forms the prediction implies how it assumes the labeled data be generated, this causal graph could also interpret the assumed generation process of the observed interaction data. Item popularity, although exerts a significant influence on data generation process, is not explicitly considered by such coarse-grained modeling methods.

We next consider how item popularity affects the process, enriching the causal graph to Figure~\ref{fig:rec_graph}(b). Let node $Z$ denote item popularity, which has two edges pointing to $C$ and $I$, respectively. First, $Z\rightarrow C$ means the item popularity directly affects the interaction probability, since many users have herd mentality (\textit{a.k.a.,} conformity), thus tend to follow the majority to consume popular items~\cite{conformity,DICE}. Second, $Z\rightarrow I$ means item popularity affects whether the item is exposed, since recommender system usually inherits the bias in the data and exposes popular items more frequently\footnote{Note that we assume users make interaction choice on the exposed items only, and do not consider other information seeking choices like search which is not the focus of this work. Thus the observed interactions are conditioned on the exposure of the interacted items before. }. Remarkably, we find $Z$ is the common cause for $I$ and $C$, acting as the \textit{confounder}~\cite{pearl2009causality} between the exposed items and the observed interactions.
It means that, item popularity $Z$ affects the observed interaction data with two causal paths: 1) $Z\rightarrow C$ and 2) $Z\rightarrow I \rightarrow C$, wherein the second path contains the bad effect of  bias amplification, since it increases observed interactions of popular items even though they may not well match user interest.

To remove the bad effect of popularity bias on model training, we need to intervene recommended items $I$ to make them immune to their popularity $Z$. Experimentally, it means we need to alter the exposure policy to make it free from the impact of item popularity and then recollect data, which is however costly and impossible to achieve for academia researchers. Thanks to the progress on causal science, we can achieve the same result with \textit{do-calculus}~\cite{pearl2009causality} without performing interventional experiments. In short, we estimate the user-item matching as $P(C| do(U, I))$ that cuts off the path $Z\rightarrow I$ during training, differing from the correlation $P(C | U, I)$ estimated by existing recommender models that confounds user interest with popularity bias. Through such deconfounded training, $P(C| do(U, I))$ estimates user interest matching on items more accurately than $P(C | U, I)$, removing the spurious correlations between $I$ and $C$ due to the confounder of $Z$.
During the inference stage, we infer the ranking score as $P(C| do(U, I), do(Z))$, intervening the item popularity $Z$ with our desired bias (\emph{e.g., } the forecasted popularity in the testing stage).

The main contributions of this work are summarized as follows:
\begin{itemize}[leftmargin=*]

    \item We analyze how popularity bias affects recommender system with causal graph, a powerful tool but is seldom used in the community; we identify that the bad impact of popularity bias stems from the confounding effect of item popularity.

    \item We propose \textit{Popularity-bias Deconfounding and Adjusting} (PDA), a novel framework that performs deconfounded training with \textit{do}-calculus and causally intervenes the popularity bias during recommendation inference.

    \item We conduct extensive experiments on three real datasets, validating the effectiveness of our causal analyses and methods
\end{itemize}

\begin{figure}
\centering
\subfigure[Causal graph of traditional methods.]{\label{fig:a}\includegraphics[width=0.15\textwidth]{figures/tradition_graph.pdf}}
\subfigure[Causal graph that considers item popularity.]{\label{fig:b} \includegraphics[width=0.14\textwidth]{figures/our_graph_pop.pdf}}
\subfigure[We cut off $Z\rightarrow I$ for model training.]{\label{fig:c}\includegraphics[width=0.13\textwidth]{figures/graph_cut_off_ZtoI.pdf}}
\vspace{-10pt}
\caption{Causal graphs to describe the recommendation process. U: user, I: exposed item, C: interaction probability, Z: item popularity.
We identify $Z$ as the confounder between $I$ and $C$, and propose deconfounded training with $P(C| do(U, I))$ as the interest matching.}
\vspace{-15pt}
\label{fig:rec_graph}
\end{figure}
\section{Primary knowledge}
We use uppercase character (\emph{e.g., } $U$) to denote a random variable and lowercase character (\emph{e.g., } $u$) to denote its specific value. We use characters in calligraphic font (\emph{e.g., } $\mathcal{U}$) to represent the sample space of the corresponding random variable, and use $P(\cdot)$ to represent probability distribution of a random variable.

Let $\mathcal{D}$ denote the historical data, which is sequentially collected through $T$ stages, \emph{i.e., } $\mathcal{D} = \{\mathcal{D}_1 \cup \cdots \cup \mathcal{D}_T\}$; $\mathcal{U}=\{u_1 , \dots, u_{|\mathcal{U}|}\}$ and $\mathcal{I}=\{i_1 , \dots, i_{{|\mathcal{I}|}}\}$ denote all users and items, respectively.
Through learning on historical data, the recommender system is expected to capture user preference and serves well for the next stage $T+1$. That is, it aims to obtain high recommendation accuracy on $\mathcal{D}_{T+1}$.
In this work, we are interested in the factor of item popularity, defining the \textit{local popularity} of item $i$ on the stage $t$ as:
 \begin{equation}\label{eq:pop_define}\small
     m_{i}^{t} = D_{i}^{t} / \sum_{j\in \mathcal{I}} D_{j}^{t},
 \end{equation}
where $D_{i}^{t}$ denotes the number of observed interactions for item $i$ in $\mathcal{D}_{t}$. We can similarly define the global popularity of an item $m_i$ based on its interaction frequency in $\mathcal{D}$, but we think the local popularity has a larger impact on the system's exposure mechanism and user decision, since the system is usually periodically retrained and the most recent data has the largest impact.

\paragraph{Popularity drift.}
Intuitively, item popularity is dynamic and changes over time, meaning that the impact of popularity bias could also be dynamic. To quantify the popularity drift,
we define a metric named \textit{Drift of Popularity} (DP) to measure the drift between two stages.
First, we represent each stage $t$ as a probability distribution over items: $\big[m_{1}^{t},\dots,m_{|\mathcal{I}|}^{t}\big]$, where each entry denotes an item's frequency in the stage.
Then, we use the Jensen-Shannon Divergence (JSD)~\cite{js-range} to measure the similarity between two stages:
\begin{equation}\small
    \begin{aligned}
    DP(t,s) = JSD\left(
        \big[m_{1}^{t},\dots,m_{|\mathcal{I}|}^{t}\big],
        \big[m_{1}^{s},\dots,m_{|\mathcal{I}|}^{s}\big]
    \right),
\end{aligned}
\end{equation}
where $t$ and $s$ are the two stages.
Similar to JSD, the range of DP is $[0,log(2)]$, and a higher value indicates a larger popularity drift.

Figure~\ref{fig:drit_showing}(a) shows the DP value of two successive stages, \emph{i.e., } $DP(t,t+1)$ where $t$ iterates from 1 to 9, on three real-world datasets (details in Section~\ref{exp:setting}). We can see that the popularity drift obviously exists in all three datasets, and different datasets exhibit different levels of popularity drift.
Figure~\ref{fig:drit_showing}(b) shows the DP value of the first and present stage, \emph{i.e., } $DP(1,t)$, which measures the accumulated popularity drift.
We can see a clear increasing trend, indicating that the longer the time interval is, the large popularity drift the data exhibits.
These results reveal that popularity bias and its impact also change with time. The popularity bias in the future stage differs from that of past stages. If we set the target of model generalization as pursuing high accuracy on the next stage data $\mathcal{D}_{T+1}$\footnote{To our knowledge, many industrial practitioners on recommender systems use the next stage data for model evaluation.}, then a viable way is to predict popularity trend and inject it into the recommendations.

\begin{figure}
    \centering

    \subfigure[$DP(t,t+1)$]{\includegraphics[width=0.23\textwidth]{figures/drift_neibor_type.pdf}}
    \subfigure[$DP(1,t)$]{\includegraphics[width=0.23\textwidth]{figures/drift_ref1_type.pdf}}
    \vspace{-15pt}

    \caption{Popularity drift between: (a) two successive stages $DP(t, t+1)$; (b) the first and present stage  $DP(1, t)$.}
    \vspace{-13pt}
    \label{fig:drit_showing}
\end{figure}
\section{Methodology}
We first analyze the impact of item popularity on recommendation from a fundamental view of causality. Then we present our PDA framework that eliminates the bad effect of popularity bias on training and intervenes the inference stage with desired bias. Lastly, we present in-depth analyses of our method.

\subsection{Causal View of Recommendation}\label{ssec:causal graph}
Figure \ref{fig:rec_graph}(b) shows the causal graph that accounts for item popularity in affecting recommendation. By definition,
causal graph~\cite{pearl2009causality} is a directed acyclic graph where a node denotes a variable and an edge denotes a causal relation between two nodes~\cite{pearl2009causality}.
It is widely used to describe the process of data generation, which can guide the design of predictive models. Next, we explain the rationality of this causal graph.

\begin{itemize}[leftmargin=*]
    \item Node $U$ represents the user node, \emph{e.g., } user profile or history feature that is used for representing a user.
    \item Node $I$ represents the item node. We assume a user can only interact with (\emph{e.g., } click) the items that have been exposed to him/her, thus $I$ denotes the exposed item.
    \item Node $C$ represents the interaction label, indicating whether the user has chosen/consumed the exposed item.
    \item Node $Z$ represents the item popularity, which can be seen as a hidden variable since it is usually not explicitly modeled by traditional methods (\emph{e.g., } \cite{mf,lightgcn} form predictions with Figure \ref{fig:rec_graph}(a)). But it has a large impact on the recommendation.
    \item Edges $\{U, I, Z\} \rightarrow C$ denote that an interaction label $C$ is determined by the three factors: user $U$, item $I$, and the item's popularity $Z$. Traditional methods only consider $\{U, I\}\rightarrow C$ which is easy to explain: the matching between user interest and item property determines whether a user consumes an item. Here we intentionally add a cause node $Z$ for capturing user conformity --- many users have the herd mentality and tend to follow the majority to consume popular items~\cite{conformity,DICE}. Thus, whether there is an interaction is the combined effect of $U, I$ and $Z$.
    \item Edge $Z\rightarrow I$ denotes that item popularity affects the exposure of items. This phenomenon has been verified on many recommender models, which are shown to favor popular items after training on biased data~\cite{Himan-thesis}.

\end{itemize}

From this causal graph, we find that item popularity $Z$ is a confounder that affects both exposed items $I$ and observed interactions $C$. This results in two causal paths starting from $Z$ that affect the observed interactions: $Z\rightarrow C$ and $Z\rightarrow I\rightarrow C$. The first path is combined with user-item matching to capture the user conformity effect, which is as expected. In contrast, the second path means that item popularity increases the exposure likelihood of popular items, making the observed interactions consist of popular items more. Such effect causes bias amplification, which should be avoided since an authentic recommender model should estimate user preference reliably and is immune to the exposure mechanism.

In the above text, we have explained the causal graph from the view of data generation. In fact, it also makes sense to interpret it from the view of model prediction. Here, we denote $U$ and $I$ as user embedding and item embedding, respectively, and traditional models apply inner product~\cite{mf,lightgcn} or neural network~\cite{ncf} above them for prediction (\emph{i.e., } Figure \ref{fig:rec_graph}(a)). Item popularity $Z$ is not explicitly considered in most models, however, it indeed affects the learning of item embedding. For example, it increases the vector length of popular items, making inner product models score popular items high for every user~\cite{DICE}. Such effect justifies the edge $Z\rightarrow I$, which exerts a negative impact to learn real user interest.

To sum up, regardless of which explanation for the causal graph, $Z\rightarrow I$ causes the bad effect and should be eliminated in formulating the predictive model.

\subsection{Deconfounded Training}\label{ssec:PD}
We now consider how to obtain a model that is immune to the impact of $Z\rightarrow I$. Intuitively, if we can intervene the exposure mechanism to make it randomly expose items to users, then the collected interaction data is free from the impact of $Z\rightarrow I$. Directly training traditional models on it will do. However, the feasibility and efficacy of such solutions are low: first, only the recommender builder can intervene the exposure mechanism, and anyone else (\emph{e.g., } academia researcher) has no permission to do it; second, even for recommender builder that can intervene the exposure mechanism, they can only use small random traffic, since random exposure hurts user experience much. Effective usage of small uniform data remains an open problem in recommendation research~\cite{2018causalE, distillation}.

Fortunately, the progress of causal science provides us a tool to achieve intervention without performing the interventional experiments~\cite{pearl2009causality}. The secret is the \textit{do}-calculus. In our context, performing $do(I)$ forces to remove the impact of $I$'s parent nodes, achieving our target. As such, we formulate the predictive model as $P(C| do(U, I) )$, rather than $P(C| U, I )$ estimated by traditional methods\footnote{Note that since there is no backdoor path between $U$ and $C$ in our causal graph, it holds that $P(C| U, do(I) ) = P(C| do(U), do(I) ) = P( C| do (U, I))$. Here we use $P(C| do(U, I) )$ as the intervention model for brevity. }.

Let the causal graph shown in Figure~\ref{fig:rec_graph}(b) be $G$, and the intervened causal graph shown in Figure~\ref{fig:rec_graph}(c) be $G'$. Then, performing \textit{do}-calculus on $G$ leads to:
\begin{equation}\label{eq:do-derivation}\small
	\begin{split}
	P(C|do(U,I)) & \overset{(1)}{=}  P_{G'} (C|U,I)\\
	& \overset{(2)}{=} \sum_z P_{G'}(C|U,I,z)P_{G'}(z|U,I)\\
	& \overset{(3)}{=} \sum_z P_{G'}(C|U,I,z)P_{G'}(z)\\
	& \overset{(4)}{=} \sum_z P(C|U,I,z) P(z),
	\end{split}
\end{equation}
where $P_{G'}(\cdot)$ denotes the probability function evaluated on $G'$. Below explains this derivation step by step:
\begin{itemize}[leftmargin=*]
    \item (1) is because of \textit{backdoor criterion}~\cite{pearl2009causality} as the only backdoor path $I\leftarrow Z \rightarrow C$ in $G$ has been blocked by $do(U, I)$;
    \item (2) is because of Bayes' theorem;
    \item (3) is because that $U$ and $I$ are independent with $Z$ in $G'$;
    \item (4) $P(C|U,I,Z) = P_{G'}(C|U,I,Z)$ is because that the causal mechanism $\{U,I,Z\}\rightarrow C$ is not changed when cutting off $Z\rightarrow I$, $P(Z) = P_{G'}(Z)$ since $Z$ has the same prior on the two graphs.
\end{itemize}

Next, we consider how to estimate $P(C|do(U, I))$ from data. Apparently, we need to first estimate $P(C|U,I,Z)$, and then estimate $\sum_z P(C|U,I,z) P(z)$. \vspace{+5pt}

\noindent \textbf{Step 1. Estimating $P(C|U,I,Z)$.} This conditional probability function evaluates that given a user-item pair $U=u, I=i$ and the item's present popularity as $Z=m_i^t$, how likely the user will consume the item. Let the parameters of the conditional probability function be $\Theta$, we can follow traditional recommendation training to learn $\Theta$, for example, optimizing the pairwise BPR objective function on historical data $\mathcal{D}$:
\begin{equation}\label{eq:bpr}\small
    \max_{\Theta} \sum_{(u, i, j) \in \mathcal{D}} \text{log} \sigma \left(
        P_\Theta\left(
            c=1 | u, i, m_{i}^{t}
        \right) - P_\Theta \left(
            c=1 | u, j, m_{j}^{t}
        \right)
    \right),
\end{equation}
where $j$ denotes the negative sample for $u$, $\sigma(\cdot)$ is the sigmoid function. The $L_2$ regularization is used but not shown for brevity.

Having established the learning framework, we next consider how to parameterize $P_\Theta(c=1 | u, i, m_{i}^{t})$. Of course, one can design it as any differentiable model, factorization machines, or neural networks. But here our major consideration is to decouple the user-item matching with item popularity. The benefits are twofold: 1) decoupling makes our framework extendable to any collaborative filtering model that concerns user-item matching; and 2) decoupling enables fast adjusting of popularity bias in the inference stage (see Section~\ref{sec:PDA}), since we do not need to re-evaluate the whole model. To this end, we design it as:
\begin{equation} \label{eq:model}\small
    P_\Theta(c=1 | u, i, m_{i}^{t}) = ELU'(f_\Theta(u, i)) \times ({m}_{i}^{t})^\gamma,
\end{equation}
where $f_\Theta(u, i)$ denotes any user-item matching model and we choose the simple \textit{Matrix Factorization} (MF) in this work; hyper-parameter $\gamma$ is to smooth the item popularity and can control the strength of conformity effect: setting $\gamma=0$ means no impact and larger value assigns larger impact. $ELU'(\cdot)$ is a variant of the Exponential Linear Unit \cite{elu} activation function that ensures the positivity of the matching score:
\begin{equation} \label{eq:alph-example}\small
	ELU'(x) = \begin{cases}
	e^{x}, &\text{if } x \leq 0 \\
	x+1, &\text{else}
		   \end{cases}
\end{equation}
This is to ensure the monotonicity of the probability function since $({m}_{i}^{t})^\gamma$ is always a positive number. Lastly, note that one needs to normalize $ P_\Theta(c=1 | u, i, m_{i}^{t})$ to make it a rigorous probability function, but we omit it since it is time-consuming and does not affect the ranking of items. \vspace{+5pt}

\noindent \textbf{Step 2. Estimating  $\sum_z P(C|U,I,z) P(z)$.} Now we move forward to estimate the interventional probability $P(C|do(U, I))$. Since the space of $Z$ is large, it is inappropriate to sum over its space for each prediction evaluation. Fortunately,  we can perform the following reduction to get rid of the sum:
\begin{equation}\label{eq:doUI}\small
   \begin{split}
    P(C|do(U,I)) &= \sum_{z} P(C|U,I,z)P(z) \notag\\
    & = \sum_z ELU'(f_{\Theta}(u,i)) \times z^{\gamma} P(z) \\
    & = ELU'(f_{\Theta}(u,i)) \sum_z z^{\gamma}P(z) \\
    & = ELU'(f_{\Theta}(u,i)) E(Z^{\gamma})
\end{split}
\end{equation}
where $E(Z^{\gamma})$ denotes the expectation of $Z^{\gamma}$. Note that the expectation of one variable is a constant. Since $P(C|do(U, I))$ is used to rank items for a user, the existence of $E(Z^{\gamma})$ does not change the ranking. We can thus use $ELU'(f_{\Theta}(u,i))$ to estimate $P(C|do(U, I)$. \vspace{+5pt}

\noindent To summarize, we fit the historical interaction data with $P_\Theta(c=1|u, i, m_i^t)$, and use the user-item matching component $ELU'(f_{\Theta}(u,i))$ for deconfounded ranking. We name this method as \textit{Popularity-bias Deconfounding} (PD).

\subsection{Adjusting Popularity Bias in Inference}\label{sec:PDA}
Owing to $P(C|do(U, I))$, we can eliminate the bad effect of popularity bias. We now seek for better usage of popularity bias, \emph{e.g., } promoting the items that have the potential to be popular.
Assume the target popularity bias is $\Tilde{z}$, and we want to endow the recommendation policy with this bias. To achieve the target, we just need to do the intervention of $Z=\Tilde{z}$ for model inference:
\begin{align}\small\label{eq:do-z}
    P(C|do(U=u,I=i),do(Z=\Tilde{z})) = P_\Theta(c=1|u,i, \Tilde{m}_i),
\end{align}
where $\Tilde{m}_i$ denotes the popularity value of $\Tilde{z}$.
This intervention probability directly equals the conditional probability since there is no backdoor path between $Z$ and $C$ in the causal graph.
Since the focus of this work is not on popularity prediction, we employ a simple time series forecasting method to set $\Tilde{m}_i$:
\begin{equation}\small \label{eq:pre_pop}
    \Tilde{m}_i = m_{i}^{T} + \alpha (m_{i}^{T} - m_{i}^{T-1}),
\end{equation}
where $m_{i}^{T}$ is the popularity value of the last stage, and $\alpha$ is the hyper-parameter to control the strength of popularity drift in predicting the future. We name this method as \textit{Popularity-bias Deconfounding and Adjusting} (PDA).

Figure \ref{fig:framework} illustrates the workflow of PDA where the training stage optimizes Equation~(\ref{eq:bpr}) and the inference stage can adjust the popularity bias as:
\begin{equation}\label{eq:mf}\small
     PDA_{ui} = ELU'(f_\Theta(u, i)) \times (\Tilde{m}_i)^{\Tilde{\gamma}}.
\end{equation}
$\Tilde{\gamma}$ denotes the popularity smoothing hyper-parameter used for model inference, which can be different from that used in training.
This is to consider that the strength of popularity bias can drift in the future.
In practice, setting $\Tilde{\gamma} = \gamma$ can achieve expected performance (\emph{cf. } Table~\ref{tab:PDA-result}).
Setting $\Tilde{\gamma}=0$  degrades the method to PD that uses only interest matching for recommendation. Considering the importance of properly stopping model training, PDA uses a model selection strategy based on the adjusted recommendation, which is slightly different from that of PD (see Algorithm \ref{alg:PD-PDA}).

\begin{figure}
    \centering
    \includegraphics[width=0.35\textwidth]{figures/framework_5.pdf}
    \vspace{-0.4cm}
    \caption{The workflow of our method. The blue arrows represent the training stage and the red arrows represent the inference stage.}
    \label{fig:framework}
    \vspace{-13pt}
\end{figure}

\subsection{Comparison with Correlation $P(C|U, I)$}
In the beginning, we argue that traditional methods estimate the correlation $P(C|U, I)$ and suffer from the bad effect of popularity bias amplification. Here we analyze the difference between $P(C|U, I)$ and $P(C|do(U, I))$ to offer more insights on this point.

We first transform $P(C|U, I)$ with the following steps:
\begin{equation}\small
    \begin{split}
        P(C|U,I) & \overset{(1)}{=} \sum_{z} P(C,z|U,I) \\
        & \overset{(2)}{=} \sum_{z} P(C|U,I,z)P(z|U,I) \\
        & \overset{(3)}{=} \sum_{z} P(C|U,I,z)P(z|I) \\
        & \overset{(4)}{\propto	} \sum_{z} P(C|U,I,z)P(I|z)P(z)
    \end{split}
\end{equation}
(1) is the definition of marginal distribution; (2) is because of the Bayes' theorem; (3) is because $U$ is independent to $Z$ according to the causal graph; (4) is because of the Bayes' theorem.

Comparing with $P(C|do(U,I)) = \sum_{z} P(C|U,I,z) P(z)$,
we can see that $P(C|U, I)$ has an additional term: $P(I|Z)$, which fundamentally changes the recommendation scoring. Suppose $I=i$ is a popular item, $P(C|U, I=i, Z=m_i^t)$ is a large value since $m_i^t$ is large, which means the term has a large contribution to the prediction score. $P(I=i | Z=m_i^t)$ is also a large value due to the popularity bias in the exposure mechanism. Multiplying the two terms further enlarges the score of $i$, which gives $i$ a higher score than it deserves. As a result, the popularity bias is improperly amplified. This comparison further justifies the rationality and necessity of $P(C|do(U,I))$ for learning user interests reliably.

\begin{algorithm}[t]
	\caption{PD/PDA}
	\LinesNumbered
	\label{alg:PD-PDA}
	\KwIn{dataset $\mathcal{D}=\{(u,i,m_{i}^{t})\}$; hyper-parameter $\gamma$; predicted popularity $\{\Tilde{m}_i\}$; \textit{mode}: PD or PDA }
	\While{stop condition is not reached}{
	Update model parameters $\Theta$ by optimizing Equation~\eqref{eq:bpr}\;
	\lIf{mode==PD}{validate model with $ELU^{'}(f_\Theta (u,i))$}
	\lElse{validate model with Equation~\eqref{eq:mf} (simplify $\Tilde{\gamma}=\gamma$)}
	}

    \lIf{mode==PD}{recommend items using $ELU^{'}(f_\Theta (u,i))$}
    \lElse{recommend items using Equation~\eqref{eq:mf} (simplify $\Tilde{\gamma}=\gamma$)}
\end{algorithm}

\section{Experiments}
In this section, we conduct experiments to answer two main questions:
\textbf{RQ1:} Does PD achieve the goal of removing the bad effect of popularity bias?
How is its performance compared with existing methods?

\textbf{RQ2:} Can PDA effectively inject desired popularity bias? To what extent leveraging popularity bias enhance the recommendation performance.

\begin{table*}[]
\caption{Recommendation performance after deconfounded training on the three datasets. ``RI'' refers to the relative
improvement of PD over the corresponding baseline. The best results are highlighted in bold and sub-optimal results are underlined.}
\vspace{-8pt}
\label{tab:PD-result}
\resizebox{0.8\textwidth}{!}{%
\begin{tabular}{cccccccccccc}
\toprule
\multirow{2}{*}{Dataset} & \multicolumn{1}{c|}{\multirow{2}{*}{Methods}} & \multicolumn{5}{c|}{Top 20} & \multicolumn{5}{c}{Top 50} \\
 & \multicolumn{1}{c|}{} & Recall & \multicolumn{1}{l}{Precision} & \multicolumn{1}{l}{HR} & \multicolumn{1}{l}{NDCG} & \multicolumn{1}{l|}{RI} & \multicolumn{1}{l}{Recall} & \multicolumn{1}{l}{Precision} & \multicolumn{1}{l}{HR} & \multicolumn{1}{l}{NDCG} & \multicolumn{1}{l}{RI} \\ \hline
\multicolumn{1}{c|}{\multirow{6}{*}{Kwai}} & MostPop & 0.0014 & 0.0019 & 0.0341 & 0.0030 & 632.4\% & 0.0040 & 0.0021 & 0.0802 & 0.0036 & 480.9\% \\
\multicolumn{1}{c|}{} & BPRMF & 0.0054 & {\ul 0.0057} & 0.0943 & 0.0067 & 146.3\% & 0.0125 & {\ul 0.0053} & 0.1866 & 0.0089 & 121.0\% \\
\multicolumn{1}{c|}{} & xQuad & 0.0054 & {\ul 0.0057} & 0.0948 & 0.0068 & 145.0\% & 0.0125 & {\ul 0.0053} & 0.1867 & 0.0090 & 120.3\% \\
\multicolumn{1}{c|}{} & BPR-PC & {\ul 0.0070} & 0.0056 & {\ul 0.0992} & {\ul 0.0072} & 125.0\% & {\ul 0.0137} & 0.0046 & 0.1813 & {\ul 0.0092} & 123.7\% \\
\multicolumn{1}{c|}{} & DICE & 0.0053 & 0.0056 & 0.0957 & 0.0067 & 147.8\% & 0.0130 & 0.0052 & {\ul 0.1872} & 0.0090 & 119.0\% \\
\multicolumn{1}{c|}{} & PD & \textbf{0.0143} & \textbf{0.0138} & \textbf{0.2018} & \textbf{0.0177} & \textbf{-} & \textbf{0.0293} & \textbf{0.0118} & \textbf{0.3397} & \textbf{0.0218} & - \\ \hline
\multicolumn{1}{c|}{\multirow{6}{*}{Douban}} & MostPop & 0.0218 & 0.0297 & 0.2373 & 0.0349 & 75.4\% & 0.0490 & 0.0256 & 0.3737 & 0.0406 & 55.9\% \\
\multicolumn{1}{c|}{} & BPRMF & 0.0274 & {\ul 0.0336} & 0.2888 & 0.0405 & 47.0\% & 0.0581 & {\ul 0.0291} & 0.4280 & {\ul 0.0475} & 34.3\% \\
\multicolumn{1}{c|}{} & xQuad & 0.0274 & {\ul 0.0336} & {\ul 0.2895} & 0.0391 & 48.3\% & 0.0581 & {\ul 0.0291} & {\ul 0.4281} & 0.0473 & 34.4\% \\
\multicolumn{1}{c|}{} & BPR-PC & {\ul 0.0282} & 0.0307 & 0.2863 & 0.0381 & 51.6\% & {\ul 0.0582} & 0.0271 & 0.4260 & 0.0457 & 38.0\% \\
\multicolumn{1}{c|}{} & DICE & 0.0273 & {\ul 0.0336} & 0.2845 & {\ul 0.0421} & 46.2\% & 0.0513 & 0.0273 & 0.4000 & 0.0460 & 44.5\% \\
\multicolumn{1}{c|}{} & PD & \textbf{0.0453} & \textbf{0.0454} & \textbf{0.3970} & \textbf{0.0607} & \textbf{-} & \textbf{0.0843} & \textbf{0.0362} & \textbf{0.5271} & \textbf{0.0686} & - \\ \hline
\multicolumn{1}{c|}{\multirow{6}{*}{Tencent}} & MostPop & 0.0145 & 0.0043 & 0.0684 & 0.0093 & 340.8\% & 0.0282 & 0.0035 & 0.1181 & 0.0135 & 345.8\% \\
\multicolumn{1}{c|}{} & BPRMF & 0.0553 & {\ul 0.0153} & 0.2005 & 0.0328 & 27.1\% & 0.1130 & {\ul 0.0129} & 0.3303 & 0.0497 & 25.3\% \\
\multicolumn{1}{c|}{} & xQuad & 0.0552 & {\ul 0.0153} & 0.2007 & 0.0326 & 27.3\% & 0.1130 & {\ul 0.0129} & 0.3302 & 0.0497 & 25.3\% \\
\multicolumn{1}{c|}{} & BPR-PC & {\ul 0.0556} & {\ul 0.0153} & {\ul 0.2018} & {\ul 0.0331} & 26.5\% & {\ul 0.1141} & 0.0128 & {\ul 0.3322} & {\ul 0.0500} & 24.9\% \\
\multicolumn{1}{c|}{} & DICE & 0.0516 & 0.0149 & 0.1948 & 0.0312 & 32.8\% & 0.1010 & 0.0132 & 0.3312 & 0.0486 & 29.0\% \\
\multicolumn{1}{c|}{} & PD & \textbf{0.0715} & \textbf{0.0195} & \textbf{0.2421} & \textbf{0.0429} & \textbf{-} & \textbf{0.1436} & \textbf{0.0165} & \textbf{0.3875} & \textbf{0.0641} & - \\ \hline
\end{tabular}%
}
\end{table*}

\subsection{Experimental Settings}
\subsubsection{Datasets}\label{exp:setting}

We conduct experiments on three datasets:

1) \textbf{Kwai}: This dataset was adopted in Kuaishou User Interest Modeling Challenge\footnote{https://www.kuaishou.com/activity/uimc}, which contains click/like/follow records between users and videos. It spans about two weeks. In this paper, we only utilize clicking data. Following previous work~\cite{ngcf}, we take 10-core filtering to make sure that each user and each item have at least 10 interacted records. After filtering, the experimented data has 7,658,510 interactions between 37,663 users and 128,879 items.

2) \textbf{Douban Movie}: This dataset is collected from a popular review website Douban\footnote{www.douban.com} in China by~\cite{douban}. It contains user ratings for movies and spans more than ten years. We take all rating records as positive samples and only utilize the data after 2010. The 10-core filtering is also adopted. After filtering, there are 7,174,218 interactions between 47,890 users and 26,047 items.

3) \textbf{Tencent}: This dataset is collected from Tencent short-video platform from Dec 27 2020 to Jan 6 2021. In this dataset, user interactions are likes, which are reflective of user satisfaction but far more sparse than clicks. Because of its high sparsity, we adopt the 5-core filtering setting, obtaining 1,816,046 interactions between 80,339 users and 27,070 items.

Following \citet{FollowCrowd}, we split the datasets chronologically. In particular, we split the datasets into 10 time stages according to the interaction time, and each stage has the same time interval. The last stage is left for validation and testing, in which $30\%$ of users are regarded as validation set and the other $70\%$ of users are regarded as testing set. Note that we do not consider new users and new items in validation and testing.

\subsubsection{Baselines}
We compare PD with the following baselines:

\textbf{- MostPop.} This method simply recommends the most popular items for all users without considering personalization.

\textbf{- BPRMF,} which optimizes the MF model with BPR loss~\cite{bpr}.

\textbf{- xQuAD~\cite{xquad2019}.} This is a Ranking Adjustment method that aims to increase the coverage of long-tail items in the recommendation list. We use the codes released by the authors, and utilize it to re-rank the result of BPRMF. Following its original setting, we tune the bias controlling hyper-parameter $\lambda$ in $[0,1]$ with step 0.1.

\textbf{- BPR-PC~\cite{BPR_PC}.} This is a state-of-the-art Ranking Adjustment method for controlling popularity bias. It has two choices: re-ranking and regularization. Here we implement the re-ranking method based on the BPRMF, because this method shows better performance in the original paper. Following the paper, we tune the bias controlling hyper-parameters $\alpha$ [0.1, 2.0] with step 0.2, and $\beta$ in [0, 1] with step 0.2.

\textbf{- DICE~\cite{DICE}.} This is a state-of-the-art method for learning causal embeddings to handle the popularity bias.
It designs a framework with causality-oriented data to disentangle user interest and conformity into two sets of embedding.
We use the codes provided by the authors, and we take the settings suggested for large datasets instead of the default settings since our datasets are relatively large.
Because DICE demonstrates superior performance over IPS-based methods~\cite{IPS_rec,rec_ips_cap_norm}, we do not include IPS-based methods as baselines.

To evaluate the effect of PDA on leveraging desired popularity bias, we compare it with the following popularity-aware methods:

\textbf{- MostRecent~\cite{MostRecent}}. This is an improved version of MostPop, which recommends the most popular items in the last stage rather than in the entire history. Note that PDA takes the last stage to forecast the popularity of next stage to make adjustment.

\textbf{- BPRMF(t)-pop~\cite{temp_pop2017}}. This method models the temporal popularity by training a set of time-specific parameters for each stage. We directly utilize the parameters corresponding to the last training stage for validation and testing.

\textbf{- BPRMF-A} and \textbf{DICE-A}. We enhance BPRMF and DICE to inject desired popularity bias during inference in the same manner as PDA by substituting $f_{\Theta}(u,i)$ with the prediction of BPRMF/DICE in Equation~\eqref{eq:mf}. We tune the smooth hyper-parameter $\Tilde{\gamma}$.

\subsubsection{Hyper-parameters and Metrics.}
For a fair comparison, all methods are optimized by BPR loss and tuned on the validation data. We optimize all models with the Adam~\cite{adam} optimizer with batch size as $2,048$ and default choice of learning rate (\emph{i.e., } 0.001) for all the experiments. We search $L_2$ regularization coefficients in the range of $\{0.01, 0.001\dots,10^{-6},0\}$ for all models.
We adopt the early stopping strategy that stops training if Recall@20 on the validation data does not increase for 100 epochs. For DICE, we use its provided stopping strategy and tune the threshold for stopping to make it comparable to our strategy.
For PD and PDA, we search $\gamma$ in $[0.02,0.25]$ with a step size of $0.02$. For BPRMF-A and DICE-A, we search $\Tilde{\gamma}$ from $0.02$ to $1$ with a step size of $0.02$ and stop searching if the evaluation results do not change for $3$ steps.

To measure the recommendation performance, we adopt four widely-used evaluation metrics: Hit Ratio (HR), Recall,  Precision, which consider whether the relevant items are retrieved within the top-$K$ positions, and NDCG that measures the relative orders among positive and negative items in the top-$K$ list. All metrics are computed by the all-ranking protocol --- all items that are not interacted by a user are the candidates. We report the results of $K=20$ and $K=50$ in our experiments.

\subsection{Deconfounding Performance (RQ1)}
In this section, we will first study the recommendation performance of our proposed algorithm PD. Then, we analyze its recommendation lists and showcase its rationality.
At last, we study the necessity of computing item popularity at different stages in training (\emph{i.e., } computing local popularity).

\subsubsection{Overall Performance} Table~\ref{tab:PD-result} shows the comparison of top-$K$ recommendation performance. From the table we can find:
\begin{itemize}[leftmargin=*]
\item The proposed PD achieves the best performance, and consistently outperforms all the baselines on all datasets.
This verifies the effectiveness of PD, which is attributed to the deconfounded training to remove the bad effect of popularity bias. In addition, the superior performance of PD reflects the rationality of our causal analysis and the potential of \textit{do}-calculus in mitigating the bias issue in recommendation.
\item PD achieves different degrees of improvement on the three datasets. In particular, PD outperforms all baselines by at least $119\%$, $34\%$, and $24\%$ on Kwai, Douban, and Tencent, respectively. We think the difference in the improvement is due to the different properties of the three datasets. For example, the level of popularity drift in the data differs substantially as shown in Figure~\ref{fig:drit_showing}: the popularity drifts of Kwai are much larger than that of the other two datasets.
Higher drifts imply that item popularity $Z$ has larger impact on the recommendation data.
That is why the proposed PD outperforms baseline methods the most on the Kwai dataset. Hence, we believe that formulating the recommender model as $P(C|do(U, I))$ has larger advantages over $P(C|U, I)$ in the scenarios exhibiting higher popularity drift

\item For BPR-PC and xQuad, which perform ranking adjustment, only BPR-PC has limited improvements over BPRMF on recall, \emph{i.e., } ranking adjustment has
little effect in our setting. The reason is that: heuristic adjustment of the original model prediction can only guarantee an increase in the proportion of long-tail items but cannot guarantee the promoted items are of real interest. BPR-PC performs better than xQuad since  xQuad only considers group-level popularity, while BPR-PC computes popularity in a finer granularity at item level

\item DICE tries to disentangle popularity and interest influence for an observed interaction to get causal embedding. It only demonstrates a similar performance to BPRMF in our setting, which indicates the rationality of leveraging popularity bias instead of blindly eliminating the bias.
Note that in the original paper of DICE, the evaluation data is manually intervened to ensure an even distribution of testing interactions over items, which removes popularity bias in the testing data. Whereas in our setting, the evaluation data is split by time to reflect real-world evaluation of recommender system, and contains popularity bias.
Although DICE boosts the scores of unpopular items over popular items, it may under-estimate users' interest matching for popular items.

\end{itemize}

\begin{table*}[]
\caption{Top-K recommendation performance with popularity adjusting on Kwai, Douban, and Tencent Datasets.}
\vspace{-8pt}
\label{tab:PDA-result}
\resizebox{0.9\textwidth}{!}{%
\begin{tabular}{c|c|cccc|cccc|cccc}
\hline
\multicolumn{2}{c|}{Datasets} & \multicolumn{4}{c|}{Kwai} & \multicolumn{4}{c|}{Douban} & \multicolumn{4}{c}{Tencent} \\ \hline
\multicolumn{2}{c|}{} & \multicolumn{2}{c|}{top 20} & \multicolumn{2}{c|}{top 50} & \multicolumn{2}{c|}{top 20} & \multicolumn{2}{c|}{top 50} & \multicolumn{2}{c|}{top 20} & \multicolumn{2}{c}{top 50} \\
\multicolumn{2}{c|}{\multirow{-2}{*}{Methods}} & Recall & \multicolumn{1}{c|}{NDCG} & Recall & NDCG & Recall & \multicolumn{1}{c|}{NDCG} & Recall & NDCG & Recall & \multicolumn{1}{c|}{NDCG} & Recall & NDCG \\ \hline
\multicolumn{2}{c|}{MostRecent} & 0.0074 & 0.0096 & {\color[HTML]{0D0D0D} 0.0139} & 0.011 & 0.0398 & 0.0582 & 0.0711 & 0.0615 & 0.0360 & 0.0222 & 0.0849 & 0.0359 \\ \cline{1-2}
\multicolumn{2}{c|}{BPRMF(t)-pop} & 0.0188 & 0.0241 & 0.0372 & 0.0286 & 0.0495 & 0.0682 & {\ul 0.0929} & 0.0760 & 0.1150 & 0.0726 & 0.2082 & 0.1001 \\ \cline{1-2}
 & (a) & 0.0191 & 0.0249 & 0.0372 & 0.0292 & 0.0482 & 0.0666 & 0.0898 & 0.0744 & 0.1021 & 0.0676 & 0.1805 & 0.0905 \\
\multirow{-2}{*}{BPRMF-A} & (b) & 0.0201 & 0.0265 & 0.0387 & 0.0306 & 0.0486 & 0.0667 & 0.0901 & 0.0746 & 0.1072 & 0.0719 & 0.1886 & 0.0953 \\ \cline{1-2}
 & (a) & 0.0242 & 0.0315 & 0.0454 & 0.0363 & 0.0494 & 0.0681 & 0.0890 & 0.0736 & 0.1227 & 0.0807 & 0.2161 & 0.1081 \\
\multirow{-2}{*}{DICE-A} & (b) & 0.0245 & 0.0323 & 0.0462 & 0.0370 & 0.0494 & 0.0680 & 0.0882 & 0.0734 & 0.1249 & 0.0839 & 0.2209 & 0.1116 \\ \cline{1-2}
 & (a) & {\ul 0.0279} & {\ul 0.0352} & {\ul 0.0531} & {\ul 0.0413} & {\ul 0.0564} & \textbf{0.0746} & \textbf{0.1066} & \textbf{0.0845} & {\ul 0.1357} & {\ul 0.0873} & {\ul 0.2378} & {\ul 0.1173} \\
\multirow{-2}{*}{PDA} & (b) & \textbf{0.0288} & \textbf{0.0364} & \textbf{0.0540} & \textbf{0.0429} & {\textbf{0.0565} } & {\ul 0.0745} & {\textbf{0.1066}} & {\ul 0.0843} & \textbf{0.1398} & \textbf{0.0912} & \textbf{0.2418} & \textbf{0.1210} \\ \hline
\end{tabular}
}
\vspace{-10pt}
\end{table*}

\begin{figure}
\centering
\subfigure[\vspace{-5pt}\textbf{ Kwai}]{\includegraphics[width=0.23\textwidth]{figures/group_rec_kwai_n_int_RR_type.pdf}}
\subfigure[ \textbf{Douban}]{ \includegraphics[width=0.235\textwidth]{figures/group_rec_douban_n_int_RR_type.pdf}}
\quad
\subfigure[ \textbf{Tencent}]{ \includegraphics[width=0.23\textwidth]{figures/group_rec_tencent_n_int_RR_type.pdf}}
\subfigure[ \textbf{std. dev.}]{ \includegraphics[width=0.23\textwidth]{figures/std_n_int_type.pdf}}
\vspace{-15pt}
\caption{Recommendation rate(RR) over item groups.}
\label{fig:group-rec}
\vspace{-18pt}
\end{figure}

\subsubsection{Recommendation Analysis}
We claim that PD can alleviate the amplification of popularity bias, \emph{i.e., } removing the bad effect of popularity bias. We conduct an analysis of the recommendation lists to show that deconfounded recommendation does not favor popular items as much as conventional recommendations.
Items are placed into $10$ different groups in two steps: (1) we sort items according to their popularity in descending order and (2) we divide the items into $10$ groups and ensure that the \emph{sum} of popularity over items in each group is the same. As such, the \emph{number} of items will increase from group 1 to group 10 and the \emph{average} popularity among items in each group decreases from group 1 to group 10. Hence we say group 1 is the most popular group and group 10 is the least popular group.

For each method, we sum the number of times of items recommended from each group. Then we divide this number by the total number of recommendations to get the recommendation rate for each group.  Formally, we define recommendation rate (RR) for group $g$ as:
\begin{equation}\nonumber \small
    RR(g,algo)= \frac{\sum_{i\in Group_g}{RC(i, algo)}}{\sum_{i\in All Groups}{RC(i,algo)}}
\end{equation}
where $RC(i, algo)$ gives the number of times the algorithm $algo$ recommends item $i$ in  top-$K$ recommendation lists.

For training, the recommendation rate of different groups is uniform (see the black line in Figure~\ref{fig:group-rec}). If more popular groups have higher recommendation rates, there is the amplification of popularity bias. We show the recommendation rates of each group by different methods in Figure~\ref{fig:group-rec}, in which (a), (b), and (c) are the results on Kwai, Douban, and Tencent, respectively. Figure~\ref{fig:group-rec} (d) is to show the standard deviation (std. dev.) of recommendation rates over all the groups: a smaller value means the recommendation rate is more uniform and is closer to the recommendation rates from training set. The main findings are as follows.
\begin{itemize}[leftmargin=*]
\item It is clear that popularity bias will be amplified by BPRMF. Figures~\ref{fig:group-rec}~(a) (b) (c) show that as the average item popularity decreases from group 1 to group 10, the recommendation rate decreases almost linearly. The gap between the chances to recommend popular items and unpopular items becomes bigger in the recommendation lists compared with that in the training set.

\item For DICE, the recommendation rate increases initially and then decreases from group 1 to group 10. This is because DICE strategically increases the scores of unpopular items to suppress the problem of amplified interest for popular items. However, this suppression brings two side effects: (1) the items from group 1 and group 2 are over-suppressed such that they have lower chances to be recommended than they have in the training set on Kwai and Tencent and (2) it over-amplifies the interest for the middle groups that contain sub-popular items. As such, blindly eliminating popularity bias would remove the beneficial patterns in the data.

\item In Figure~\ref{fig:group-rec}~(a) (b) (c), the RRs of PD exhibit lines that are flatter than other methods, which are the most similar to the training set. Moreover, as we can see in Figure~\ref{fig:group-rec}~(d), the RR of PD has the smallest standard deviation across groups, indicating the least difference across groups. These results verify that PD provides relatively fair recommendations across groups, and will not over-amplify or over-suppress interest for popular items. Along with the recommendation performance in Table~\ref{tab:PD-result}, the results demonstrate the rationality and effectiveness of only removing the bad effect of popularity bias by deconfounded training

\item For the group $10$, which has the least popular items, PD has a little higher recommendation rate than other methods. However, PD also cannot provide enough chances of recommendation because of the sparse interactions of items in this group and the
insufficient learning of item embedding.
\end{itemize}

Note that we omit the RRs of the ranking adjustment methods (\emph{i.e., } BPR-PC and xQuad) since their recommendation rates for unpopular items are manually increased. Instead of tuning the recommendation rate for groups of different popularity, we concentrate more on whether the learned embedding can precisely reveal the real interests of users.

\subsubsection{Global Popularity \textit{V.S.} Local Popularity.} To investigate the necessity of local popularity, we evaluate a variant of PD that uses global popularity defined on the interaction frequency in the entire training set $\mathcal{D}$. We keep the other settings the same as PD and name this method as PD-Global (PD-G). Figure~\ref{fig:PD-G} shows the recommendation performance of PD and PD-G on the three datasets. Due to the limited space, we only show the results on recall and NDCG for top-$50$ recommendation, other metrics show the similar results. According to the figure, PD-G is far worse than PD, which verifies the importance of calculating item popularity locally. We believe PD-G performs poorly because there is only one instance of $P(C|u, i, Z)$ on $Z$, making $P(C|u, i, Z)$ hard to be learned well.

\begin{figure}
    \centering
    \subfigure[\textbf{ Recall}]{\label{fig:a}\includegraphics[width=0.23\textwidth]{figures/PD_PDG_Recall_50_type.pdf}}
    \subfigure[ \textbf{NDCG}]{\label{fig:b} \includegraphics[width=0.23\textwidth]{figures/PD_PDG_NDCG_50_type.pdf}}
    \vspace{-15pt}
    \caption{Comparisons between PD and PD-G on Recall@50 and NDCG@50. PD-G is a version of PD that computes item popularity over the total training set.}
    \vspace{-15pt}
    \label{fig:PD-G}
\end{figure}

\subsection{Performance of Adjusting Popularity (RQ2)}
In this section, we verify the performance of our proposed method PDA, which introduces some desired popularity bias into recommendation. Firstly, we compare our method with the baselines that inject some popularity bias as well --- MostRecent, BPRMF (t)-pop, BPRMF-A, and DICE-A. Secondly, we study the influence of refined prediction of item popularity for PDA. Lastly, we demonstrate overall improvements compared with the baseline model BPRMF.

\subsubsection{Comparisons with Baselines.} We utilize two methods to predict the popularity of the testing stage. Specially, method (a) takes the popularity of items from the last stage $\mathcal{D}_T$ as the predicted popularity, while method (b) computes populairity as Equation~(\ref{eq:pre_pop}).
We apply both (a) and (b) to PDA, DICE-A, and BPRMF-A. The results are shown in Table~\ref{tab:PDA-result}, which lead to the following observations:
\begin{itemize}[leftmargin=*]
\item By comparing the results of PDA, DICE-A, BPRMF-A, and MostRecent with the results of PD, DICE, BPRMF, and MostPop in Table~\ref{tab:PD-result}, we find that introducing desired popularity bias into recommendations can improve the recommendation performance. This validates the rationality of leveraging popularity bias, \emph{i.e., }  promoting the items that have the potential to be popular (desired popularity bias) benefits recommendation performance.

\item Compared with other methods, PDA can outperform both the method BPRMF(t)-pop, which is specially designed for utilizing temporal popularity, and the modified methods BPRMF-A and DICE-A. This is attributed to PDA introducing the desired popularity bias by an intervention, preventing the learning model to amplify the effect of popularity and injecting desired popularity bias in inference, whereas other models are influenced by popularity bias in training.
 The results validate the advantage of causal intervention in more accurately estimating true interests.
\item On Kwai and Tencent, method (b) outperforms method (a), indicating that introducing the linear predicted popularity is more effective than straightforwardly utilizing the popularity from the previous stage; On Douban, these two methods achieve almost the same performance. This is because the Douban dataset spans over a longer period of time than the other two datasets, making it harder to predict the popularity from the trends of popularity in training stages.
\end{itemize}

\begin{figure}
    \centering
    \includegraphics[width=0.49\textwidth]{figures/PD_sub_2_type.pdf}
    \vspace{-15pt}
    \caption{PDA recommendation performance regarding Recall@50 and NDCG@50 on Kwai and Douban, with the last training stage  split into $N$ sub-stages. }
    \vspace{-15pt}
    \label{fig:PD_sub}
\end{figure}

\subsubsection{More Refined Predictions for Popularity.} We conduct an experiment to show that the performance of PDA can be improved with more refined predicted popularity.
To get more refined predictions, we uniformly split the data in the last training stage into $N$ sub-stages by time, and achieve a more refined prediction of popularity in the last sub-stage, which is closest to the testing stage. We set $N=1,2,4,8,16$ and employ method (a) to predict the desired popularity. The results are shown in Figure~\ref{fig:PD_sub}. The performance of PDA can be improved with better-predicted popularity for testing. On the Kwai dataset, the performance first increases and then slightly decreases as $N$ increases. However, on the Douban dataset, with a larger $N$, the performance does not drop evidently.
This is mainly because the item space of Kwai is far larger than Douban, but they have a similar number of interactions. When $N$ is set to a large value, the interaction frequency is insufficient to accurately compute item popularity for Kwai. This experiment shows that the performance of PDA can be further improved with a more precise prediction of popularity.

\subsubsection{Total Improvements of PDA} To quantify the strengths of PDA, we show the relative improvements of PDA over the basic model that PDA is built on, \emph{i.e., } BPRMF in Table~\ref{tab:total-improvement}. On the Kwai dataset, we get more than $241\%$ improvements over all metrics. On the Douban and Tencent datasets, for the recall and NDCG, we get at least $97\%$ improvements, and for precision and HR, we get at least $57\%$ improvements. The results clearly indicate the superiority of PDA in recommendation and the huge potential of causal analysis and inference for leveraging bias in recommendation.

\section{Related Works}
This work is closely related to two topics: popularity bias in recommender system and causal recommendation.
\begin{table}[]
\caption{The total relative improvements of PDA compared with the base recommendation model, \emph{i.e., } BPRMF. Each result is averaged on top 20 and top 50.}
\vspace{-10pt}
\label{tab:total-improvement}
\resizebox{0.4\textwidth}{!}{%
\begin{tabular}{c|cccc}
\toprule
Dataset & Recall & Precision & HR & NDCG \\ \hline
Kwai & 488\% & 441\% & 241\% & 532\% \\
Douban & 124\% & 68\% & 57\% & 97\% \\
Tencent & 133\% & 133\% & 79\% & 161\% \\ \hline
\end{tabular}
}
\vspace{-15pt}
\end{table}

\subsection{Popularity Bias}
Popularity bias has long been recognized as an important issue in recommendation systems.
One line of research~\cite{connection2020,UMAP_recbias,FollowCrowd,FeedbackLoop,PopAcc} studies the phenomenon of popularity bias.
\cite{UMAP_recbias} studies how popularity bias influences  different recommendation models and concludes that popularity has a larger impact on BPR-based models \cite{DICE} than point-wise models \cite{ncf}; \cite{connection2020} tries to study the relation between popularity bias and other problems such as fairness and calibration.
Another line of research~\cite{DICE,BPR_PC,xquad2019,Himan_reg,Himan-thesis,IPS_rec,2018causalE} attempts to improve  recommendation performance by removing popularity bias. These works mainly consider three types of methods: IPS-based methods~\cite{IPS_rec,rec_ips_cap_norm}, causal embedding-based methods~\cite{2018causalE,distillation,DICE}, and the methods based on post-hoc re-ranking or model regularization~\cite{BPR_PC,xquad2019,Himan_reg,correcting-pop-bias,2011-PPT}.
Similar to our work, a very recent work \cite{DICE}  also tries to compute the influence of popularity via causal approaches.
The difference is that we analyze the causal relations from a confounding perspective whearas this important aspect was not considered in their work.
Some other works have attempted to remove the impact of other bias problems~\cite{chen_survey},  position bias~\cite{position-bias}, selection bias~\cite{wang2019doubly}, and exposure bias~\cite{exposure-modeling,chen2018modeling} \textit{etc}. The most widely considered method for handling aforementioned biases is IPS-based method~\cite{IPS_rec,ips_cjl,unbiased_PS,postclick-ips}.
However, there is no  guaranteed accuracy for estimating propensity scores. Utilizing unbiased data is another option for solving bias problems~\cite{2018causalE,distillation},   however the unbiased data is difficult to obtain and may hurt user experience~\cite{wang2019doubly}.

Some other  methods recognize the importance of utilizing popularity and consider    drifted popularity to improve the model performance~\cite{temp_pop2017,MostRecent,POP_music,ECIR-localPop}.
\cite{POP_music} tries to utilize temporal popularity  in music recommendation and \cite{POP-news}  utilizes local popularity for recommending  news. \cite{temp_pop2017,ECIR-localPop,MostRecent}  also try to utilize temporal popularity for general recommendation. The biggest difference between our work and other works is that we precisely differentiate popularity bias to be removed and popularity drift to be taken advantage of.

\vspace{-8pt}
\subsection{Causal Recommendation}
Some efforts have considered causality in recommendation. The first type of work is on confounding effects~\cite{CauInferRec,CE_for_rec,alg-confounding,deconfounding-response-rate}.  For example, \cite{CauInferRec}  takes the  de-confounding technique in linear models to learn real interest influenced by unobserved confounders. They learn substitutes for unobserved confounders by fitting exposure data.
\cite{CE_for_rec} tries  to estimate the causal effect of the  recommending operation. They consider the features of users and items as confounders, and reweight training samples to solve the confounding problem.
\cite{alg-confounding} explores the impact of algorithmic confounding on simulated recommender systems. \cite{deconfounding-response-rate} identifies  expose rate as a confounder for user satisfaction estimation and uses IPS to handle the confounding problem.
In this work, we introduce a different type of confounder that is brought by item popularity.

The second type is counterfactual learning methods~\cite{counterfactual-music,double_unbias_data,CVIB}. \cite{counterfactual-music} utilizes a quasi-experimental Bayesian framework to quantify the effect of treatments on an outcome to generate counterfactual data, where counterfactual data refers to the outcomes of unobserved treatments.
\cite{double_unbias_data} utilizes a model that is learned with few unbiased data to generate labels for unobserved data.
\cite{CVIB} extends  Information Bottleneck~\cite{CVIB} to counterfactual learning with the idea that the counterfactual world should be as informative as the factual world to deal with the  Missing-Not-At-Random problem~\cite{wang2019doubly}.
Though these works study the problem of causal recommendation, they do not specifically consider popularity bias as well as popularity drift for boosting recommendation performance.

\vspace{-5pt}
\section{conclusion}
In this work, we study how to model and leverage popularity bias in recommender system. By analyzing the recommendation generation process with the causal graph, we find that item popularity is a confounder that misleads the estimation of $P(C|U,I)$ as done by traditional recommendation models.
Different from existing work that eliminates popularity bias, we point out that some popularity bias should be decoupled for better usage.
We propose a Popularity-bias Deconfounding and Adjusting  framework to deconfound and leverage popularity bias for better recommendation.
We conduct experiments on three real-world datasets, providing insightful analyses on the rationality of PDA.

This work showcased the limitation of pure data-driven models for recommendation, despite their dominant role in recommendation research and industry. We express the prior knowledge of how the data is generated for recommendation training as a causal graph, performing deconfounded training and intervened inference based on the causal graph.
We believe this new paradigm is general for addressing the bias issues in recommendation and search systems, and will extend its impact to jointly consider more factors, such as position bias and selection bias.
Moreover, we are interested in extending our method to the state-of-the-art graph-based recommendation models and incorporate the content features of users and items.


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

\subsection{Related Work}
\label{sec:relatedwork}
\input{sec-relatedwork}
\section{Methodology}
\label{sec:method}
\input{sec-method}
\section{Experimental Results}
\label{sec:experiment}
\input{sec-experiment}
\section{Conclusion}
This paper presents a novel topological biomarker, \emph{TopoTxR}, that leverages the rich geometric information embedded in structural MRI and enables improvement in downstream CNN processing. In particular, we compute 1D cycles and 2D bubbles from breast DCE-MRIs with the theory of persistent homology; these structures are then used to direct the attention of neural networks. We further demonstrate that \emph{TopoTxR} on treatment-naive imaging is predictive of pCR.

\bibliographystyle{abbrv}

\subsection{Discussion and TopoTxR Feature Interpretation}
The topological structures extracted by TopoTxR capture the breast tissue structures.
Learning directly on these structures and their vicinity provides the opportunity for interpreting the learning outcomes and drawing novel biological insights.
Here we provide some visual analysis as a proof of concept.

Fig.~\ref{fig:qualitative} shows the TopoTxR topographical structures from different views for a representative DCE-MRI scan from each group. We observe that the structures (1D and 2D) are sparse for the case exhibiting pCR, and are relatively dense for the non-pCR case. In the corresponding MRI images, we note that the pCR breast has scattered fibroglandular breast density with minimal background parenchymal enhancement. The non-pCR breast has a more heterogenous fibroglandular breast density with moderate background parenchymal enhancement. This possibly suggests that the \textit{TopoTxR} features capture the complex fibrogladular structure which can be a potential indicator of treatment response.

We also compare the topological behavior of the two populations. Recall the birth time of a topological structure is the threshold at which a cycle appears. In our experiments, since we use the inverse image $f=-I$, the birth time essentially captures -1 times the brightness of a structure. In Fig.~\ref{fig:qualitative} (right), we plot the cumulative density function (CDF) of the birth time of topological structures for pCR (red) and non-pCR (blue) patients. The CDFs suggest that pCR patients' tissue structures are generally less bright (or less visible) compared with that of non-pCR patients. This is consistent with our observation on qualitative examples.
A Kolmogorov-Smirnov test \cite{massey1951kolmogorov} is performed to compare these CDFs. The computed \textit{p-value} is 0.0002, indicating a significant difference between the distributions of birth times of the pCR and non-pCR patient groups.

\subsection{Background: Persistent Homology}
We review basics of persistent homology in this section. Interested readers may refer to \cite{edelsbrunner2010computational} for more details.
Persistent homology extracts the topological information of data observed via a scalar function. Given an image domain, $X$, and a real-valued function $f:X\rightarrow \mathbb{R}$, we can construct a sublevel set $X_{t} =
\{x \in X: f(x) \leq t\}$ where $t$ is a threshold controlling the ``progress"
of sublevel sets. The family of sublevel sets $\mathcal{X} = \{X_{t}\}_{t \in \mathbb{R}}$ defines a filtration, i.e., a family of subsets of $X$ nested with respect to the inclusion: $X_{\alpha} \subseteq X_{\beta}$ if $\alpha \leq \beta$. As the threshold $t$ increases from $-\infty$ to $+\infty$, topological structures such as connected components, handles, and voids appear and disappear. The birth time of a topological structure is the threshold $t$ at which the structure appears in the filtration. Similarly, the death time is the threshold at which the structure disappears. Persistent homology tracks the topological changes of sublevel sets $X_{t}$ and encodes them in a \emph{persistence diagram}, i.e., a point set in which each point $(b,d)$ represents a topological structure born at $b$ and killed at $d$.

\begin{figure}[btph]
\begin{center}
\includegraphics[width=.96\textwidth]{figs/ipmi21-phf.pdf}
\caption[ph]{From left to right: a synthetic image $f$, sublevel sets at thresholds $b_{1} < b_{2} < d_{2} < d_{1}$ (in black), and the persistence diagram of dimension 1. The red loop represents a 1D structure born at $b_1$ and killed at $d_1$. The green loop represents a 1D structure born at $b_{2}$ and killed at $d_2$. They correspond to the red and green dots in the diagram. }
\label{fig:ph}
\end{center}
\end{figure}

See Fig.~\ref{fig:ph} for an example function $f$ and its sublevel sets at different thresholds. At time $b_1$, a new handle (delineated by the red cycle $c_{1}$) is created. This handle is later destroyed at time $d_{1}$. Another handle delineated by the green cycle $c_{2}$ is created and killed at $b_{2}$ and $d_{2}$ respectively. The topological changes are summarized in a persistence diagram on the right. Each handle corresponds to a 2D dot in $\mathbb{R}^{2}$, whose $x$ and $y$ coordinates are birth and death times. The difference between a dot's birth and death times is called its \emph{persistence}.

\subsection{Persistence Cycles and Their Computation}
In this section, we introduce cycles that represent topological structures discovered by persistent  homology. We also present an algorithm to compute these cycles.
Intuitively, a topological cycle of dimension $p$ is a $p$-manifold without boundary. A 1-dimensional (1D) cycle is a loop (or a union of a set of loops). A 2-dimensional (2D) cycle is a bubble (or a union of a set of bubbles).
A cycle $z$ represents a persistent homology structure if it delineates the structure at its birth. For example, in Fig.~\ref{fig:ph}, the red and the green loops represent the handles born at time $b_1$ and $b_2$, respectively.

We assume a discretization of the image domain into distinct elements, i.e., vertices (corresponding to voxels), edges connecting adjacent vertices, squares, and cubes.
These elements are 0-, 1-, 2-, and 3-dimensional cells.
Any set of $p$-cells is called a \emph{$p$-chain}.
We define the \emph{boundary operator} of a $p$-cell, $\sigma$, as the set of its $(p-1)$ faces.
The boundary of an edge is its two vertices. The boundary of a square consists of the 4 edges enclosing it. The boundary of a cube consists of the 6 squares enclosing it.
For any $p$-chain, $c$, its boundary is
$\partial(c) = \sum_{\sigma\in c} \partial(\sigma)$, under mod-2 sum.
For a set of edges forming a path, its boundary are the two end vertices. For any set of squares forming a patch, its boundary is the loop enclosing the patch. The boundary of a set of cubes is the bubble enclosing them.

\begin{figure}[btph]
\begin{center}
\includegraphics[width=.88\textwidth]{figs/ipmi21-cycle.pdf}
\caption[cycle]{(a) Example of a cubical complex whose cells are sorted according to the function values. (b) 2D Boundary matrix $\partial$. (c) Reduced boundary matrix. (d) Persistence diagram and resulting cycles corresponding to $\partial$. (e) 1D boundary matrix.}
\label{fig:cycle}
\end{center}
\end{figure}

A $p$-chain is a \emph{$p$-cycle} if its boundary is empty.
All $p$-cycles form the null space of the boundary operator, i.e. $\{c : \partial(c)=\emptyset\}$. Any topological structure, formally defined as a \emph{homology class}, can be delineated by different cycles, which are equivalent to each other in terms of topology. We can choose any of them to represent this class. In persistent homology, for each dot in the persistence diagram, we can represent it with one representative cycle at its birth. In Fig.~\ref{fig:ph}, the red and green cycles represent the two corresponding handles.
Note that the choice of representative cycle is not unique. A relevant question is to choose the shortest representative cycle (i.e., one with the least number of edges) for each dot in the diagram \cite{wu2017optimal,zhang2019heuristic}. In this paper, we focus on choosing a standard representative cycle, leaving the optimal cycle for future work.

\myparagraph{Computation of persistent homology and representative cycles.}
We assume a filtration function on a discretization of the image domain.
An example discretization of a 2D image is given in Fig.~\ref{fig:cycle}(a). We first sort all cells in increasing order according to their function values. The computation of persistence diagrams is then performed by encoding the $p$-dimensional boundary operator in binary matrices named \emph{boundary matrices}, $\partial_p$. $\partial_p$ maps $p$-cells to their boundaries. Fig.~\ref{fig:cycle} shows the 1D and 2D boundary matrices of the given complex and its filtration. The 1D boundary matrix is essentially the incidence matrix of the underlying graph (Fig.~\ref{fig:cycle}(e)).
High dimensional boundary matrices are defined similarly (e.g., a 2D boundary matrix in Fig.~\ref{fig:cycle}(b)).

Persistence diagram is computed by reducing the boundary matrix similar to Gaussian elimination, but without row or column perturbation. We reduce through column operations performed on $\partial$ from left to right.
Fig.~\ref{fig:cycle}(c) shows the reduced 2D boundary matrix.
Once the boundary matrices are reduced. Each non-zero column corresponds to a persistent dot in the diagram. The reduced column itself is the cycle representing the corresponding topological structure. In this paper, we pay attention to both 1D and 2D cycles, corresponding to loops and bubbles. The extracted cycles will be used to explicitly guide 3D CNNs for analysis.
The computation of topological cycles is of the same complexity as the computation of persistent homology. In theory, it takes $O(n^\omega)$ time ($\omega \approx 2.37$ is the exponent in the matrix multiplication time, i.e., time to multiply two $n\times n$ matrices) \cite{milosavljevic2011zigzag}. Here $n$ is the number of voxels in an image. In practice, to compute all cycles of an input image ($256^3$), it takes approximately 5 minutes.
\subsection{Topological-Cycle-Driven 3D CNN}
An overview of our topological-cycle-driven 3D CNN has been provided in the beginning of Section \ref{sec:method}.
Here we describe the technical details.

To compute persistence and topological cycles, we invert the MRI image $f = -I$ so that the tissue structures correspond to low intensity.
After the computation, we select topological cycles representing dots of the diagram with high persistence.
The general belief is that low-persistence dots tend to be due to noise. Thus we only select high-persistence cycles, which are considered more salient structures, and more likely to represent true tissue structures. The threshold is a hyperparamter tuned in practice.

Next, we create two binary 3D masks representing 1D and 2D topological cycles, respectively. Both masks are dilated slightly in order to cover both the structures and their vicinity.
Instead of directly using these binary masks for pCR prediction, we fill the foreground voxels with their original image intensity values. In other words, we mask the input image with the complement of the cycle mask. See\emph{Topological Structure Masking} step in Fig.~\ref{fig:overview} for the masked MRI image. We generate masked images for both 1D and 2D, and provide them to two CNNs.
All masked MRIs are padded to the same size of 256 $\times$ 256 $\times$ 256.

We use separate networks with the same architecture for cycles and bubbles. The CNN consists of 5 3D convolution layers, each followed by a batch normalization layer and a LeakyReLU. The output feature maps from these two 3D CNNs are reshaped and concatenated into a feature vector. This feature vector is sent into a fully connected (FC) network with three FC layers for final pCR prediction. Besides ReLU and batch normalization, a dropout layer is added to the second FC layer. The final output is a vector of size 2. All three networks are trained together in an end-to-end fashion with stochastic gradient descent (SGD) as the optimizer and cross-entropy as loss.


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

Understanding the origin and evolution of the solar system is one long-standing problem.
Meteoritic data provide unique constraints on the physical and chemical environments of the solar nebula and on the formation of terrestrial planets and small bodies, such as asteroids.

Chondrules are one of the main components of chondrites and contain key information about solar system formation \citep[e.g.,][]{Scott2007}.
They are roughly millimeter-sized particles,
and some of them are coated with dust grains, which are called rims.
There are two types of rims: fine-grained rims \citep[FGRs, e.g.,][]{King_TVV&King_EA1981} and coarse-grained igneous rims \citep[e.g.,][]{Rubin1984,Matsuda+2019}.
One plausible explanation for the formation of FGRs is the accretion of tiny dust grains by chondrules in the solar nebula \citep[e.g., ][]{Morfill+1998,Liffman2019}.
After FGRs formed, some of them may be transformed into igneous rims by remelting events \citep[][]{Rubin2010}.
Rims and their structure recorded the history of dust accretion and remelting events that the host chondrules experienced in the first few million years of the solar system's history.

Studies of FGRs have revealed that the thickness of rims ($a_{\rm rim}$) increases linearly as the radius of chondrules ($a_{\rm ch}$) increases:
\begin{eqnarray}
	\label{eq:a_rim_obs}
	a_{\rm rim} \approx K a_{\rm ch} +c,
\end{eqnarray}
where $K$ and $c$ range from 0.02 to 0.42 and from 8~$\mu$m to 29~$\mu$m, respectively \citep[e.g.,][{see also Table \ref{table:data_Kc}}]{Metzler+1992,Sears+1993, Friend+2018, Hanna&Ketcham2018}
\footnote{
The linear relation is confirmed in carbonaceous chondrites, especially CM chondrites.
\cite{Bigolski2017_PhD} showed that the correlation is weak in the Semarkona chondrules and less evident in those found in the unequilibrated ordinary chondrites.
}.

Theoretical estimations of the FGR thickness reproduce the aforementioned relationship between $a_{\rm rim}$ and $a_{\rm ch}$.
\cite{Morfill+1998} considered that accretion of FGRs in a solar nebula is characterized by the local mass densities of chondrule and dust, which are given by $\rho_{\rm ch}$ and $\rho_{\rm d}$, respectively, and derived
\begin{eqnarray}
	\label{eq:a_rim_Morfill}
	a_{\rm rim} = a_{\rm ch} \left[ \left( \frac{\rho_{\rm ch}+\rho_{\rm d} }{ \rho_{\rm ch} + \rho_{\rm d} \exp{ (-t/\tau_{\rm acc} )}} \right)^{1/3} -1\right].
\end{eqnarray}
In the above equation, $\tau_{\rm acc}$ is the typical accretion timescale and can be written as
\begin{eqnarray}
		\tau_{\rm acc} = \frac{\rho_{\rm d}}{\rho_{\rm ch} + \rho_{\rm d} } \frac{m_{\rm ch} }{ Q \rho_{\rm d} \sigma \Delta v },
\end{eqnarray}
where $m_{\rm ch}$ is a chondrule mass, $Q$ is the sticking efficiency, $\sigma$ is the collisional cross-section, and $\Delta v$ is the relative velocity between dust and chondrules.
Equation (\ref{eq:a_rim_Morfill}) is effectively the reflection of mass conservation of the total dust mass around chondrules in the nebula.
Equation (\ref{eq:a_rim_Morfill}) indicates that when $t=0$, no dust grains are accreted, and hence $a_{\rm rim}=0$.
After $t\gg \tau_{\rm acc}$, almost all dust grains are accreted, and the rim thickness becomes
\begin{equation}
\label{eq:chi_local}
a_{\rm rim}=a_{\rm ch}( \chi_{\rm local}^{-1/3} -1 ),
\end{equation}
where $\chi_{\rm local}=\rho_{\rm ch}/(\rho_{\rm ch}+\rho_{\rm d})$ is the initial value of the local chondrule mass fraction.
The $K$ values of $0.1-0.42$ yield the values of $\chi_{\rm local}$ ranging from 0.35 to 0.75.

The above estimation has been improved by several studies, where the more realistic treatments of turbulence, the sticking efficiency, the internal structure of rims, and/or erosive collisions were adopted \citep[][]{Cuzzi+2004,Carballido2011, Xiang_C+2019,Xiang_C+2019b,Liffman2019}.
However, such improvements still do not allow for accurate estimation of the rim thickness.
This is because, as shown by the previous studies \citep{Ormel+2008, Arakawa2017,Matsumoto+2019}, the growth timescale of dust grains is shorter than the accretion timescale of dust grains onto chondrules.
Moreover, when the internal density evolution of growing dust aggregates is properly taken into account, dust growth is further accelerated; dust aggregates initially grow via hit-and-stick collisions, which decreases their bulk densities \citep[e.g.,][]{Blum2004, Wada+2008} and increases their collisional cross-sections.
In such a picture of dust growth, the rim accretion onto chondrules is expected to be more efficient.
Also, the rim thickness should become a function of the chondrule size, chondrule fraction, and disk parameters (c.f., Equation (\ref{eq:a_rim_Morfill})).

In this paper, we adopt a porous dust growth model and estimate the thickness of rims around chondrules.
We compare our results to the thicknesses of both FGRs and igneous rims to consider whether these rim sizes can be explained by accretion.
Our estimation and results are summarized in Section \ref{sec:rim}.
In Section \ref{sec:dis}, we present the conclusion of this work and discuss the formation processes of FGRs and igneous rims.

\section{Estimation of rim thickness} \label{sec:rim}

\subsection{Overview}\label{sect:model}

We consider a local region in a gas disk where chondrules and dust grains co-exist.
The growth of these chondrules and dust grains may be modeled as follows \citep[][]{Arakawa2017,Matsumoto+2019}:
(i) dust grains grow via their collisions, forming fluffy dust aggregates;
(ii) these dust aggregates then accrete onto chondrules;
and (iii) after chondrules accrete enough amount of these dust aggregates, they begin to collide with each other.
Hereafter, chondrules that accreted dust aggregates are referred to as compound aggregates.
The accreted dust particles are regarded as the source material of rims.
Hence, the rim thickness around chondrules can be estimated by specifying the total mass of accreted dust aggregates at a time when the collision between compound aggregates begins.
We hereafter refer this time as $\tau_{\rm rim}$ in the following discussion\footnote{It should be noted that $\tau_{\rm rim}$ is used only for descriptive purposes;
in our model, dust growth is computed, based on the timescale argument, and hence we do not follow the actual time evolution.
Moreover, our $\tau_{\rm rim}$ and $t$ in Equation (\ref{eq:a_rim_Morfill}) have different meanings since the rim accretion pictures are different.}.

We formulate here the chondrule mass fraction ($\chi$) of a compound aggregate as
\begin{equation}
	\chi = \frac{3 m_{\rm ch}}{
	4\pi \left[ \rho_{\rm rim} \left(a_{\rm rim}+ a_{\rm ch} \right)^3 -\rho_{\rm rim} a_{\rm ch}^3 + \rho_{\rm ch} a_{\rm ch}^3 \right]},
\end{equation}
where $m_{\rm ch}$, $\rho_{\rm ch}$, and $a_{\rm ch}$ are the mass, internal density, and radius of a chondrule, respectively, and $\rho_{\rm rim}$ is the density of the rim.
Under the assumption that dust grains around a chondrule are isotropically distributed \citep[][]{Bland+2011} and the densities of chondrules and rims are $\rho_{\rm ch}/\rho_{\rm rim}\sim1$ in chondrites, the above equation is re-written as
\begin{eqnarray}
	a_{\rm rim}
	&=& a_{\rm ch} \left( \left[ \left(\chi^{-1} -1\right) \frac{\rho_{\rm ch}}{\rho_{\rm rim}}+1 \right]^{1/3} -1 \right)
	\nonumber \\
	&\simeq& a_{\rm ch} (\chi^{-1/3} -1 ).
	\label{eq:estimation_ach}
\end{eqnarray}
One should note that while rims are porous in the accretion stage,
the assumption that $\rho_{\rm ch}/\rho_{\rm rim}\sim1$ is reasonable in estimating the rim thickness.
This is because fluffy dust rims are compacted in the subsequent evolution processes: the dynamical compression via high-velocity collisions \citep[e.g.,][]{Wada+2008}; static compression via the ram pressure of the disk gas and via self-gravity \citep[e.g.,][]{Seizinger+2012, Kataoka+2013a}.
The porosities of observed FGRs are up to $\sim$50\% \citep{Daly+2018LPSC, Zanetta+2021}, and the realistic density ratio $\rho_{\rm ch}/\rho_{\rm rim}$ is 1 -- 2.
Although the rim thickness is an increasing function of $\rho_{\rm ch}/\rho_{\rm rim}$, the modest variation of the density ratio does not affect the rim thickness significantly.
When $\chi=$0.1 -- 0.8, the rim thickness for the case that $\rho_{\rm ch}/\rho_{\rm rim}=2$ becomes 1.4 -- 1.9 times thicker than that with $\rho_{\rm ch}/\rho_{\rm rim}=1$.

One observes that the rim thickness is now expressed by simple functions of both $a_{\rm ch}$ and $\chi$.
As described below, the value of $\chi$ at $\tau_{\rm rim}$ is not necessarily the same as the value of $\chi_{\rm local}$ in Equation (\ref{eq:chi_local}).
Note that in our model, it is idealized that compound aggregates stop the accretion of dust aggregates immediately at $\tau_{\rm rim}$ and start collisions between compound aggregates.
It is more realistic to consider that some of the compound aggregates would continue to accrete dust aggregates, that is, chondrules would acquire thicker rims.
In order to quantify this effect, we also estimate the rim thickness using $\chi_{\rm local}$.
This estimation assumes that all the chondrules accrete all the surrounding dust aggregates at the comparable rate, and hence the resulting rim thickness takes the maximum value;
if a small number of chondrules accrete all dust aggregates, then the thicker rim can be formed around these chondrules
\footnote{
	We note that both FGRs and igneous rims are not found around all chondrules \citep[e.g.,][]{Rubin1984, Simon+2018}.
	}.
Rims are, however, stripped in subsequent collisions.
Rim grains are lost in collisions between chondrules with rims \citep[][]{Umstatter+2019}, and high-velocity collisions induce chondrule ejection from aggregates \citep[][]{Arakawa2017}.
We ignore these subsequent events and give an upper limit estimation of rim thicknesses.

\subsection{Model calculations}

We perform simulations of the growth of chondrules and dust in the solar nebula to estimate $a_{\rm rim}$ using our dust growth model \citep{Matsumoto+2019}.
Our model computes $\chi$ as a function of $a_{\rm ch}$.

A complete explanation of our model can be found in \cite{Matsumoto+2019}, and hence only a short description is provided here.
Following the model, we assume that the growth stage is determined by the shortest timescale among the growth timescale of dust aggregates, the accretion timescale of dust aggregates onto compound aggregates (or chondrules), and the growth timescale of compound aggregates by their collisions.
The growth and accretion timescales are computed, using the disk and dust growth models described below.

We adopt a power-law disk model similar to the Minimum Mass Solar Nebula model \citep{Hayashi1981} to describe the solar nebula.
The surface density of the gas is given by
$\Sigma_{\rm g} = 2400\mbox{~g~cm}^{-2} \left( r/1~\mbox{au} \right)^{-3/2}$
\citep[e.g.,][]{Ida&Lin2004}.
The temperature is $T=280~\mbox{K} (r/1~\mbox{au})^{-1/2}$.
We consider a local region at 2~au where the accretion of dust onto chondrules took place.
The size of a chondrule ($a_{\rm ch}$) is a free parameter in this work.
The (initial) radius of a dust grain is assumed to be 0.25~$\mu$m.
The surface densities of chondrules and dust are given by
\begin{eqnarray}
	\Sigma_{\rm ch} &=& 10\mbox{~g~cm}^{-2} \chi_{\rm local} \left( \frac{r}{1~\mbox{au}} \right)^{-3/2},\\
	\Sigma_{\rm d} &=& 10\mbox{~g~cm}^{-2} (1-\chi_{\rm local}) \left( \frac{r}{1~\mbox{au}} \right)^{-3/2},
\end{eqnarray}
respectively.
The motions of chondrules, dust, dust aggregates and compound aggregates are governed by Brownian motion and turbulence \citep{Ormel&Cuzzi2007}.
Their (relative) velocities induced by the turbulence are expressed as a function of the turbulent $\alpha$, which is assumed to be $10^{-4}$, according to the paleomagnetic measurements of chondrules \citep[][]{Fu_R+2014, Hasegawa+2016b}.
Their Stokes numbers (St), which are the dimensionless stopping times, are $\St\lesssim \alpha \ll 1$.
These motions are well coupled with the gas since we consider the early growth phase,
where dust aggregates and compound aggregates are still small.
Their vertical scale heights are therefore equal to that of gas \citep[][]{Youdin&Lithwick2007}, which is given as $h_{\rm g}=c_{\rm s}/\Omega_{\rm K}$, where $c_{\rm s}$ is the sound speed and $\Omega_{\rm K}$ is the Keplerian frequency.
These small $\St$ aggregates undergo low velocity ($\lesssim 10\mbox{~cm~s}^{-1}$) collisions.
Comparing these collision velocities to the fragmentation velocity \citep[][]{Wada+2009,Arakawa2017}, these collisions are classified as merging events.
We, therefore, assume the perfect sticking, while lower sticking efficiencies are assumed in previous studies \citep[][]{Cuzzi+2004, Carballido2011}.
Realistically, the sticking efficiency depends on the internal structure of aggregates and is different in each collision.

Our previous work confirmed that the rim thickness does not strongly depend on disk parameters if all collisions are merging events.
This is because the rim thickness is given by the ratio of the growth timescale between the compound aggregates to the accretion timescale.
These timescales are given by their masses divided by their mass densities, collisional cross-sections, and their collision velocities.
In these quantities, only the collision velocities depend on the disk parameters.
The timescale ratio is not affected by most disk parameters when the collision velocities between any aggregates are given by the turbulence.

\subsection{Cross-section bias in the rim thickness measurements} \label{sect:bias}

\begin{figure}
	\centering
	\includegraphics[scale=.75]{section.pdf}
	\caption{
		Schematic sketch of the cross-section bias.
		It is assumed that a spherical chondrule with the radius of $a_{\rm ch}$ has an isotropic rim whose thickness is $a_{\rm rim}$.
		When this chondrule is cut at a horizontal plane given by an angle $\theta$, the observed chondrule is $a_{\rm ch,obs}$ in radius and has a rim with the thickness of $a_{\rm rim, obs}$.
		}
	\label{Fig:section}
\end{figure}

The size estimates of chondrules and of the surrounding rims are obtained mainly from measurements of two-dimensional cross-sections of meteorites.
It is known that the rim thickness is overestimated due to the cross-section bias.
We here consider this bias in the following simple way.

Figure \ref{Fig:section} shows the schematic picture of our estimation of the cross-section bias.
We assume that chondrules are spherical and homogeneously surrounded by rims.
When we observe a chondrule with a rim at a cutting horizontal plane that forms an angle $\theta$ with the vertical axis, both the observed chondrules radius ($a_{\rm ch,obs}$) and rim thickness ($a_{\rm rim, obs}$) are expressed as a function of $\theta$.
Due to curvatures, the chondrule radius is underestimated and the rim thickness is overestimated as $\theta$ becomes smaller.
A possible range of $a_{\rm ch,obs}$ is between 0 and $a_{\rm ch}$.
We take an average on $a_{\rm ch,obs}$ and $a_{\rm rim, obs}$ in the range of $ \cos^{-1}{(a_{\rm ch}/(a_{\rm ch}+a_{\rm rim}))} \leq\theta\leq 0.5\pi$.
The averaged values of $a_{\rm ch,obs}$ and $a_{\rm rim, obs}$ are not significantly different from the actual values of $a_{\rm ch}$ and $a_{\rm rim}$ since $a_{\rm ch,obs}$ decreases rapidly at $\theta\simeq\cos^{-1}{(a_{\rm ch}/(a_{\rm ch}+a_{\rm rim}))} $.

\subsection{Results}\label{sect:results}

\begin{figure*}
	\centering
	\begin{minipage}{0.45\hsize}
		\includegraphics[width=\textwidth]{ach_arim_F.pdf}
	\end{minipage}
	\begin{minipage}{0.45\hsize}
		\includegraphics[width=\textwidth]{chi_arim_F.pdf}
	\end{minipage}
	\\
	\begin{minipage}{0.45\hsize}
		\includegraphics[width=\textwidth]{ach_arim_C.pdf}
	\end{minipage}
	\begin{minipage}{0.45\hsize}
		\includegraphics[width=\textwidth]{chi_arim_C.pdf}
	\end{minipage}
	\caption{
	Thickness of the rim ($a_{\rm rim}$) is shown as functions of the chondrule diameter ($2a_{\rm ch}$, the left panels) and the chondrule mass fraction ($\chi$, the right panels).
	We show thicknesses of fine-grained dust rims in the upper panels and those of igneous rims in the bottom panels.
	The solid lines with the open squares are our results, which are the same between the upper and bottom panels.
	The filled symbols are measured rims found in the literature (see Table \ref{table:data}, we note that the \cite{Metzler+1992} data in this Table are not plotted since these values are not exact ones).
	When the measured values have ranges, the averaged values are plotted as symbols and the ranges are shown by the solid line segments.
	In the upper panels, we also plot the measured linear relationships between $a_{\rm rim}$ and $2a_{\rm ch}$ (see the dotted lines, Table \ref{table:data_Kc}).
	The linear relationships by \cite{Metzler+1992} are plotted between $2a_{\rm ch}=100$~$\mu$m and the larger edge of the chondrule diameter in Figure 17 of \cite{Metzler+1992}.
	We note that these measured linear relationships become curved lines on the log-log plot since they have interceptions.
	\newline
	The details of our results are as follows.
	{\em In the left panels}, the solid lines with the open squares are our simulation results where $\chi_{\rm local}=0.1, 0.3, 0.5$, and 0.7 from the top.
	Different colors represent different chondrule fractions at $\tau_{\rm rim}$ ($\chi$, see the colorbars).
	{\em In the right panels}, the hatched regions bounded by the open squares are given by our numerical results where $a_{\rm ch}=1000, 300, 100$, and 30~$\mu$m from the top.
	These hatched regions are defined by the results with $\chi= \chi_{\rm local}$ (the left edge) and those with $\chi$ at $\tau_{\rm rim}$ (the right edge).
	The dotted lines are obtained from Equation (\ref{eq:estimation_ach}).
	The dashed lines are drawn by applying the cross-section bias to Equation (\ref{eq:estimation_ach}).
	The colors of the symbols represent the values of $a_{\rm ch}$ (see the color bar on the right).
	}
	\label{Fig:a_ch_chi_arim}
\end{figure*}

\begin{figure}
	\centering
	\includegraphics[scale=.75]{ach_arim_line_mat.pdf}
	\caption{
		Each panel shows the measured linear relationships between the rim thickness ($a_{\rm rim}$) and the chondrule diameter ($2a_{\rm ch}$) of each meteorite (Table \ref{table:data_Kc}) in the solid line.
		The hatched regions are given by the rim thicknesses with $\chi= \chi_{\rm local}$ (the top edge) and those with $\chi$ at $\tau_{\rm rim}$ (the bottom edge).
	}
	\label{Fig:ach_arim_line_mat}
\end{figure}

The dependencies of the thickness of rims ($a_{\rm rim}$) on $a_{\rm ch}$ and $\chi$ are shown in Figure \ref{Fig:a_ch_chi_arim}.
In the left panels, the numerical results show that the thickness of rims is almost proportional to the chondrule diameter ($2a_{\rm ch}$), although there are some tiny fluctuations.
These fluctuations arise because the thickness also depends weakly on $\chi$ ($\chi^{-1/3}$, Equation (\ref{eq:estimation_ach})).

The relations between $a_{\rm rim}$ and $\chi$ are shown in the right panels (see the right edge of the hatched region).
For the comparison purpose, we also plot these thicknesses at $\chi = \chi_{\rm local}$ (see the left edge of the hatched region).
Since chondrules do not accrete all dust particles at the time of $\tau_{\rm rim}$ \citep[][]{Matsumoto+2019}, $\chi$ is larger than $\chi_{\rm local}$.
These two values are approximately the upper and lower limits of the chondrule abundance in finally formed chondrites (Section \ref{sect:model}).
These hatched regions become wider as $\chi_{\rm local}$ increases;
In large $\chi_{\rm local}$ cases, chondrules are the dominant component (over dust particles), and hence compound aggregates quickly undergo collisions between them.
As a result, compound aggregates can keep a high value of $\chi$ (i.e., chondrule-rich).
We also show the measured thicknesses of FGRs in the top panels and igneous rims in the bottom panels.
These observed rim thicknesses of both FGRs and igneous rims tend to increase as the chondrule size increases.
These data are summarized in Appendix (Table \ref{table:data}), where the estimated rim thicknesses are also summarized.

We consider that the measured rim thicknesses agree with our numerical results if the measured values (i.e., $\langle a_{\rm rim}\rangle$, $\langle 2a_{\rm ch}\rangle$, and $\langle \chi \rangle$) are located in the hatched regions (see also Table \ref{table:data} for specific values).
The thicknesses of FGRs in ordinary chondrites \citep{Bigolski2017_PhD} are well reproduced by our results.
The comparisons between our estimation and the observed linear relation are performed in Figure \ref{Fig:ach_arim_line_mat}, where the estimated rim thickness ranges are plotted on the chondrule diameter and rim thickness plane as the hatched regions.
The observed linear relations are consistent with our estimation in CM Y791198 and CM Y74662.
Most of their rim thicknesses are close to the estimated rim thicknesses from $\chi$ at $\tau_{\rm rim}$.
The rim thicknesses become larger than the estimated rim thicknesses from $\chi$ at $\tau_{\rm rim}$ as chondrule diameters decrease due to the interceptions of the observed linear relations.
In CM Kivesvaara, the measured rim thicknesses are consistent with our estimation in most of their chondrule diameter range, although the measured rim thicknesses are slightly larger than the estimated thicknesses around $100~{\rm \mu}$m diameters and slightly smaller than the estimated thicknesses at $1800~{\rm \mu}$m diameter.
However, the rim thicknesses around Murray chondrules are consistent only at $2a_{\rm ch}< 200~{\rm \mu m}$ and those around Murchison chondrules are always thinner than our estimation.

In contrast, the thicknesses of most igneous rims are larger than our results (Figure \ref{Fig:a_ch_chi_arim} and Table \ref{table:data}).
Only two data points of the rims in ALHA 77034 are consistent with our estimation, and the other rims are thicker than our estimation.
Because our estimation of the maximum rim thickness is based on the $\chi=\chi_{\rm local}$, there are two possible explanations: the cross-section bias; the subsequent acquisition of rim materials by mass transfer.
Large thicknesses of igneous rims can be reproduced by artificially adjusting the cross-section bias;
however, such adjustment makes the apparent chondrule size very small,
that is, 87\% of chondrules should be measured at $<0.4$ times smaller than their actual sizes.
This suggests that igneous rims are not simply formed through dust aggregate accretion since it is unlikely that almost all chondrules are cut at such small $\theta$ values.
Our results thus suggest that igneous rims would be made thicker by the subsequent processes such as remelting events.

\section{Discussion \& Conclusions}\label{sec:dis}

Rims around chondrules are often found in chondrites.
The thicknesses of fine-grained rims in CM chondrites is known to be proportional to the size of its host chondrule.
This linear relation was explained by theoretical studies.
However, the estimation of the thickness of rims in previous studies did not include the growth of dust, although the growth timescale of dust is shorter than the accretion timescale of dust onto chondrules.
We investigate the growth of chondrules and dust and estimate the thickness of the rims.
Our estimation is an upper limit estimation of rim thicknesses since we do not consider the subsequent rim erosion events in aggregate growth processes.
Our estimation agrees with the thicknesses of observed fine-grained rims in ordinary chondrites.
Some of the linear relationships observed in CM chondrites are consistent with our estimation, but the rim thicknesses around CM Murray and CM Murchison are smaller than our estimation.
The chondrules in these CM chondrites would lose their surrounding dust rims in the latter accretion processes such as collisions between chondrules with rims \citep[][]{Umstatter+2019} and collision-induced chondrule ejection \citep[][]{Arakawa2017}.
The thickness of the rims would be affected by the high-velocity collisions and subsequent fragments accretion \citep{Liffman2019}.
In addition, our model considers identical-sized chondrules.
The size distribution of the initial chondrules would affect their rim thicknesses.
These processes should be taken into account to consider the rim thicknesses around CM Murray and CM Murchison.

Besides, the rim thickness of CM chondrite is affected by parent body processes \citep[][]{Sears+1993, Trigo-Rodriguez+2006, Takayama&Tomeoka2012}.
Our results may imply that both the accretion of dust aggregates and the subsequent parent body processes are important to fully understand the rim thickness around CM chondrites.
More rim measurements and future simulations, where the above processes are included, provide more detailed views of the rim formation.

In contrast, the thicknesses of igneous rims are larger than our estimation.
This suggests that igneous rims would not be simply converted from fine-grained rims but get larger in their formation or subsequent processes.
Igneous rims are formed via remelting events, and they may be formed from the materials in chondrule surfaces and their surrounding rims \citep[e.g.,][]{Krot&Wasson1995}.
The accretion or collisions in these remelting events is one possible explanation of why igneous rims are thick.
The accretion of molten aggregates onto surficial melting chondrules would make thicker igneous rims.
These accretion events would be similar to the formation events of compound chondrules \citep[e.g.,][]{Wasson+1995}.
The size evolution of rims in the remelting events would be the key to understanding the thickness of igneous rims.
A more comprehensive study of the thickness of chondrule rims (both fine-grained and igneous) will help improve our modeling and understanding of the formation mechanism of the rims, which shed light on the origin of the solar system.

\section*{Acknowledgments}

We thank John N. Bigolski and the anonymous referee for constructive comments.
Y.H. is supported by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.