\section{Introduction}
\label{sec:Introduction}
\IEEEPARstart{E}{lectronic} Health Records (EHRs) are increasingly large and varied collections of health information
about patients. However, it is difficult for clinicians to retrieve information from EHRs efficiently and effectively given today's
limitations of user interfaces. Clinicians often operate under time pressure, and must invest significant amounts of time and
effort in retrieving information from EHRs. While using a search function can provide a useful alternative to browsing through
a patient's record, searching repeatedly for the same or similar information on similar patients is repetitive and cumbersome.
Under such circumstances,  there is a critical need to build effective systems that can generate accurate search
term recommendations for clinicians.

In this manuscript, we tackle the problem of next search term recommendation.
Given each patient's previous encounters and previous search terms that clinicians conducted on the patient, the objective of the next search
term recommendation problem is to recommend information items to clinicians that are most relevant to the patient to help diagnosis.
We assume clinicians conduct search over the patient EHR data in order to manually identify such information
items, while our method is able to proactively identify such information items and present them to the clinicians so as to
save the time and effort of manual search.

Recommender systems (RS) aim to recommend the most relevant items to a user  by
learning user preferences from their previously interacted items (e.g., products, information items).
Traditional RS techniques, such as collaborative filtering (CF)
~\cite{ricci2011introduction, Ning2015},
have been widely used to recommend the top-$N$ most relevant items to each user -- the so-called top-$N$ recommendation.
In this manuscript, we adapt traditional CF techniques and consider patients as users,
search terms as items, and propose a RS model to recommend the next search term on each patient, and to present the
terms to the clinicians.
While traditional RS techniques usually recommend the next items only based on users' previously
interacted items, the next search term recommendation can be a different problem as a search term on a patient
may have  strong relations with two factors: previous search terms (e.g., for chronic diseases)
and the patient's previous encounters (e.g., for most recent health concerns).

Based on this intuition, we propose a model named
\underline{H}ybrid \underline{C}ollaborative \underline{F}iltering \underline{M}odel using encounter information, denoted as \mbox{$\mathop{\text{HCFM}}\limits$}\xspace, to recommend
the next search term on a patient based on
previous search terms on the patient and his/her previous encounters.
Given the ICD codes in a patient's encounters and the previous search terms on each patient,
in \mbox{$\mathop{\text{HCFM}}\limits$}\xspace, we first calculate the co-occurrence frequency between each ICD code and each search term, where
a co-occurrence is considered when a term is searched right after the ICD  code has been
assigned to the patient.
The \mbox{$\mathop{\text{HCFM}}\limits$}\xspace model then recommends terms that have high co-occurrence frequencies with the most recent ICD  codes
and are highly relevant to the most recent search terms.
In addition, we propose another model that
uses only ICD-search term co-occurrence patterns to recommend the next search term for clinicians.
This method  recommends the terms that are highly relevant to all the previous encounters of a patient,
and is named
\underline{C}o-\underline{o}ccurrence \underline{P}attern based recommendation \underline{M}ethod, denoted as \mbox{$\mathop{\text{CoPM}}\limits$}\xspace.
The experimental results show that the proposed models outperform all the state-of-the-art
methods on top-$N$ search term recommendation on different real datasets.
We have also conducted comprehensive parameter studies
that enable important insights on the two factors.

\section{Related Work}
\label{sec:related}
The most relevant work to our method is from recommender systems (RS) research field.
Recently, a lot of research in RS has been conducted to tackle top-$N$ recommendation problems (i.e., to recommend
$N$ most relevant items).
Specifically, factorized personalized Markov chains
(FPMC)~\cite{rendle2010factorizing} generate recommendations for the next item using Markov
chains, which capture item-item transition relations.
Hidasi \mbox{\emph{et al.}}\xspace~\cite{hidasi2015session,hidasi2018recurrent} adapted deep learning
techniques, and used gated recurrent  units (GRUs) to model the dynamics of a
user's preferences.
Kang \mbox{\emph{et al.}}\xspace~\cite{kang2018self} developed a self-attention based sequential
model (SASRec) to capture the  most informative items in a user's previously
interacted items for top-$N$ recommendation.
Tang \mbox{\emph{et al.}}\xspace~{\cite{tang2018personalized}} developed
a convolutional sequence embedding recommendation model (\mbox{$\mathop{\text{Caser}}\limits$}\xspace)
by adapting multiple convolutional filters on the most recent interacted items of each user
to model sequential features in a user's historical interactions.
Ma \mbox{\emph{et al.}}\xspace~\cite{ma2019hierarchical} developed a hierarchical gating network (\mbox{$\mathop{\text{HGN}}\limits$}\xspace)
that adapts gating mechanisms over users' historical interacted items to identify
important items and their features  for top-$N$ recommendation.

Very limited work has been conducted on next search term recommendation problem for clinical applications.
Fan \mbox{\emph{et al.}}\xspace~\cite{fan2019improving} have developed several CF
models to tackle this problem.
Among these models, Transition-Involved Patient-Term-Similarity-based CF Scoring (\mbox{$\mathop{\text{TptCF}}\limits$}\xspace)
recommends the next search terms for clinicians on a patient based on the patient-patient similarities,
term-term similarities and term-term transition relations.
The difference between \mbox{$\mathop{\text{TptCF}}\limits$}\xspace and our model is that our model generates top-$N$ recommendations
using both previous search terms and previous encounters, while \mbox{$\mathop{\text{TptCF}}\limits$}\xspace only uses
previous search terms.

\section{Definitions and Preprocessing}
\label{sec:notations}

\subsection{Definitions and Notations}
\label{subsec:definitions}
Figure \ref{fig:matching} presents the  data preprocessing protocols and
Table~\ref{tbl:notations} presents the key notations used in this manuscript.
Formally,
the terms searched on each patient will be sorted chronologically.
The sequence of patient $p$'s sorted search terms is denoted as $S_p$, and
the subsequence of $S_p$ from
the $i$-th search to the $j$-th search is denoted as $S_p(i, j)$.
For simplicity, we store in $S_p$ the indices of search terms from a universal dictionary instead of search terms themselves.
Similar to search terms, the encounters of each patient will be sorted chronologically.
The sequence of patient $p$'s sorted encounters is denoted as $C_p$, and
the subsequence of $C_p$ from
the $i$-th encounter to the $j$-th encounter are denoted as $C_p(i, j)$.
For each patient, each search will be matched using timestamps to his/her most recent, previous encounter that the search
immediately follows (indicated by the green arrows in Figure~\ref{fig:matching}).
Please note that such matching does not necessarily mean that the searches happen during
the matched encounters,
or they are triggered or induced by the encounters, but it only indicates the temporal proximity.

\begin{figure}[!t]
\vspace{-5pt}
\centering
    \includegraphics[width=1\linewidth]{plots/matching_updated.pdf}
    \caption{{Data Preprocessing Protocols}}
    \label{fig:matching}
    \vspace{-5pt}
\end{figure}

For each patient, an ICD code may appear in his/her multiple encounters.
We denote the encounters of patient $p$ that contain ICD code $c$ as $C_p(c)$.
Similarly, a term can be searched for multiple times on a patient.
We denote the encounters of patient $p$ that each search term $s$ is matched to as $C_p(s)$.
We use the indices of ICD codes or search terms to refer them correspondingly in the sequences.

\subsection{Identifying Sessions}

In our dataset (discussed later in Section~\ref{sec:materials:dataset}), the time difference
between consecutive searches in a search sequence may vary from minutes to years, but
unfortunately session information is not always logged.
Therefore, we segment searches into sessions based on their timestamps using a sliding window
with size of three months. Therefore,
if the time interval between two consecutive searches are less than three months, the two search terms will be
grouped into one session.
To generate recommendations for a search,
we  only use the searches that are in a same session at the time when the recommendation
is to be made.
This is analogous to the real patient visit scenario in which clinicians search within the context of the
 current visit or recent searches.
Figure \ref{fig:matching} shows the search sessions with blue parentheses.
Please note that the sessions here are defined as a natural segmentation of search/encounter clusters, and they are
not the same as the sessions in internet connections.

\input{tables/notations.tex}

\section{Methods}
\label{sec:method}

Our method has the following steps.
First, we learn to numerically represent each ICD code and each search term so that the representations
can be used to calculate recommendation scores. Such representations will encode how ICD codes and
search terms co-occur in our data (Section~\ref{sec:method:learning}).
Then, we calculate a recommendation score for each term for each patient.
The recommendation score will be based on the previous search terms on the patients,
and the previous encounters of the patients.
We rank the recommendation scores of all terms for each patient and recommend top-$N$ terms
to the clinician seeing the patient (Section~\ref{sec:method:recommendation:EiCF}).

\subsection{Learning from Previous Search Terms and Encounters}
\label{sec:method:learning}

\subsubsection{Constructing an ICD Code-Search Term Co-Occurrence Matrix}
\label{sec:method:learning:construct}
We hypothesize that each search term on a patient is highly related to the patient's most recent
encounters, that is, given an ICD  code that is assigned to a patient, and thus the corresponding diagnoses,
certain terms that are related to the encounters/ICD codes are more likely to be searched next.
For example, if a patient was assigned an ICD code ``588.81: secondary hyperparathyroidism
(of renal origin)" in his/her recent encounters, terms such as potassium levels that are highly related to hyperparathyroidism
have high probability to follow compared to the case of ICD code ``786.2: Cough"
with which potassium level provides little information.
Thus, co-occurrence frequencies between ICD codes and search terms can provide
useful information to predict search terms that will follow;
given the ICD  codes in a patient's recent encounters, the terms with high co-occurrence
frequencies with these ICD  codes among many patients are more likely to be searched next, and thus should be
recommended.

Following the above intuition, we first calculate the co-occurrence frequency between
each ICD  code and each search term
by counting how many times in total the term has been searched after the ICD code in encounters.
We construct a matrix $A \in \mathbb{R}^{n \times m}$ to store such co-occurrence frequencies,
where $n$ is the number of all involved ICD  codes and $m$ is the number of all search terms.
We hypothesize that clinicians tend to search information from recent encounters of a
patient. Thus,
a search term is more likely to be induced from its recent ICD codes compared to
those from long time ago.
Based on this, we emphasize the recent encounters using
a time-decay parameter, and calculate the ICD code-search term co-occurrence frequencies $a_{cs}$ as follows:
\begin{equation}
    \label{eqn:cooccur}
    \begin{aligned}
    a_{cs} = \sum_{p=1}^{l}  \   \sum_{{e_c} \in C_p(c)} \   \sum_{{e_s} \in C_p(s)}  \lambda^{i(e_s)-i(e_c)} \mathbbm{1} (i(e_s) \geq i(e_c)) ,
    \end{aligned}
\end{equation}
where $e_s$ and $e_c$ are two encounters;
$l$ is the total number of patients,
$\lambda \in (0,1)$ is the time-decay parameter (in our experiments, $\lambda=0.5$);
$\mathbbm{1}(x)$ is the indicator function ($\mathbbm{1}(x) = 1$ if $x$ is true, otherwise, $\mathbbm{1}(x) = 0$);
$i(e_s)$ and $i(e_c)$ are the indices of encounter $e_s$ and encounter $e_c$, respectively, in patient $p$'s encounter
sequence $C_p$.
When calculating the co-occurrence frequencies between ICD code $c$ and term $s$, we only consider the
cases in which  term $s$ is searched after ICD code $c$ or in the same encounter as $c$
(i.e., $\mathbbm{1} (i(e_s) \geq i(e_c))$).
Please note that $a_{cs}$ is not a probability value and can have values greater than 1; larger $a_{cs}$ indicates
more likely that ICD code $c$ and search term $s$ co-occur.

\subsubsection{Learning ICD Code and Search Term Representations}
\label{sec:method:learning:latent}
Note that the co-occurrence matrix $A$ as constructed above is typically sparse because most of the ICD codes are not
observed as co-occurring with most of the search terms.
In order to capture the underlying relations between each ICD code and each search term that are not directly observed
in $A$,
we use a matrix factorization method~\cite{koren2009matrix} to learn the representations of ICD codes and search terms
which together produce what we observe in matrix $A$ and recover what we do not observe in $A$.
Specifically, we factorize $A$ into two low-rank matrices,
$U \in \mathbb{R}^{n \times d}$ and $V \in \mathbb{R}^{m \times d}$ ($d < \min(n, m)$),
representing the ICD  codes and search terms, respectively.
Particularly, each row of matrix $U$, denoted as $\mathbf{u}_c$, represents
the ICD  code $c$,
and each row of matrix $V$, denoted as $\mathbf{v}_s$, represents
the search term $s$.
Thus, all ICD  codes and search terms are represented by size-$d$ latent vectors and can be learned
from matrix $A$.
The co-occurrence ``chance" between ICD  code $c$ and search term $s$ can be estimated as follows:
\begin{equation}
    \label{eqn:factorize_cooccur}
    \hat{a}_{cs} = \mathbf{u}_c \mathbf{v}_s^{\mathsf{T}},
\end{equation}
where $\hat{a}_{cs}$ is the estimation of ${a}_{cs}$.
To learn the representations of each ICD code and each search term,
we formulate the following optimization problem:
\begin{equation}
    \label{eqn:factorize_opt}
    \underset{U, V}{\min}\  \|A - UV^{\mathsf{T}}\|_F^2 + \frac{\gamma}{2}(\|U\|_F^2 + \|V\|_F^2),

\end{equation}
where $U = [\mathbf{u}_1; \mathbf{u}_2; \cdots; \mathbf{u}_n]$,
$V = [\mathbf{v}_1; \mathbf{v}_2; \cdots; \mathbf{v}_m]$,
$\gamma$ is the weight for regularization term;
$\|.\|_F$ is the Frobenius norm, and the regularization on the Frobenius norm restricts large values in $U$ and $V$.
We  solve this problem using alternative gradient descent method \cite{koren2009matrix}.
Other representation learning methods include deep learning based methods~\cite{tang2018personalized,kang2018self}.
However, given the very sparse nature of our data (Section~\ref{sec:materials:dataset}), we do not use the deep learning
methods which typically require a lot of data for model training. Instead, matrix factorization based method as in
Equation~\ref{eqn:factorize_cooccur} has been demonstrated very effective in learning from sparse data~\cite{koren2009matrix}.

\subsection{\mbox{$\mathop{\text{HCFM}}\limits$}\xspace: Hybrid Collaborative Filtering Model}
\label{sec:method:recommendation:EiCF}

Our hybrid collaborative filtering model using encounter information,
denoted as \mbox{$\mathop{\text{HCFM}}\limits$}\xspace, generates recommendations for the next search term for clinicians on each patient
using two factors:
(i) the previous search terms on the patient, and
(ii) the previous encounters of the patient.
Each of two factors will be used to calculate a recommendation score for each term for the next
search on a patient (Section~\ref{sec:method:recommendation:EiCF:prevTerm} and Section~\ref{sec:method:recommendation:EiCF:prevEnc}),
and the two recommendation scores will be combined into a final recommendation score (Section~\ref{sec:method:recommendation:EiCF:rec}).

\subsubsection{{Recommendation Score from Previous Search Terms}}
\label{sec:method:recommendation:EiCF:prevTerm}

Previous studies from RS field have shown that more recent items will provide more
pertinent information to recommend the next item \cite{tang2018personalized,ma2019hierarchical}.
We hypothesize similarly that most recent search terms provide more pertinent information to recommend the
next search term.
Based on this, we recommend the next term using the
most recent search terms on a patient.
We aggregate the information from the most recent $m_s$ search terms in the current search session
by calculating the mean values of their
latent feature representations as follows:
\begin{equation}
\begin{aligned}
    \label{eqn:meanprevterm}
    \mathbf{m}_p = \frac{1}{m_s}\sum_{i \in S_p(n_p-m_s,n_p)} \mathbf{v}_i,
\end{aligned}
\end{equation}
where
$n_p$ is the number of all  search terms on patient $p$ at the time when the recommendation is to be made;
$m_s$ is the number of the most recent search terms that are used for recommendation
($m_s$ is a fixed number in our experiments);
$\mathbf{m}_p  \in \mathbb{R}^{1 \times d}$ is the aggregated representation from previous $m_s$
search terms on patient $p$;
The recommendation score of term $s$ on patient $p$ based on previous search terms
is calculated as the dot-product similarity between $\mathbf{m}_p$ and $\mathbf{v}_s$ as follows:
\begin{equation}
    \label{eqn:prevterm}
    {x}_{ps} = \mathbf{m}_p \mathbf{v}_s^{\mathsf{T}}.
\end{equation}

\subsubsection{{Recommendation Score from Previous Encounters}}
\label{sec:method:recommendation:EiCF:prevEnc}
To calculate the recommendation score for each term, we also use  the information from the most recent
$m_c$ encounters of each patient.
Here, we hypothesize that clinicians tend to search within the context of the recent searches;
ICD codes that induce the most recent searches  are more likely to induce the next search, and thus
should be emphasized to recommend the next search term.
Based on this,
we learn an importance weight on each ICD code $c$ for each patient $p$.
The importance weight is calculated as
the normalized dot-product similarity between each ICD  code and the most recent $m_s$ search terms
as follows:
\begin{equation}
\begin{aligned}
    \label{eqn:wjs}
     w_{pc}
             =  \frac{\exp(\mathbf{u}_c \mathbf{m}_p^{\mathsf{T}} )}{{\sum_{e' \in C_p(l_p-m_c,l_p)}  \  \sum_{c' \in e'} \exp(\mathbf{u}_c' \mathbf{m}_p^{\mathsf{T}} )}}
\end{aligned}
\end{equation}
where $\mathbf{m}_p$ is calculated as in Equation~\ref{eqn:meanprevterm};
$l_p$ is the number of all encounters of patient $p$ at the time when the recommendation is to be made;
$m_c$ is the number of the most recent encounters that are used for recommendation
($m_c$ is a fixed number in our experiments);
$e'$ is an encounter in $C_p(l_p-m_c, l_p)$, and $c'$ is an ICD code in $e'$.
The recommendation score of term $s$ on patient $p$ based on previous encounters is calculated as follows:
\begin{equation}
    \label{eqn:prevenc}
    {y}_{ps} = \sum_{e \in C_p(l_p-m_c, l_p)} \  \sum_{c \in e}  w_{pc} \mathbf{u}_c \  \mathbf{v}_s^{\mathsf{T}},
\end{equation}
where $e$ is an encounter in $C_p(l_p-m_c, l_p)$, and $c$ is an ICD code in $e$.

\subsubsection{{Combination of Recommendation Scores} }
\label{sec:method:recommendation:EiCF:rec}
The recommendation scores of term $s$ on patient $p$ calculated as above are then combined into a final recommendation
score as follows:
\begin{equation}
    \label{eqn:recommendation}
    {r}_{ps} = \alpha \  {x}_{ps} + (1-\alpha) \  {y}_{ps}
\end{equation}
where $\alpha \in [0,1]$ is a pre-defined weight for the two factors.
In Equation~\ref{eqn:recommendation}, $\alpha$=$1$ indicates
that only previous search terms are used for recommendation,
and $\alpha$=$0$ indicates that only  previous encounters are used for recommendation.
The recommendation scores of all  terms will be sorted, and the  terms with
top-$N$ scores will be recommended for the next search.

\subsection{\mbox{$\mathop{\text{CoPM}}\limits$}\xspace: Co-occurrence Pattern-based Model}
\label{sec:method:recommendation:CMF}

We propose another method, denoted as \mbox{$\mathop{\text{CoPM}}\limits$}\xspace, that
uses only ICD-search term co-occurrence patterns (Section~\ref{sec:method:learning})
to recommend the next search term for clinicians.
The difference between \mbox{$\mathop{\text{HCFM}}\limits$}\xspace and \mbox{$\mathop{\text{CoPM}}\limits$}\xspace is that
\mbox{$\mathop{\text{HCFM}}\limits$}\xspace recommends the search terms that are most relevant to the most recent
search terms and the most recent encounters,
whereas \mbox{$\mathop{\text{CoPM}}\limits$}\xspace recommends the search terms that are most relevant to
all the previous encounters of a patient.
Both methods represent ICD codes and search terms using the representation matrices, $U$ and $V$, respectively,
as learned based on Equation \ref{eqn:factorize_opt}.
The \mbox{$\mathop{\text{CoPM}}\limits$}\xspace model will serve as a baseline method in our experiments.

Specifically, we aggregate the information of all ICD codes
in the patient's previous encounters  to calculate the recommendation score for each search term.
We hypothesize that more recent ICD codes are more likely to induce the next search term. Therefore,
we emphasize the encounters/ICD codes that are closer to a search to generate the term
recommendations for this search  using a time-decay parameter.
Specifically,
the recommendation score for term $s$ on  patient $p$ is calculated as follows:
\begin{equation}
    \label{eqn:factorize_cooccur_pred}

     r_{ps} =  \sum_{e \in C_p(1,l_p)} \sum_{c \in e} \sigma^{i(e_s)-i(e)} \mathbf{u}_c \mathbf{v}_s^{\mathsf{T}},
\end{equation}
where $e$ is an encounter in $C_p(1, l_p)$, and $c$ is an ICD code in $e$;
$e_s$ is the encounter that current search is matched to, that is, the most recent, previous encounter
at the time that the recommendation is to be made;
$i(e_s)$ and $i(e)$ are the indices of encounter $e_s$ and encounter $e$, respectively;
$\sigma \in (0,1)$ is the time-decay parameter (in our experiments, $\sigma$=$0.5$).
Note that the time-decay parameter $\sigma$ here indicates how long ago each encounter occurred before the
time of recommendation, whereas the time-decay weights $\lambda$ in Equation~\ref{eqn:cooccur} indicates
the temporal proximity between an encounter and a search term.
Thus, the two time-decay parameters represent different information in the model.
Finally, we  sort the recommendation scores of all the terms and recommend the
top-$N$ terms for the next search.

\section{Materials}
\label{sec:materials}

\subsection{Baseline Methods}
\label{sec:materials:baseline}
We compare the \mbox{$\mathop{\text{HCFM}}\limits$}\xspace and \mbox{$\mathop{\text{CoPM}}\limits$}\xspace with the following state-of-the-art baseline methods.
Some of the methods are from recommender systems (RS) research, in which
they have users and items.
To use these methods on the next search term recommendation task,
we consider patients analogous to users, and search terms analogous to items in our
experiments.

\begin{itemize}[noitemsep,nolistsep,leftmargin=*
\item \textbf{Transition-Involved Patient-Term-Similarity-based CF Scoring (\mbox{$\mathop{\text{TptCF}}\limits$}\xspace)}
\cite{fan2019improving}.
The \mbox{$\mathop{\text{TptCF}}\limits$}\xspace model generates recommendations for the next search term using two factors:
(1) patients' similarities and search terms' similarities (similarity-based scoring) and
(2) search term dynamic transitions (dynamics-based scoring).
For each patient, \mbox{$\mathop{\text{TptCF}}\limits$}\xspace will calculate a recommendation score for each search term using
similarity-based scoring and  dynamics-based scoring, and recommend the top search terms
with the highest scores.
The detailed explanation of the \mbox{$\mathop{\text{TptCF}}\limits$}\xspace method will be presented in the supporting document.

\item \textbf{Personalized top-$N$ (\mbox{$\mathop{\text{PTN}}\limits$}\xspace)}
The \mbox{$\mathop{\text{PTN}}\limits$}\xspace model
recommends the next item using the most frequently interacted items in each user's history.
In our experiments,
the terms are ranked based on how many times they are searched on each patient, and the patient's
top-$N$ most frequently searched terms will be recommended to the patient.
\item \textbf{Hierarchical Gating Network (\mbox{$\mathop{\text{HGN}}\limits$}\xspace)}~\cite{ma2019hierarchical}.
\mbox{$\mathop{\text{HGN}}\limits$}\xspace selects important items and their important features by adapting gating
mechanisms \cite{ma2019hierarchical} over users' historical interacted items.
To generate recommendations for the next items,
\mbox{$\mathop{\text{HGN}}\limits$}\xspace uses the identified important items  and important features, and
users' general preferences to calculate a recommendation score of each item for a user.
\mbox{$\mathop{\text{HGN}}\limits$}\xspace is a state-of-the-art sequential recommendation method.
\item \textbf{Hybrid Associations Model (\mbox{$\mathop{\text{HAM}}\limits$}\xspace)}~\cite{Peng2020HAM}.
\mbox{$\mathop{\text{HAM}}\limits$}\xspace models users' long-term preferences from all their historical items,
and models users' short-term preferences from their most recent items.
The short-term preferences contain both high-order and low-order association patterns among items.
Both the long-term and short-term preferences are used to recommend the next
items.
\mbox{$\mathop{\text{HAM}}\limits$}\xspace has been demonstrated as the state-of-the-art algorithm for sequential top-$N$ recommendation.

\end{itemize}

\subsection{Datasets}
\label{sec:materials:dataset}

The data used in our experiments are collected from Eskenazi Health organization in Indiana, US,
from 04/2013 to 05/2016.
This dataset contains 13,934 patients, their 1,377,381 encounters, 9,565 valid ICD 9 codes
and 7,215 unique search terms.
Among the 7,215 unique search terms, we remove the irregular search terms,
such as numbers and punctuations, and also infrequent search terms that appear only once,
and map the misspelled terms to their most similar
terms (e.g., ``adimssion" is mapped to ``admission").
We also remove the infrequent patients and keep those who have at least two search terms
and at least three encounters.
Table \ref{tbl:dataset} presents the statistics of the dataset after data preprocessing.
Overall, 2,955 patients and 2,101 unique search terms are retained in the
dataset.
On average, each
patient has 10.22  searches
and 173.26 encounters.
Each term has been searched for 14.37 times on average over all the patients, and each encounter has 2.09 ICD  codes on average.
Our study is conducted under Protocol
\#1612682149 ``Supporting information retrieval in the ED through collaborative filtering"
approved by the Indiana University Institutional Review Board (IRB).

\input{tables/statistics}

We present the distribution of search sequences' lengths in the dataset in Figure~\ref{fig:seq_len}.
As Figure~\ref{fig:seq_len} shows, there are many more short sequences in
the dataset  than long sequences.
Figure~\ref{fig:term_freq} presents the distribution of the number of unique terms for each patient. As it indicates,
most of the search terms are not frequent.
On average, each patient has 7.00 unique search terms.
Moreover, as presented in Section \ref{subsec:definitions}, we group search terms into
 sessions.
Overall, there are 3,488 sessions identified from the dataset.
On average, each patient has 1.18 sessions, and each session has 8.66 search terms.
It is notable that search sequences are typically very short, and the number of unique search terms
per patient is very small. This will make the recommendation problem difficult because of the
sparsity of the available data.
Table \ref{tbl:freq-terms} lists the most frequently searched terms,
where frequency is calculated based on how many times each term is searched
on all patients.

\begin{figure}[h!]
\centering
    \centering
    \includegraphics[width=0.75\linewidth]{plots/seq_len.png}

    \caption{Distribution of Search Sequence Lengths}
    \label{fig:seq_len}
\end{figure}

\begin{figure}[!h]
    \centering
    \includegraphics[width=0.75\linewidth]{plots/uniq_term_patient.png}

   \caption{Distribution of Unique Terms per Patient}
    \label{fig:term_freq}
\end{figure}

\input{tables/top_20_freq_terms}

\subsection{Experimental Protocols}
\label{sec:materials:setting}
\begin{figure}[!h]
    \includegraphics[width=\linewidth]{plots/setting.pdf}
    \caption{CUTOFF Experimental Protocol}
    \label{fig:setting}
\end{figure}
\input{tables/cutoff_stats}

We use the following experimental protocol to evaluate our methods on the dataset.
In this protocol, all the search sequences are split by a same cut-off time.
Any searches before the cut-off time are split into training set, and any searches after the cut-off time
are split into testing set. The models are trained using only the training set, for example, the co-occurrence
matrix (Equation \ref{eqn:cooccur}) is constructed only using the search terms and encounters in the training set.
This protocol is referred to as cut-off cross validation, denoted as CUTOFF.
Figure \ref{fig:setting} demonstrates the CUTOFF experimental protocol.

We selected four cut-off times: 10/01/2013, 04/01/2014, 10/01/2014 and 04/01/2015.
These cut-off times are selected because sufficient search terms and encounters from a majority of the search sequences
are retained in training set before the cut-off time and meanwhile sufficient search sequences have testing terms after the
cut-off time. After the split, the statistics for the training and testing data of four cut-off times that we choose
are presented in Table \ref{tbl:dataset_cutoff}. This CUTOFF setting is close
to the realistic scenario, that is, all the data before a certain time should be used to predict information after that time.
However, a shortcoming of CUTOFF is that many early search sequences may not have testing terms, and many late
search sequences will not have anything in training set. Sequences that do not have testing terms are still used
to train models. Sequences that do not have training terms are not used. For those sequences which have terms
after the cut-off time, only the first one of the terms after the cut-off time will be used for evaluation.

\subsection{Evaluation Metrics}
\label{sec:materials:metric}
We use hit rate at $k$ (HR@$k$) to evaluate the different methods.
HR@$k$ measures the ratio of the patients whose ground-truth next search term is within the
top-$k$ recommendations, over all patients.
Higher HR@$k$ indicates better recommendation performance.

\subsection{Parameter Tuning}
\label{sec:materials:parameter}
We use grid search to tune four parameters:
$m_s$ (Equation \ref{eqn:meanprevterm}), $m_c$ (Equation \ref{eqn:wjs}, \ref{eqn:prevenc}),
 $\alpha$ (Equation \ref{eqn:recommendation}) and latent dimension $d$.
We will present the best results among different parameter settings in Section \ref{sec:results}.

\section{Results and Discussions}
\label{sec:results}

\subsection{Overall Performance}
\label{sec:results:overall}

We compare \mbox{$\mathop{\text{HGN}}\limits$}\xspace, \mbox{$\mathop{\text{HAM}}\limits$}\xspace,  \mbox{$\mathop{\text{TptCF}}\limits$}\xspace, \mbox{$\mathop{\text{PTN}}\limits$}\xspace, \mbox{$\mathop{\text{CoPM}}\limits$}\xspace and \mbox{$\mathop{\text{HCFM}}\limits$}\xspace in our experiments.
Table \ref{tbl:results20150401} , \ref{tbl:results20141001}, \ref{tbl:results20140401} and \ref{tbl:results20131001}
present the best performance of each method in terms of HR@$k$ ($k \in \{1,2,3,4,5,10,20\}$)
for four cut-off times, respectively.
In Table \ref{tbl:results20150401}, we present different parameter settings of each method that achieve the best HR@$k$ results for
each $k$ value.
For simplicity, in Table \ref{tbl:results20141001}, \ref{tbl:results20140401} and \ref{tbl:results20131001},
we only present the parameter setting of each method that achieve the best HR@$5$.
We also present the improvement of \mbox{$\mathop{\text{HCFM}}\limits$}\xspace over the second best results in terms of all HR metrics in the tables.
In addition, we add the experimental results for the \mbox{$\mathop{\text{HCFM}}\limits$}\xspace model when $m_s$ equals to the number of all previous
search terms in the current search session for each patient
(i.e., the row for \mbox{$\mathop{\text{HCFM}}\limits$}\xspace with parameter ``all" in Table~\ref{tbl:results20150401}).

\subsubsection{Overall Performance on Cut-off Time 04/01/2015}
\label{sec:results:overall:20150401}

\input{tables/results20150401}

Table \ref{tbl:results20150401} shows that the \mbox{$\mathop{\text{HCFM}}\limits$}\xspace model outperforms all the baseline methods
in terms of HR@$k$ with $k \in \{1,2,3,4,5,10\}$, and achieves the second best
result in terms of HR@$20$.
\mbox{$\mathop{\text{HCFM}}\limits$}\xspace achieves HR@$1$=$0.0917$, meaning that about 9\% of the patients have their ground-truth search terms correctly
recommended at the very top.
Although this HR@$1$ value is not high by itself,
it is actually substantially higher than random guessing, which has an expected HR@$1$=$\frac{1}{1,915}$=$5.22e^{-4}$
(1,915 search terms in training set for cut-off time 04/01/2015 as in Table~\ref{tbl:dataset_cutoff}).
That is, \mbox{$\mathop{\text{HCFM}}\limits$}\xspace is 175-fold better than random guessing on HR@$1$.
\mbox{$\mathop{\text{HCFM}}\limits$}\xspace achieves HR@$5$=$0.2569$, meaning that about 26\% of the patients have their ground-truth search terms
among top-5 recommendations. This is also substantially higher than random guessing, which has an expected
HR@$5$=$5.22e^{-4}\times 5$=$2.6e^{-3}$. That is, \mbox{$\mathop{\text{HCFM}}\limits$}\xspace is 98-fold better than random guessing on HR@$5$.

The second best method is \mbox{$\mathop{\text{CoPM}}\limits$}\xspace as it achieves the second best results  on
HR@$2$, HR@$3$, HR@$4$, HR@$5$ and HR@$10$, and the best results on HR@$20$.
The difference between \mbox{$\mathop{\text{HCFM}}\limits$}\xspace and \mbox{$\mathop{\text{CoPM}}\limits$}\xspace is that \mbox{$\mathop{\text{HCFM}}\limits$}\xspace uses the most recent search terms and the most recent encounters to
recommend the next search term, whereas \mbox{$\mathop{\text{CoPM}}\limits$}\xspace uses all the previous encounters to recommend the next search term.
However, as Table \ref{tbl:dataset_cutoff} shows, each patient has 167.88 encounters on average.
Since previous encounters that occurred long time ago may not contain the information that the clinicians want to search at the time of recommendation,
therefore, \mbox{$\mathop{\text{CoPM}}\limits$}\xspace may generate inaccurate search term recommendations using all previous encounters.

\mbox{$\mathop{\text{PTN}}\limits$}\xspace is slightly worse than \mbox{$\mathop{\text{CoPM}}\limits$}\xspace and \mbox{$\mathop{\text{HCFM}}\limits$}\xspace. This is probably because \mbox{$\mathop{\text{PTN}}\limits$}\xspace uses the most popular search terms
in a patient's entire searching history for recommendation,
and the popular terms that occurred long time ago may not be of interest to clinicians at the time of recommendation,
given the progression of the patient's health condition.
\mbox{$\mathop{\text{TptCF}}\limits$}\xspace has slightly worse performance than \mbox{$\mathop{\text{PTN}}\limits$}\xspace. This might be due to the data sparsity issue.
As shown in Table \ref{tbl:dataset_cutoff}, the dataset contains
2,627 patients and 1,915 unique search terms in the training set. However, each patient has only 8.22 search terms on average.
The data sparsity issue may cause \mbox{$\mathop{\text{TptCF}}\limits$}\xspace to learn inaccurate patient-patient similarities and term-term similarities, both of
which can lead to inaccurate search term recommendations.

\mbox{$\mathop{\text{HAM}}\limits$}\xspace and \mbox{$\mathop{\text{HGN}}\limits$}\xspace have the worst performance in this experiment.
In the RS settings, \mbox{$\mathop{\text{HGN}}\limits$}\xspace assumes that each item contributes differently to  next item recommendations, and  therefore
learns  importance weights for  items using gating mechanisms.
However, the data sparsity issue in our dataset may cause \mbox{$\mathop{\text{HGN}}\limits$}\xspace to learn inaccurate weights for items, and thus
leads to poor recommendation performance.
On the other hand, \mbox{$\mathop{\text{HAM}}\limits$}\xspace generates recommendations using three factors:
users' long-term preferences modeled from all previous items,
high-order association patterns modeled from a number of most recent previous items
and low-order association patterns modeled from few most recent previous items.
Our dataset may not contain the long-term preferences and the association patterns, and
thus \mbox{$\mathop{\text{HAM}}\limits$}\xspace may generate less meaningful recommendations and thus have
poor recommendation performance.

Table \ref{tbl:results20150401} also shows that \mbox{$\mathop{\text{HCFM}}\limits$}\xspace consistently improves over the second best method on all the
evaluation metrics except  HR@$20$. In real applications, we prefer good performance with small values $k$ in HR@$k$,
indicating that the correct recommendations are on very top.
Therefore, good performance on HR@$k$, $k<20 $ is of more significance
than that on HR@$20$.
In addition, \mbox{$\mathop{\text{HCFM}}\limits$}\xspace
has relatively high improvement (more than $10\%$) over the second best results in terms of
HR@$4$ and HR@$5$. This indicates that \mbox{$\mathop{\text{HCFM}}\limits$}\xspace can push the most relevant search terms on very top of the recommendation list.

\subsubsection{Overall Performance on Other Cut-off Times}
\label{sec:results:overall:other}

\input{tables/results20141001}

\input{tables/results20140401}

\input{tables/results20131001}

Table \ref{tbl:results20141001} and Table \ref{tbl:results20140401} presents the performance of the different methods
on cut-off time 10/01/2014 and 04/01/2014, respectively.
They show similar trends as those in
Table \ref{tbl:results20150401}, and \mbox{$\mathop{\text{HCFM}}\limits$}\xspace outperforms the baseline methods in terms of all HR metrics.
Further comparing Table \ref{tbl:results20150401}, Table \ref{tbl:results20141001} and
Table \ref{tbl:results20140401}, we
notice that all the methods tend to have higher HR results in Table  \ref{tbl:results20140401} (cut-off time 04/01/2014) than those in
Table \ref{tbl:results20141001} (cut-off time 10/01/2014) and Table \ref{tbl:results20150401} (cut-off time 04/01/2015).
This is probably because there are fewer search terms (i.e., 1,065) in the training set with cut-off time
04/01/2014 compared to those with cut-off-time 10/01/2014 (i.e., 1,493) and those
with cut-off time 04/01/2015 (i.e., 1,915), as shown in Table \ref{tbl:dataset}.
Note that in our experiments, the model will only recommend the search terms that have appeared in the training set,
and will not recommend new search terms.
Since generating recommendations from 1,065 search terms is
more likely to hit the accurate term by chance
than generating recommendations from 1,915 search terms.
Thus,
\mbox{$\mathop{\text{HCFM}}\limits$}\xspace is more likely to recommend search terms accurately with cut-off time 04/01/2014
than with cut-off time 04/01/2015.

Table \ref{tbl:results20131001} presents the performance of the different methods
on cut-off time 10/01/2013. It shows that  \mbox{$\mathop{\text{HCFM}}\limits$}\xspace achieves the best results
in terms of HR@$3$, HR@$4$, HR@$5$ and HR@$10$,
and achieves the second best results in terms of HR@$1$ and HR@$20$.
In this experiment, \mbox{$\mathop{\text{PTN}}\limits$}\xspace achieves the best HR@$1$, HR@$2$ and HR@$3$ results.
Furthermore, \mbox{$\mathop{\text{HAM}}\limits$}\xspace achieves the best result in terms of HR@$20$, and achieves the second best
results in terms of HR@$2$, HR@$3$, HR@$5$ and HR@$10$.
In this experiment,
there are very few patients and search terms in the training set and testing set, as shown in Table \ref{tbl:dataset}.
Due to the data sparsity, the representations of search terms and ICD 9 codes may not be well learned, and thus
it may cause both \mbox{$\mathop{\text{HCFM}}\limits$}\xspace and \mbox{$\mathop{\text{CoPM}}\limits$}\xspace  to have relatively poor recommendation performance, and do not achieve
significant improvement over the baseline methods.
Table \ref{tbl:results20131001} also shows that \mbox{$\mathop{\text{PTN}}\limits$}\xspace has the same values for all the HR metrics.
This is probably because there are sessions in which testing search terms do not appear in the training part
of the same sessions. Therefore, \mbox{$\mathop{\text{PTN}}\limits$}\xspace is not able to recommend the accurate search terms for such sessions
despite of the value of $k$.

Further comparing Table \ref{tbl:results20150401}, \ref{tbl:results20141001}, \ref{tbl:results20140401}
and \ref{tbl:results20131001},
we notice that \mbox{$\mathop{\text{HCFM}}\limits$}\xspace tends to significantly outperform the baseline methods when there is sufficient data
for training, and would
not perform very well when data is very sparse.
Overall, for all three datasets,  \mbox{$\mathop{\text{HCFM}}\limits$}\xspace consistently outperforms the baseline methods in terms of
HR@$3$, HR@$4$, HR@$5$ and HR@$10$, and
achieves improvement of at least 4.56\%, 7.28\%, 7.15\% and
5.54\% over the second best method, respectively.

\subsection{Parameter Study}
\label{sec:results:parameter}

For a parameter study, we evaluate the models with cut-off time 04/01/2015
since there is sufficient training and testing data with this cut-off time.
We choose the set of parameter as follows: $m_s$=$6$, $m_c$=$2$ and  $\alpha$=$0.2$,
since \mbox{$\mathop{\text{HCFM}}\limits$}\xspace
achieves the best HR@$2$, HR@$3$ and HR@$5$ with this set of parameters.
To conduct parameter study, we fix two out of three parameters and evaluate the model with
different values of the third parameter.
Table \ref{tbl:parameter-alpha}, \ref{tbl:parameter-lc} and \ref{tbl:parameter-ls} present the results
for parameter study on \mbox{$\mathop{\text{HCFM}}\limits$}\xspace on $m_s$, $m_c$ and $\alpha$, respectively.

\input{tables/parameter}

Table \ref{tbl:parameter-ls} shows that \mbox{$\mathop{\text{HCFM}}\limits$}\xspace achieves the best HR@$5$ result with $m_s \in [5,9]$
and the value of HR@$5$ decreases when $m_s$ increases or decreases from the range.
Recall that $m_s$ is the number of previous search terms that are used to recommend the next search term.
The result here indicates that
too few previous search terms may not provide sufficient information
and too many pervious search terms may provide irrelevant information for recommendation.
Table \ref{tbl:parameter-ls}  also shows that \mbox{$\mathop{\text{HCFM}}\limits$}\xspace achieves  the same HR@$1$ result with different
$m_s$ values. This indicates that increasing or decreasing $m_s$ values may not be able to push the most relevant
search terms to the very top of the recommendation list, but can rank more relevant search terms on top
since HR@$5$ is improved when choosing a proper $m_s$.

Table \ref{tbl:parameter-lc}  shows that \mbox{$\mathop{\text{HCFM}}\limits$}\xspace achieves  the same HR@$1$ result with different
$m_c$ values, achieves the best HR@$5$ with $m_c$$=$$2$ and achieves the best HR@$5$
with $m_c$$=$$2$ and $3$.
Recall that $m_c$ is the number of previous encounters that are used to recommend the next search term.
Table \ref{tbl:parameter-lc}   also shows that both HR@$5$ and HR@$10$
decrease as $m_c$ increases or decreases from 2 (or 3).
This indicates that changing $m_c$ values may not help \mbox{$\mathop{\text{HCFM}}\limits$}\xspace to push the most relevant search terms
on the very top of the recommendation list, but when $m_c$ is set to a small number (e.g., $m_c$=$2$),
\mbox{$\mathop{\text{HCFM}}\limits$}\xspace can achieve its best results in terms of HR@$5$ and HR@$10$.

Table \ref{tbl:parameter-alpha} shows that when $\alpha$ (Equation \ref{eqn:recommendation}) increases, HR@$1$ increases,
and HR@$10$ tends to decrease.
Furthermore, \mbox{$\mathop{\text{HCFM}}\limits$}\xspace achieves the best HR@$5$ results with $\alpha$ equal to $0.2$, $0.4$, $0.8$.
Recall that with $\alpha$ set to $1.0$, \mbox{$\mathop{\text{HCFM}}\limits$}\xspace only uses previous search terms to generate recommendations, and
with  $\alpha$ as $0$, \mbox{$\mathop{\text{HCFM}}\limits$}\xspace only uses previous encounters for recommendation.
This indicates that having more weights on previous search terms can increase HR@$1$ performance,
but hurt HR@$10$ performance,
whereas adding more weights to previous encounters can hurt HR@$1$ performance, but improve HR@$10$ performance.
This is probably because using previous search terms can push the most relevant
search terms to the very top of the recommendation list, thus, \mbox{$\mathop{\text{HCFM}}\limits$}\xspace can achieve
the best HR@$1$ result with $m_s$ equal to $1.0$.
However, only using previous search terms overlooks encounters' information, which may be highly relevant
to clinicians' searches.
Therefore,  \mbox{$\mathop{\text{HCFM}}\limits$}\xspace  can push more relevant search terms on top, and achieves the best HR@$10$ with $m_s$ equal to $0$.
When using both previous search terms and previous encounters, \mbox{$\mathop{\text{HCFM}}\limits$}\xspace is able to achieves the best HR@$5$ results.

\subsection{Analysis on Length of Search Sessions}
\label{sec:results:session_length}

\begin{comment}
\begin{table}
\vspace{-10pt}
\centering
\captionof{table}{statistics for lengths of sessions}
  \label{tbl:session_length}
      \begin{tabular}{
        @{\hspace{0pt}}l@{\hspace{2pt}}
        @{\hspace{2pt}}r@{\hspace{2pt}}
        @{\hspace{2pt}}r@{\hspace{2pt}}
        @{\hspace{2pt}}r@{\hspace{2pt}}
        }
        \toprule
        length of sessions & min   & max & mean \\
        \midrule
        top 20\% shortest sessions       & 2 & 2 & 2 \\
        top 20\%-40\% shortest sessions  & 3 & 5 & 3.95 \\
        mid 40\%-60\% sessions           & 6 & 8 & 6.76 \\
        top 20\%-40\% longest sessions   & 9 & 24 & 15.04 \\
        top 20\% longest sessions        & 26 & 202 & 62.75 \\
        \bottomrule
      \end{tabular}
\end{table}
\end{comment}

\input{plots/session_length}

We evaluate \mbox{$\mathop{\text{HCFM}}\limits$}\xspace on sessions of different sequence lengths
in order to better understand the  influence of the session lengths and information content on the
recommendation performance.
Specifically, we first divide the sessions in the testing set into five groups based on the session lengths.
Table \ref{tbl:session_length} presents the statistics of the lengths of search sessions in our experiments.
For each group, we present the min length, max length and average length of the sessions.
In this experiment, we include all the previous searches for the next search term recommendation,
therefore, we set $m_s$ as ``all" in \mbox{$\mathop{\text{HCFM}}\limits$}\xspace.
Given $m_s$ as ``all", $m_c$=$2$ and $\alpha$=$0.8$ enable the best
performance of \mbox{$\mathop{\text{HCFM}}\limits$}\xspace, as shown in Table \ref{tbl:results20150401}.
Thus, we select $m_s$=``all", $m_c$=$2$ and $\alpha$=$0.8$ as the parameters in our experiment.
In addition, $m_s$=6, $m_c$=$2$ and $\alpha$=$0.2$ in \mbox{$\mathop{\text{HCFM}}\limits$}\xspace achieve the best HR@2,
HR@3 and HR@5 results, as shown in Table \ref{tbl:results20150401}.
For comparison purposes, we also add this set of parameters in our experiment.
Finally, in order to understand the effects of parameter $\alpha$,
we add another two parameter settings for \mbox{$\mathop{\text{HCFM}}\limits$}\xspace including $m_s$=``all", $m_c$=$2$, $\alpha$=$0.2$; and
$m_s$=6, $m_c$=$2$, $\alpha$=$0.8$, that is, with $\alpha$ varying and the other two fixed from the above two parameter
settings, respectively.
We present the results of these four parameter settings in Figure \ref{fig:session_length_all_2_0.8},
\ref{fig:session_length_all_2_0.2}, \ref{fig:session_length_6_2_0.8} and \ref{fig:session_length_6_2_0.2},
respectively.

Figure \ref{fig:session_length_all_2_0.8},
\ref{fig:session_length_all_2_0.2}, \ref{fig:session_length_6_2_0.8} and \ref{fig:session_length_6_2_0.2}
show the results in terms of HR@$1$, HR@$5$ and HR@$10$
for the five groups of search sessions.
The four figures show similar trends, that is,
top-20\%-40\% longest sessions achieve the best performance,
 %
and the top-20\% longest sessions achieve the worst.
This indicates that previous search terms from long time ago  may not represent
the information that the clinicians intend to search at the time of recommendation.
Furthermore, groups with top-20\% shortest sessions, top-20\%-40\% shortest sessions
and mid-40\%-60\% sessions achieve similar HR@$1$ and HR@$10$, and are all much
worse than those with top-20\%-40\% longest sessions.
This indicates that  too few search terms may not contain sufficient information to produce
accurate search term recommendations, while the top-20\%-40\% longest sessions,
that is, the sessions with length from 9 to 24 (Table~\ref{tbl:session_length}), can achieve accurate recommendations
with adequate search term information.

Moreover, \mbox{$\mathop{\text{HCFM}}\limits$}\xspace with $m_s$ as ``all" and $m_s$=$6$ show a similar trend.
Specifically, for the sessions with length less or equal to $6$
(i.e., top-20\% shortest sessions and top-20\%-40\% shortest sessions), \mbox{$\mathop{\text{HCFM}}\limits$}\xspace with the
two parameter settings  have identical results.
This is because when $m_s$ is larger than the number of all previous search terms in the current session,
we will use all the search terms to recommend the next search term.
In this case, $m_s$=$6$ will have the same results as $m_s$=``all".
For the sessions of mid-40\%-60\% lengths, given that their average length is 6.76 as shown in
Table \ref{fig:session_length},
the 6 most recent searches constitute 88.76\% of a search session on average.
Thus, using the 6 most recent searches in a session shows similar trends
as using all previous searches for these sessions.
For the top-20\%-40\% longest sessions, the average length is 15.04 as shown in Table \ref{fig:session_length},
and  the 6 most recent searches constitute 39.90\% of a search session on average.
Furthermore, comparing Figure~\ref{fig:session_length_all_2_0.8} and \ref{fig:session_length_6_2_0.8},
Figure \ref{fig:session_length_all_2_0.2} and \ref{fig:session_length_6_2_0.2}, we notice that
using previous 6 searches achieves significant improvement
compared to using all previous searches in term of HR@$5$ for top-20\%-40\% longest sessions.
This indicates that although using the most recent 6 searches discards a lot of information in a session,
the most recent, pertinent information still retains in these sessions and thus enables
superior recommendation performance of \mbox{$\mathop{\text{HCFM}}\limits$}\xspace.
For top-20\% longest sessions, the average length is 62.75 as shown in Table \ref{fig:session_length}.
For these sessions,  the most recent 6 searches on average constitute 9.56\% of a search session.
Therefore, using the most recent 6 searches for recommendation may lose too much information
and may have equally bad recommendation results as using all the previous searches, which, on the other hand,
may contain too much irrelevant information.

\subsection{Case Study}
\label{sec:results:casestudy}

\input{tables/case_study}

We present a case study based on the results in Section~\ref{sec:results:session_length}.
Specifically, we extract top-20\%-40\% longest sessions (as shown in Table \ref{fig:session_length})
and divide the sessions into two groups: top-5 hit group and top-5 miss group,
based on whether their testing search terms are correctly  recommended among top 5 or not.
It turns out that each group has 12 sessions.
We examine the testing search terms of each group and investigate the following two questions:
(1) is the testing search term frequently searched? and
(2) has the testing search term appeared earlier in a same session?
Table \ref{tbl:case-study-freq-terms} presents the testing search terms and
how many times they are searched on all patients
in our dataset (i.e., including training set and testing set).
Table~{\ref{tbl:case-study-freq-terms}}  also shows the percentage of search sessions that have a
corresponding testing search term.
Table \ref{tbl:case-study-freq-terms} shows that for
top-$5$ hit group, all the search terms in the testing set are frequently searched, and are searched for at least 200 times.
However, for top-$5$ miss group, only 3 search terms in the testing set are searched for over 200 times, and 7 search terms
are searched for less than 100 times.
This shows that if the testing search term is frequently searched, the model is more likely to generate accurate recommendations
for this search term.
On the other hand, if the testing search term is infrequently searched, it is hard for the model to generate accurate recommendations.
Moreover, in top-$5$ hit group, 75.00\% of the sessions have their testing search terms also appear earlier in the same sessions,
whereas in top-$5$ miss group, only 41.67\% of the sessions have the testing search terms in the earlier part of the same sessions.
This indicates that for each session, if the testing search term has appeared earlier in the same session,
\mbox{$\mathop{\text{HCFM}}\limits$}\xspace is much likely to generate accurate recommendations for this search term.

\section{Conclusion}
\label{sec:conclusion}

In this manuscript, we developed a model named
Hybrid Collaborative Filtering model using encounter information for search term
recommendation  (\mbox{$\mathop{\text{HCFM}}\limits$}\xspace) for clinicians.
The \mbox{$\mathop{\text{HCFM}}\limits$}\xspace model generates recommendations for the next search using two factors: previous search terms
and previous encounters.
Specifically, \mbox{$\mathop{\text{HCFM}}\limits$}\xspace recommends the terms that have high co-occurrence frequencies with the most recent ICD  codes
and are highly relevant to the most recent search terms for the next search.
We conduct comprehensive experiments on different datasets to compare the proposed model with the state-of-the-art
baseline methods.
The experimental results demonstrate that the proposed model can outperform all the
baseline methods with improvement of at least 4.56\%, 7.28\%, 7.15\% and
5.54\% over the second best results
 in terms of  HR@$3$, HR@$4$ , HR@$5$ and HR@$10$, respectively, on different datasets with different sparsities.

 %
We have also conducted comprehensive parameter studies to analyze the impact of the two factors.
The experimental results show that when using the most recent previous 5 to 9 search terms and
the most recent previous two encounters, with the weight of previous search term factor between 0.2 and 0.8,
\mbox{$\mathop{\text{HCFM}}\limits$}\xspace is most likely to generate accurate search term recommendations.
Furthermore,
we evaluated the \mbox{$\mathop{\text{HCFM}}\limits$}\xspace model on individual session groups divided by session length.
The experimental results show that the search sessions with length from 9 to 24 are
more likely to have better performance than a shorter or longer session.
Finally, we conducted a case study to better understand the performance of \mbox{$\mathop{\text{HCFM}}\limits$}\xspace, and concluded that
the \mbox{$\mathop{\text{HCFM}}\limits$}\xspace model tends to have better performance if the testing search term is  frequently searched
and has appeared in a same search session before.

\section*{Acknowledgements}

This project was made possible, in part, by support from the National
Science Foundation under Grant Number IIS-1855501 and
IIS-1827472, and from National Library of Medicine under Grant Number 1R01LM012605-01A1.
Any opinions, findings, and conclusions or recommendations
expressed in this material are those of the authors and do
not necessarily reflect the views of the funding agencies.

\section{Transition-Involved Patient-Term-Similarity-based CF Scoring}
\label{sec:app:baseline}

The \mbox{$\mathop{\text{TptCF}}\limits$}\xspace model generates recommendations for the next search term using two factors:
(1) patients' similarities and search terms' similarities (similarity-based scoring) and
(2) search term dynamic transitions (dynamics-based scoring).
Specifically, \mbox{$\mathop{\text{TptCF}}\limits$}\xspace calculates the patient-patient similarities and the term-term similarities using
all patients' search histories, and defines similar patients and similar search terms accordingly.
For each patient, the similar search terms of her last search term, and the search terms from her
similar patients will be used to calculate a recommendation score of each term.
This similarity-based recommendation score will be
combined with the dynamics-based score that is calculated based on term-term transition relations.
Particularly, the recommendation score of  term $s$ on patient $p$ is calculated as follows:

\begin{equation}
\label{eqn:TptCF}
    \begin{aligned}
\operatorname{Score}(t | \vec{T}(p))   & \   =  (1-\alpha) \cdot \operatorname{Score}_{\mathrm{DYN}}(t | \vec{T}(p)) \\
             					          & \   + \alpha \cdot \operatorname{Score}_{\mathrm{CF}}(t | \vec{T}(p)),
    \end{aligned}
\end{equation}
where
\begin{equation}
    \label{eqn:DYN}
\operatorname{Score}_{\mathrm{DYN}} (t | \vec{T}(p))=P\left(t | t_j\right),
\end{equation}
and
\vspace{-8pt}
\begin{equation}
    \label{eqn:TptCF}
    \begin{aligned}
\operatorname{Score}_{\mathrm{CF}} (t | \vec{T}(p))  &  =  \sum_{p^{\prime} \in \mathcal{S}_{p}(p)}\{\frac{\operatorname{sim}_{p}\left(p, p^{\prime}\right)}{\sum_{p^{\prime \prime} \in \mathcal{S}_{p}(p)} \operatorname{sim}_{p}\left(p, p^{\prime \prime}\right)}  \\
                   &  \times \sum_{t^{\prime} \in \mathcal{S}_{t}\left(t_j\right)} \frac{g\left(t^{\prime} \rightarrow t | p^{\prime}\right) \operatorname{sim}_{t}\left(t_j, t^{\prime}\right)}{\sum_{t^{\prime \prime} \in \mathcal{S}_{t}\left(t_j\right)} g\left(t^{\prime \prime} \rightarrow t | p^{\prime}\right)}\}
    \end{aligned}
\end{equation}
where $\operatorname{Score}_{\mathrm{DYN}}$ and $\operatorname{Score}_{\mathrm{CF}}$ are
dynamics-based recommendation score and similarity-based score, respectively;
$P\left(t | t_j\right)$ is the transition probability from term $t_j$ to term $t$;
$t_j$ is the last search term prior to time $t$;
$\vec{T}(p)$ denotes the sequence of patient $p$'s search terms;
$\mathcal{S}_{p}(p)$ denotes the set of patients similar to $p$;
$\mathcal{S}_{t}(t)$ denotes the set of search terms similar to $t$;
$g\left(t^{\prime} \rightarrow t | p^{\prime}\right)$ is the frequency of
transitions from term $t^{\prime}$ to term $t$ for patient $p^{\prime}$.
$\operatorname{sim}_{p}\left(p, p^{\prime}\right)$ is the similarity
between patient $p$ and $p^{\prime}$.
$\operatorname{sim}_{t}\left(t_j, t^{\prime}\right)$ is the similarity
between term $t_j$ and $t^{\prime}$.
The \mbox{$\mathop{\text{TptCF}}\limits$}\xspace model has three pre-defined parameters: $S_p$, $S_t$ and $\alpha$, among
which $S_t$, $S_p$ are similarity thresholds to define similar patients and similar terms,
respectively; $\alpha \in [0,1]$ is the weight of similarity-based scoring for recommendation.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi\

\end{document}


\section{INTRODUCTION}

Learning robot control policies in simulation \cite{Cutler2015EfficientRL} is typically safer, cheaper, and faster than learning in the real world, but it also introduces a \textit{reality gap} \cite{Kober2012ReinforcementLI} between the training environment (the simulator) and the deployment environment (the real world). {\em Sim-to-real} algorithms,  which focus on overcoming this gap,  have recently received a great deal of attention.

This paper focuses on a class of solutions referred to as grounding algorithms \cite{AAMAS13-Farchy} which use a small amount of real world data to improve (or \textit{ground}) the simulator.
We can assume that the ungrounded simulator at least approximates the correct real world dynamics, so many grounding approaches \cite{AAMAS13-Farchy, Tan_2018, nas_golemo, AAAI17-Hanna} learn parameters to correct for these differences.
Since it may not always be feasible to modify the internal parameters of a simulator, we focus our attention on grounding algorithms that treat the simulator as a black-box, such as the Grounded Action Transformation (\textsc{gat}) algorithm \cite{AAAI17-Hanna}.

The dynamics of robots are often best modeled as stochastic processes due to effects like friction, gear backlash, uneven terrain, and other sources of noise in the environment; however many earlier sim-to-real approaches do not explicitly account for stochasticity in the real world. These approaches learn a single value for the correction terms whereas learning a distribution would more accurately reflect the real world.

We hypothesize that accounting for the stochasticity in the real world improves sim-to-real transfer of policies learned in simulation. In this work, we analyze the effects of target environment stochasticity in the {\em sim-to-real} problem. We show several domains where \textsc{gat} fails to adequately ground the simulator, and we propose a new algorithm, {\em Stochastic Grounded Action Transformation} (\textsc{sgat}), that handles this issue gracefully by learning the stochasticity in the environment.
We first show on the \textit{Cliff Walking} domain that using \textsc{sgat} achieves better mean returns than \textsc{gat} on the target environment. We further conduct sim-to-``real" transfer experiments on OpenAI gym MuJoCo environments \textit{InvertedPendulum} and \textit{HalfCheetah}, using a modified version of the environment as a surrogate for the real world.

To test the efficacy of \textsc{sgat} in handling real world stochasticity, we set up a sim-to-real experiment where a humanoid \textsc{nao} robot learns to walk on uneven ground. Confirming our hypothesis, we indeed found that policies trained with \textsc{sgat}, unlike those learned by \textsc{gat}, learned to cope with environment stochasticity, leading to better real world performance. Using \textsc{sgat}, the \textsc{nao} completed the course 9 out of 10 times (Fig. \ref{fig:lumpy_field_setup}), whereas with \textsc{gat}, it fell down every time.

\begin{figure*}[!tb]
    \centering
    \includegraphics[scale=1.1, trim={0 0.3cm 0 0.3cm} ,clip]{robotwalk.png}
    \caption{Experiment setup showing a  robot walking on the uneven ground. The \textsc{nao} begins walking 40 cms behind the center of the circle and walks 300 cms towards the white penalty cross. This image shows a successful walk executed by the robot at 2 sec intervals, learned using the proposed \textsc{sgat} algorithm.}
    \label{fig:lumpy_field_setup}
\end{figure*}

\section{BACKGROUND}

In this section, we review existing literature on state-of-the-art sim-to-real algorithms that are relevant to our approach, and we formalize the problem addressed in this paper.

\subsection{Related Work}

In the past, many sim-to-real techniques \cite{Tan_2018, DBLP:journals/corr/abs-1710-06537, Jakobi1995, miglino1995, rajeswaran2016epopt} have achieved policy learning for the real world using a simulator.
We can roughly divide these techniques into black-box techniques (which do not modify the simulator internally) and white-box techniques (which do). Another useful distinction is robustness methods vs. grounding methods.

Robustness methods learn a policy in simulation that is robust to variations in the environment. Dynamics Randomization (\textsc{DR}) is a relatively simple approach that has been shown to be successful at robot control tasks.
Peng et al. \cite{DBLP:journals/corr/abs-1710-06537} introduce a white-box \textsc{DR} algorithm that can learn robust RL control policies for robotic manipulation tasks.
DR approaches where the effect of actuators (actions) or sensor information (observations) are perturbed with an ``envelope" of noise can be considered black-box approaches \cite{Jakobi1995, miglino1995}. Other DR methods where the internal simulation parameters are perturbed fall under white-box approaches \cite{Tan_2018, DBLP:journals/corr/abs-1710-06537, Molchanov_2019, Yu_2017}. Robust Adversarial Reinforcement Learning (\textsc{rarl}) \cite{pinto2017robust} randomizes the training environment using an adversarial agent. While DR methods randomly sample simulation parameters before a trajectory is generated, the simulator dynamics may still be deterministic. Our approach models real world stochasticity during individual state transitions.

Orthogonal to robustness methods are approaches that focus on grounding a simulator to behave like a target real world environment.
Unlike DR methods that learn a more general robust policy to perform well in any environment, grounding methods optimize performance in that specific target environment. Grounded Simulation Learning (\textsc{gsl}) \cite{AAMAS13-Farchy} is a learning framework in which data from the real world is used to modify (ground) the simulator. After this grounding step, transitions in the grounded simulator look similar to transitions in the real world environment, hence minimizing the ``reality gap". In \textsc{gsl}, all of the policy learning happens in the simulator; the robot is used only to evaluate policies and to collect transition data. Farchy et al. demonstrate its success at transferring a walk policy from a simulated \textsc{nao} in the SimSpark simulator to a real SoftBank \textsc{nao} and achieved 26.7\% faster walk than baseline methods  \cite{AAMAS13-Farchy}.

Grounded Action Transformation (\textsc{gat}) \cite{AAAI17-Hanna} is a recent black-box algorithm based on the \textsc{gsl} framework.
\textsc{gat} was experimentally shown to be effective at learning a fast walk policy by grounding a high fidelity simulator, which when deployed in the real world achieved the fastest known walk on the \textsc{nao} robot: over 43\% faster than the state-of-the-art hand-coded gait upon which it was based.

\subsection{Preliminaries}

We consider the real world (real) and simulation (sim) domains to be two different Markov Decision Processes (\textsc{mdp}s) \cite{Suttonbook}, $\mathcal{M}_{real}$ and $\mathcal{M}_{sim}$ respectively. An \textsc{mdp} comprises a set of states, $\mathcal{S}$, a set of actions, $\mathcal{A}$, the transition dynamics associated with those actions, $\mathcal{T}$, and a reward function, $R$. At each time step, $t$, an agent observes the current state, $s_t \in \mathcal{S}$, and chooses an action, $a_t \in \mathcal{A}$, sampled from its policy, $a_t \sim \pi(\cdot|s_t)$.
For a given state and action, there is a distribution of possible next states, from which the environment samples a state---$s_{t+1} \sim \mathcal{T}(\cdot|s_t,a_t)$.
The reward, in this work, is a function of the action and the next state, $r_{t+1} = R(a_t, s_{t+1})$.

In our formulation of the sim-to-real problem, only the transition dynamics $\mathcal{T}_{sim}$ and $\mathcal{T}_{real}$ differ between the environments. Furthermore, we only consider deterministic simulators in this work, but the real environment may have stochastic transitions. The reinforcement learning (RL) objective is to find a policy that maximizes the expected sum of rewards on the real domain, $\mathbb{E}[\sum_{t=0}^T R(a_t, s_{t+1})]$.

\subsection{Grounded Action Transformation (\textsc{gat})}

\begin{figure}[!b]
    \centering
    \begin{tikzpicture}
    \node[process, text width=20mm] (pol) {Policy};
    \filldraw[fill=black!10, draw=black!50, dashed] (2.3,-2.4) rectangle (5.2,.7);
    \node[process, right=16mm of pol, fill=blue!20] (fwd) {$f_{real}(s_t,a_t)$};
    \node[process, below=8mm of fwd, fill=red!20] (inv) {$f^{-1}_{sim}(s_t,\hat{s}_{t+1})$};
    \node[process, below=8mm of pol, text width=20mm] (env) {Simulated Environment};
    \draw[arrow] (pol) -- node [above, xshift=-2mm] {$a_t$} (fwd);
    \draw[arrow] (fwd) -- node [right] {$\hat{s}_{t+1}$} (inv);
    \draw[arrow] (inv) -- node [above] {$\hat{a}_t$} (env);
    \draw[arrow] (env) -- node[right] {$s_{t+1}$} (pol);
    \end{tikzpicture}
    \caption{\textsc{gat} Diagram\cite{AAAI17-Hanna}}
    \label{fig:gat_diagram}
\end{figure}

\textsc{gat} follows the \textsc{gsl} framework, alternating between a \textit{grounding step} and a \textit{policy improvement step}. During the grounding step, the policy remains unchanged and is deployed on both sim and real environments to collect state transition data. The \textsc{gat} algorithm trains two neural networks---a deterministic forward model $\hat{s}_{t+1}=f_{real}(s_t,a_t)$ and an inverse dynamics model $\hat{a}_t=f^{-1}_{sim}(s_t,{s}_{t+1})$  of the real and sim environments respectively, using supervised learning. The learned forward and inverse models are composed to form the action transformer function $g(s_t,a_t)=f^{-1}_{sim}(s_t,f_{real}(s_t,a_t))$ that grounds the simulator. The inputs to the action transformer are the current state, $s_t$, and the action chosen by the agent, $a_t \sim \pi(\cdot|s_t)$. The output is a transformed action, $\hat{a}_t$, that when taken in the simulator, produces a transition similar to the expected transition when executing $a_t$ in the real world. The combination of the action transformer and the original simulator is known as the grounded simulator. The grounding step is followed by a policy improvement step, in which the parameters of the action transformer are fixed. Thus, the policy improvement step is a standard reinforcement learning problem in the grounded simulator. These two steps are repeated until the optimized policy performs well on the real environment. The block diagram depicting the \textsc{gat} framework is shown in Fig. \ref{fig:gat_diagram}. While \textsc{gat} works well on fairly deterministic environments, as was shown by Hanna and Stone \cite{AAAI17-Hanna}, in our experimentation, we find that policies learned using \textsc{gat} perform poorly when transferring to highly stochastic environments.

\subsection{Transferring to a Stochastic Real World}

When the real world domain is deterministic, learning a deterministic forward model, as \textsc{gat} does, works well.
We can demonstrate using a toy example. Consider the environments shown in Figs. \ref{fig:toy1a} and \ref{fig:toy1b}. The agent starts in the initial state, $s_0$, and chooses between $a_1$ and $a_2$. In the simulator, $a_1$ leads to the higher reward of +1, but the transitions are flipped in reality. Hence, the action $a_2$ is the optimal action. The \textsc{gat} action transformer learns to transform $a_1$ into $a_2$ and $a_2$ into $a_1$. Thus from the agent's perspective the grounded simulator behaves like the real world.

\begin{figure}[!tb]
\centering
\begin{subfigure}{.2\textwidth}
    \centering
    \begin{tikzpicture}[auto,>=latex]
        \node[state] (s0) {$s_0$};
        \node[action,above right=2mm and 4mm of s0,label={$a_1$}] (a1) {};
        \node[action,below right=2mm and 4mm of s0,label={$a_2$}] (a2) {};
        \draw[arrow] (s0) -- (a1);
        \draw[arrow] (s0) -- (a2);
        \node[state, right=8mm of a1, label=right:{+1}] (s1) {$s_1$};
        \node[state, right=8mm of a2, label=right:{-1}] (s2) {$s_2$};
        \draw[arrow] (a1) -- (s1);
        \draw[arrow] (a2) -- (s2);
    \end{tikzpicture}
    \caption{Simulator}
    \label{fig:toy1a}
\end{subfigure}
\begin{subfigure}{.2\textwidth}
    \centering
    \begin{tikzpicture}[auto,>=latex]
        \node[state] (s0) {$s_0$};
        \node[action,above right=2mm and 4mm of s0,label={$a_1$}] (a1) {};
        \node[action,below right=2mm and 4mm of s0,label={$a_2$}] (a2) {};
        \draw[arrow] (s0) -- (a1);
        \draw[arrow] (s0) -- (a2);
        \node[state, right=8mm of a1, label=right:{+1}] (s1) {$s_1$};
        \node[state, right=8mm of a2, label=right:{-1}] (s2) {$s_2$};
        \draw[arrow] (a1) -- (s2);
        \draw[arrow] (a2) -- (s1);
    \end{tikzpicture}
    \caption{Real World}
    \label{fig:toy1b}
\end{subfigure}
\begin{subfigure}{.4\textwidth}
    \centering
    \begin{tikzpicture}[auto,>=latex]
        \node[state] (s0) {$s_0$};
        \filldraw[fill=black!10, draw=black!50, dashed] (0.5,-1) rectangle (3,1);
        \node[action,above right=2mm and 4mm of s0,label={$a_1$}] (a1) {};
        \node[action,below right=2mm and 4mm of s0,label={$a_2$}] (a2) {};
        \draw[arrow] (s0) -- (a1);
        \draw[arrow] (s0) -- (a2);
        \node[state, dashed, right=8mm of a1] (s1i) {$\hat{s}_1$};
        \node[state, dashed, right=8mm of a2] (s2i) {$\hat{s}_2$};
        \draw[arrow, color=blue] (a1) -- (s2i);
        \draw[arrow, color=blue] (a2) -- (s1i);
        \node[action, right=8mm of s1i,label={$\hat{a}_1$}] (a1i) {};
        \node[action, right=8mm of s2i,label={$\hat{a}_2$}] (a2i) {};
        \draw[arrow, color={red}] (s1i) -- (a1i);
        \draw[arrow, color={red}] (s2i) -- (a2i);
        \node[state, right=8mm of a1i, label=right:{+1}] (s1) {$s_1$};
        \node[state, right=8mm of a2i, label=right:{-1}] (s2) {$s_2$};
        \draw[arrow] (a1i) -- (s1);
        \draw[arrow] (a2i) -- (s2);
    \end{tikzpicture}
    \caption{\textsc{gat} Grounded Simulator}
    \label{fig:toy1c}
\end{subfigure}
\caption{For a deterministic domain, the forward model (blue) and the inverse model (red) can make the grounded simulator behave like the real world. \textsc{gat} works well here.}
\label{fig:toy1}
\end{figure}

\begin{figure}[!tb]
\centering
\begin{subfigure}{.2\textwidth}
    \centering
    \begin{tikzpicture}[auto,>=latex]
        \node[state] (s0) {$s_0$};
        \node[action, right=4mm of s0,label={$a_2$}] (a2) {};
        \node[action, above=8mm of a2,label={$a_1$}] (a1) {};
        \node[action, below=8mm of a2,label={$a_3$}] (a3) {};
        \draw[arrow] (s0) -- (a1);
        \draw[arrow] (s0) -- (a2);
        \draw[arrow] (s0) -- (a3);
        \node[state, right=8mm of a1, label=right:{+1}] (s1) {$s_1$};
        \node[state, right=8mm of a2, label=right:{-1}] (s2) {$s_2$};
        \node[state, right=8mm of a3, label=right:{+10}] (s3) {$s_3$};
        \draw[arrow] (a1) -- (s1);
        \draw[arrow] (a2) -- (s2);
        \draw[arrow] (a3) -- (s3);
    \end{tikzpicture}
    \caption{Simulator}
    \label{fig:toy2a}
\end{subfigure}
\begin{subfigure}{.2\textwidth}
    \centering
    \begin{tikzpicture}[auto,>=latex]
        \node[state] (s0) {$s_0$};
        \node[action, right=4mm of s0,label={$a_2$}] (a2) {};
        \node[action, above=8mm of a2,label={$a_1$}] (a1) {};
        \node[action, below=8mm of a2,label={$a_3$}] (a3) {};
        \draw[arrow] (s0) -- (a1);
        \draw[arrow] (s0) -- (a2);
        \draw[arrow] (s0) -- (a3);
        \node[state, right=10mm of a1, label=right:{+1}] (s1) {$s_1$};
        \node[state, right=10mm of a2, label=right:{-1}] (s2) {$s_2$};
        \node[state, right=10mm of a3, label=right:{+10}] (s3) {$s_3$};
        \coordinate [left=2.5mm of s3, label=below:{20\%}] (c3);
        \draw[arrow] (a1) -- (s1);
        \draw[arrow] (a2) -- node [above, xshift=2mm] {80\%} (s2);
        \draw[arrow] (a2) -- (c3) -- (s3);
        \draw[arrow] (a3) -- (s2);
    \end{tikzpicture}
    \caption{Real World}
    \label{fig:toy2b}
\end{subfigure}
\begin{subfigure}{.4\textwidth}
    \centering
    \begin{tikzpicture}[auto,>=latex]
        \node[state] (s0) {$s_0$};
        \filldraw[fill=black!10, draw=black!50, dashed] (0.5,-1.45) rectangle (3.1,1.45);
        \node[action, right=4mm of s0,label={$a_2$}] (a2) {};
        \node[action, above=8mm of a2,label={$a_1$}] (a1) {};
        \node[action, below=8mm of a2,label={$a_3$}] (a3) {};
        \draw[arrow] (s0) -- (a1);
        \draw[arrow] (s0) -- (a2);
        \draw[arrow] (s0) -- (a3);
        \node[state, dashed, right=8mm of a1] (s1i) {$\hat{s}_1$};
        \node[state, dashed, right=8mm of a2] (s2i) {$\hat{s}_2$};
        \node[state, dashed, right=8mm of a3] (s3i) {$\hat{s}_3$};
        \draw[arrow, color=blue] (a1) -- (s1i);
        \draw[arrow, color=blue] (a2) -- (s2i);
        \draw[arrow, color=blue] (a3) -- (s2i);
        \node[action, right=8mm of s1i,label={$\hat{a}_1$}] (a1i) {};
        \node[action, right=8mm of s2i,label={$\hat{a}_2$}] (a2i) {};
        \node[action, right=8mm of s3i,label={$\hat{a}_3$}] (a3i) {};
        \draw[arrow, color=red] (s1i) -- (a1i);
        \draw[arrow, color=red] (s2i) -- (a2i);
        \draw[arrow, color=red] (s3i) -- (a3i);
        \node[state, right=8mm of a1i, label=right:{+1}] (s1) {$s_1$};
        \node[state, right=8mm of a2i, label=right:{-1}] (s2) {$s_2$};
        \node[state, right=8mm of a3i, label=right:{+10}] (s3) {$s_3$};
        \draw[arrow] (a1i) -- (s1);
        \draw[arrow] (a2i) -- (s2);
        \draw[arrow] (a3i) -- (s3);
    \end{tikzpicture}
    \caption{\textsc{gat} Grounded Simulator}
    \label{fig:toy2c}
\end{subfigure}
\caption{When the real world has stochastic transitions, the \textsc{gat} forward model (blue) only captures the most likely next state. \textsc{gat} may fail here, since the grounded simulator behaves very differently from the real environment.}
\label{fig:toy2}
\end{figure}

However, when we add stochastic transitions, the two diagrams do not match. In Fig. \ref{fig:toy2}, the optimal action in the simulator is $a_3$, and in the real world, it is $a_2$; however, in the grounded simulator, it is $a_1$.  Since \textsc{gat}'s forward model is deterministic, it predicts only the most likely next state, but other, less likely transitions are also important when computing an action's value.

\section{STOCHASTIC GROUNDED ACTION TRANSFORMATION (\textsc{SGAT})}

\begin{figure}[!tb]
    \centering
    \begin{tikzpicture}[auto,>=latex]
        \node[state] (s0) {$s_0$};
        \filldraw[fill=black!10, draw=black!50, dashed] (0.5,-1.45) rectangle (3.1,1.45);
        \node[action, right=4mm of s0,label={$a_2$}] (a2) {};
        \node[action, above=8mm of a2,label={$a_1$}] (a1) {};
        \node[action, below=8mm of a2,label={$a_3$}] (a3) {};
        \draw[arrow] (s0) -- (a1);
        \draw[arrow] (s0) -- (a2);
        \draw[arrow] (s0) -- (a3);
        \node[state, dashed, right=8mm of a1] (s1i) {$\hat{s}_1$};
        \node[state, dashed, right=8mm of a2] (s2i) {$\hat{s}_2$};
        \node[state, dashed, right=8mm of a3] (s3i) {$\hat{s}_3$};
        \coordinate [left=2.5mm of s3i, label=below:{20\%}] (c3);
        \draw[arrow, color=blue] (a1) -- (s1i);
        \draw[arrow, color=blue] (a2) -- node [above, xshift=2mm, color=black] {80\%} (s2i);
        \draw[arrow, color=blue] (a2) -- (c3) -- (s3i);
        \draw[arrow, color=blue] (a3) -- (s2i);
        \node[action, right=8mm of s1i,label={$\hat{a}_1$}] (a1i) {};
        \node[action, right=8mm of s2i,label={$\hat{a}_2$}] (a2i) {};
        \node[action, right=8mm of s3i,label={$\hat{a}_3$}] (a3i) {};
        \draw[arrow, color=red] (s1i) -- (a1i);
        \draw[arrow, color=red] (s2i) -- (a2i);
        \draw[arrow, color=red] (s3i) -- (a3i);
        \node[state, right=8mm of a1i, label=right:{+1}] (s1) {$s_1$};
        \node[state, right=8mm of a2i, label=right:{-1}] (s2) {$s_2$};
        \node[state, right=8mm of a3i, label=right:{+10}] (s3) {$s_3$};
        \draw[arrow] (a1i) -- (s1);
        \draw[arrow] (a2i) -- (s2);
        \draw[arrow] (a3i) -- (s3);
    \end{tikzpicture}
    \caption{In the \textsc{sgat} Grounded Simulator, the transitions match the real environment (Fig. \ref{fig:toy2b}).}
    \label{fig:toy3}
\end{figure}

To address real world stochasticity, we introduce Stochastic Grounded Action Transformation (\textsc{sgat}), which learns a stochastic model of the forward dynamics. In other words, the learned forward model, $f_{real}$, predicts a distribution over next states (from which we sample) rather than the most likely next state. The sampling operation within the action transformer makes the overall process stochastic.
Fig. \ref{fig:toy3} illustrates the simulator from Fig. \ref{fig:toy2} grounded using \textsc{sgat}. Since the forward model accounts for stochasticity in the real world, the actions in the grounded simulator have the same effect as in the real world.

In continuous state and action domains, we model the next state as a multivariate Gaussian distribution and train the forward model using negative log likelihood (\textsc{nll}) loss $\mathcal{L} = -\log p(s_{t+1}|s_t,a_t)$.
Similar to \textsc{gat}, we use a neural network function approximator with two fully connected hidden layers of 64 neurons to represent the forward and inverse models, but unlike \textsc{gat}, the forward model in \textsc{sgat} outputs the parameters of a Gaussian distribution from which we sample the predicted next state.\footnote{A Mixture Density Network might be more suitable when the environment's transition dynamics exhibit multimodal behavior.}
In our implementation, the final dense layer outputs the mean, $\mu$, and the log standard deviation, $log(\sigma)$, for each element of the state vector.
The complete algorithm of \textsc{sgat} is shown in Algorithm \ref{alg:SGAT}.

\begin{algorithm}[!tb]
\caption{Stochastic Grounded Action Transformation}
\label{alg:SGAT}
\textbf{Input}: initial parameters $\theta_0$, $\phi$, and $\psi$ for the policy $\pi_\theta$, forward dynamics model $f_{{\phi}_{real}}$ and inverse dynamics model $f^{-1}_{{\psi}_{sim}}$; policy improvement method, \texttt{optimize}
\begin{algorithmic}[1]
\WHILE{policy $\pi_\theta$ improves on real}
 \STATE Collect real world trajectories \\$\tau_{real} \leftarrow \{((s_0,a_0),s_1), ((s_1, a_1), s_2),...\}_{real}$
 \STATE Train forward dynamics function $f_\phi$ with $\tau_{real}$, using the \textsc{nll} loss $\mathcal{L} = -log$ $p(s_{t+1}|s_t,a_t)$.
 \STATE Collect simulated trajectories \\$\tau_{sim} \leftarrow \{((s_0,a_0),s_1), ((s_1, a_1), s_2),...\}_{sim}$
 \STATE Train inverse dynamics function $f^{-1}_\phi$ with $\tau_{sim}$, using mean squared error loss.
 \STATE Update $\pi_\theta$ on the grounded simulator using \texttt{optimize} and the reward from the grounded simulator
\ENDWHILE
\end{algorithmic}
\end{algorithm}

\section{EXPERIMENTS}

This section reports on an empirical study of transfer from simulation with \textsc{sgat} compared to \textsc{gat}.
We begin with a toy RL domain and progress to sim-to-real transfer of a bipedal walking controller for a \textsc{nao} robot.
Our empirical results show the benefit of modelling stochasticity when grounding a simulator for transfer to a stochastic real world environment.

\subsection{Cliff Walking}
\label{cliffwalkingdomainsection}

We verify the benefit of \textsc{sgat} using a classical reinforcement learning domain, the Cliff Walking grid world shown in Fig. \ref{fig:gridworld_diagram}. In this domain, an agent must navigate around a cliff to reach a goal. The episode terminates when it either reaches the goal (reward of $+100$) or falls into the cliff (reward of $-10$). There is also a small time penalty ($-0.1$ per time step), so the agent is incentivized to find the shortest path. There is no discounting, so the agent's objective is to maximize the sum of rewards over an episode. The agent can move up, down, left, or right. If it tries to move into a wall, the action has no effect. In our version of the problem, we assume we have a deterministic simulator, but in the ``real" environment, there is a small chance at every time step that the agent moves in a random direction instead of the direction it chose.

Fig. \ref{fig:grid_plot} shows \textsc{gat} and \textsc{sgat} evaluated for different values of the environment noise parameter. Both the grounding steps and policy improvement steps (using policy iteration \cite{Suttonbook}) are repeated until convergence for both algorithms. To evaluate the resulting policy, we estimate the expected return from averaging 10,000 episodes. At a value of zero, the ``real" environment is completely deterministic. At a value of one, every transition is random. Thus, at both of these endpoints, there is no distinction between the expected return gained by the two algorithms.

For every intermediate value, \textsc{sgat} outperforms \textsc{gat}.
The policy trained using \textsc{gat} is unaware of the stochastic transitions, so it always takes the shortest and most dangerous path. Meanwhile the \textsc{sgat} agent learns as if it were training directly on the real environment in the presence of stochasticity.

\begin{figure}[!tb]
    \centering
    \includegraphics[width=\columnwidth]{gridsimple.png}
    \caption{The agent starts in the bottom left and must reach the goal in the bottom right. Stepping into the red region penalizes the robot and ends the episode. The purple path is the most direct, but the blue path is safer when the transitions are stochastic.}
    \label{fig:gridworld_diagram}
\end{figure}

\begin{figure}[tb]
    \centering
    \includegraphics[width=\columnwidth]{gridplot.png}
    \caption{The y-axis is the average performance of a policy evaluated on the ``real" domain. The x-axis is the chance at each time step for the transition to be random. \textsc{sgat} outperforms \textsc{gat} for any noise value. Error bars not shown since standard error is smaller than 1 pixel.}
    \label{fig:grid_plot}
\end{figure}

\begin{figure}[tb]
    \centering
    \includegraphics[width=\columnwidth]{sgatfig8.png}
    \caption{Sim-to-NoisyReal experiment on \textbf{InvertedPendulum}. The ``real" pendulum is 10 times heavier than the sim pendulum and has added Gaussian noise of different values. Error bars show standard error over ten independent training runs. Algorithms with striped bars used no real world data during training. \textsc{sgat} performs comparatively better in noisier target environments.}
    \label{fig:invp_s2r}
\end{figure}

\subsection {MuJoCo Domains}

\begin{figure*}[!tb]
\begin{subfigure}{\columnwidth}
    \includegraphics[width=\columnwidth]{sgatfig9a.png}
    \caption{Sim-to-NoisySim}
    \label{fig:hc_s2s}
\end{subfigure}
\begin{subfigure}{\columnwidth}
\includegraphics[width=\columnwidth]{sgatfig9b.png}
    \caption{Sim-to-NoisyReal}
    \label{fig:hc_s2r}
\end{subfigure}
\caption{Sim-to-NoisySim and Sim-to-NoisyReal experiments on \textbf{HalfCheetah}. In the NoisyReal environment, the ``real" HalfCheetah's mass is 43\% greater than the sim HalfCheetah. Error bars show show standard error over ten independent training runs. Algorithms with striped bars used no real world data during training. When the ``real" environment is highly stochastic, \textsc{sgat} performs better than \textsc{gat}. Meanwhile, \textsc{ane} does poorly on less noisy scenarios.
}
\label{hc_exps}
\end{figure*}

Having shown the efficacy of \textsc{sgat} in a tabular domain, we now evaluate its performance in continuous control domains that are closer to real world robotics settings.
We perform experiments on the OpenAI Gym MuJoCo environments to compare the effectiveness of \textsc{sgat} and \textsc{gat} when there is added noise in the target domain. We consider the case with just added noise and the case with both noise and domain mismatch between the source and target environments. We call the former Sim-to-NoisySim and the latter Sim-to-NoisyReal.
We chose the \textit{InvertedPendulum} and \textit{HalfCheetah} domains to test \textsc{sgat} in environments with both low and high dimensional state and action spaces. For policy improvement, we use an implementation of Trust Region Policy Optimization (\textsc{trpo}) \cite{TRPO}, from the stable-baselines repository \cite{stable-baselines}, using their hyperparameters for the respective domains.

We also compare against the action-noise-envelope (\textsc{ane}) approach \cite{Jakobi1995}; however, it is not a perfectly fair comparison in the sense that robustness approaches are sensitive to user-defined hyperparameters that predict the variation in the environment---in this case, the magnitude of the added noise. In grounding approaches, these parameters can be automatically learned from real world data. In these experiments, we chose the ANE hyperparameters (e.g., noise value) with a coarse grid search.

We simulate stochasticity in the target domains by adding Gaussian noise with different standard deviation values to the actions input into the environment.
We omit the results of Sim-to-NoisySim experiments for \textit{InvertedPendulum} because both algorithms performed well on the transfer task. Fig. \ref{fig:invp_s2r} shows the performance on the ``real" environment of policies trained four ways---naively on the ungrounded simulator, with \textsc{sgat}, with \textsc{gat}, and with \textsc{ane}.
In this Sim-to-NoisyReal experiment, \textsc{sgat} performs much better than \textsc{gat} when the stochasticity in the target domain increases.
Fig. \ref{hc_exps} shows the same experiment on \textit{HalfCheetah}, both with and without domain mismatch. Both these environments have an action space of $[-1, 1]$.

The red dashed lines show the performance of a policy trained directly on the ``real" environment until convergence, approximately the best possible performance.
The axes are scaled respective to this line.
The error bars show the standard error across 10 trials with different initialization weights. As the stochasticity increases, \textsc{sgat} policies perform better than those learned using \textsc{gat}. Meanwhile, \textsc{ane} does well only for particular noise values, depending on its training hyperparameters.

\subsection{\textsc{Nao} Robot Experiments}

Until this point in our analysis, we have used a modified version of the simulator in place of the “real” world so as to isolate the effect of stochasticity (as opposed to domain mismatch).  However, the true objective of this research is to enable transfer to real robots, which may exhibit very different noise profiles than the simulated environments. Thus, in this section, we validate \textsc{sgat} on a real humanoid robot learning to walk on uneven terrain.

Our robot experiments were conducted on the SoftBank \textsc{nao} bipedal robot, used in the RoboCup robot soccer competitions. For learning in simulation, we use the SimSpark physics simulator, used in the RoboCup 3D Simulation league. We compared \textsc{gat} and \textsc{sgat} by independently learning control policies using these algorithms to walk on uneven terrain, as shown in Fig. \ref{fig:lumpy_field_setup}. To create uneven ground in the lab, we placed foam packing material under the turf of the robot soccer field. On this uneven ground, the walking dynamics become more random, since the forces acting on the foot are slightly different every time the robot takes a step. For the walk policy, we use the rUNSWift walk engine and the initial published parameter set $\theta_0$ \cite{walkengine2014}. This initial unoptimized policy achieves a speed of $14.66 \pm 1.65$ cm/s on the uneven terrain. The walk engine contains 16 parameters that we optimize on the grounded simulator using Covariance Matrix Adaptation - Evolutionary Strategy (\textsc{cma-es}), with a population size of 150. Each trajectory lasts for 7.5 seconds on the simulator or terminates when the robot falls down. The reward function is a sum of the forward velocity of the robot (in cm/s) and an early termination penalty of -15 if the robot falls down. Thus, \textsc{cma-es} optimizes the policy for walks that are faster and more stable in the grounded simulator.

On flat ground, both methods produced very similar policies, but on the uneven ground, the policy learned using \textsc{sgat} was more successful than a policy learned using \textsc{gat}.
We performed ten trial runs of the best policy learned using \textsc{sgat} and \textsc{gat} after each grounding step, and the average speed of the robot on the uneven terrain is shown in Table \ref{table:nao}.
The policy learned using \textsc{sgat} takes shorter steps and stays upright, thereby maintaining its balance on the uneven terrain, whereas the policy produced using \textsc{gat} learned to lean forward and walk faster, but fell down more often due to the uneven terrain. This result is best visualized in the supplementary video attached with this paper. Both algorithms produce policies that improve the walking speed across grounding steps. The \textsc{gat} policy after the second grounding step always falls over, whereas the \textsc{sgat} policy was more stable and finished the course 9 out of 10 times.

\begin{table}[!tb]
\caption{Speed and stability of \textsc{nao} robot walking on uneven ground. The initial policy $\theta_0$ walks at $14.66 \pm 1.65$ cm/s and always falls down. Both \textsc{sgat} and \textsc{gat} find policies that are faster, but \textsc{sgat} policies are more stable than policies learned using \textsc{gat}.}
\centering
\label{table:nao}
\begin{tabular}{| c | c | c | c | c |}
    \hline
    & \multicolumn{2}{|c|}{Grounding Step 1} & \multicolumn{2}{|c|}{Grounding Step 2} \\
    \hline
    & Speed (cm/s) & Falls & Speed (cm/s) & Falls \\
    \hline
    \textsc{gat} & $15.7 \pm 2.98$ & 6/10 & $18.5 \pm 3.63$ & 10/10 \\
    \hline
    \textsc{sgat} & $16.9 \pm 0.678$ & \textbf{0/10} & $18.0 \pm 2.15$ & \textbf{1/10} \\
    \hline
\end{tabular}
\end{table}

\section{DISCUSSION AND LIMITATIONS}

In our experiments, we observe that stochastic transitions can have a detrimental effect on how \textsc{gat} grounds the simulator; however both algorithms perform similarly on deterministic environments. In real world scenarios, we cannot know how stochastic an environment is before testing. This fact suggests that we should default to \textsc{sgat}.

In this work, we have only considered deterministic simulators, but simulators may have stochastic transitions as well, especially if the simulator was designed to anticipate process noise. However, when using an action transformer grounding approach, stochastic simulators make the learning problem more difficult. We can no longer sample from the distribution provided by the forward model. Instead, the inverse model must take in a distribution over states and output a distribution over actions.

\section{CONCLUSION AND FUTURE WORK}

In this work, we introduced  \textit{Stochastic Grounded Action Transformation} (\textsc{sgat}), a sim-to-real algorithm, which gracefully handles policy learning in simulation when the target domain is stochastic.

First, in \textit{Cliff Walking}, we empirically showed that by accounting for target domain stochasticity, \textsc{sgat} policies can produce much better average returns on the real environment than \textsc{gat} policies. Then, we evaluated performance on two MuJoCo environments and showed that \textsc{sgat} outperforms other algorithms when the target environment is highly stochastic. To verify the algorithm in a real world transfer setting, we trained a \textsc{nao} humanoid robot to walk over uneven terrain. The empirical results show that \textsc{sgat} policies are more robust to process noise and walk more robustly over the uneven terrain compared to \textsc{GAT} and hand-coded polices.

Though we have focused on GAT, similar distributional approaches are directly applicable to other black-box sim-to-real techniques \cite{rgat, higuera17}, and possibly to some white-box techniques as well. Many sim-to-real algorithms are designed to find correct environment parameters, but the most accurate model might be one in which these parameters are sampled from a distribution. Investigating this scenario is a promising direction for future work.

\section{ACKNOWLEDGMENT}
\small{This work has taken place in the Learning Agents Research Group (LARG) at UT Austin.  LARG research is supported in part by NSF (CPS-1739964, IIS-1724157, NRI-1925082), ONR (N00014-18-2243), FLI (RFP2-000), ARL, DARPA, Lockheed Martin, GM, and Bosch.  Peter Stone serves as the Executive Director of Sony AI America and receives financial compensation for this work.  The terms of this arrangement have been reviewed and approved by the University of Texas at Austin in accordance with its policy on objectivity in research.}

\bibliographystyle{IEEEtran}