\section{Introduction}
Reinforcement Learning (RL) is a powerful tool to learn control policies from experience where the learning agent does not have explicit information about the task but is only allowed to interact with the environment and observe rewards \cite{sutton2018reinforcement}. Despite the trial-and-error nature is a peculiarity of RL, it is often a concern and an obstacle for its actual deployment in real-world scenarios \cite{dulac2019challenges}. A solution to this problem is offered by the offline-RL framework \cite{fu2020d4rl} where control policies are entirely learned from static datasets which make unnecessary any further interaction with the environment. However, collecting datasets of ready-to-use demonstrations (states, actions and rewards) represents another source of challenge and, even when both states and actions are available, domain and structural differences between source agent, i.e. the agent from which we collect data, and the target agent (the learner) could hinder the actual usefulness of the dataset. More commonly, offline datasets come in the form of videos and therefore methods for the extraction of usable states, actions and rewards from videos are needed.\\
In this paper, we want to leverage videos of animals navigating the environment as offline datasets. We propose an hybrid approach which leverages agent interactions to extract ready-to-use demonstrations from videos and uses the extracted data to improve RL efficiency and performance. In order to accomplish our goal we address three main problems throughout the paper: $(i)$ dealing with domain adaptation (DA) issues, $(ii)$ dealing with videos as offline datasets and $(iii)$ theoretically motivating off-policy learning algorithms.

\paragraph*{Contributions} We start by formulating an adversarial DA algorithm \cite{tzeng2017adversarial} which maps both the offline dataset and the online interactions of the learning agent into a common subspace, allowing for the use of data from different domains. Further, we show a simple method based on supervised learning to recover demonstrations (states, actions and rewards) from observations (states only) and finally, we re-derive weighted regression type of algorithms such as Relative Entropy Policy Search (REPS) \cite{peters2010relative}, Maximum a posteriori Policy Optimization (MPO) \cite{abdolmaleki2018maximum}, Advantage Weighted Regression (AWR) \cite{peng2019advantage} and Advantage Weighted Actor Critic \cite{nair2020awac} (AWAC) in the light of a recent work on principled off-policy learning \cite{queeney2021generalized}. Note that, with respect to \cite{queeney2021generalized}, we cannot use importance sampling (IS) in our formulation as the behavioral policy generating the data is unknown when learning from videos of animals. Our final algorithm Advantage Weighted Policy Optimization (AWPO) would indeed not require IS. We test our method in navigation tasks for both offline and online RL showing encouraging results.

\paragraph*{Notation}Unless indicated, we use uppercase letters (e.g., $S_t$) for random variables, lowercase letters (e.g., $s_t$) for values of random variables, script letters (e.g., $\mathcal{S}$) for sets, and bold lowercase letters (e.g., $\bm{\theta}$) for vectors. Let $[t_1 : t_2]$ be the set of integers $t$ such that $t_1 \leq t \leq t_2$; we write $S_t$ such that $t_1 \leq t \leq t_2$ as $S_{t_1 : t_2}$. We denote with $\mathbb{E}[\cdot]$ expectation, with $\mathbb{P}(\cdot)$ probability and with $\mathbb{D}_{\text{KL}}(\cdot||\cdot)$ the KL-divergence.

\section{Preliminaries}
We consider an infinite-horizon discounted Markov Decision Process (MDP) defined by the tuple $(\mathcal{S}, \mathcal{A}, P, r, D, \gamma)$ where $\mathcal{S}$ is the set of states and $\mathcal{A}$ is the set of actions, both possibly infinite. $P:\mathcal{S}\times \mathcal{A} \rightarrow \Delta_{\mathcal{S}}$ is the transition probability function and $\Delta_{\mathcal{S}}$ denotes the space of probability distributions over $\mathcal{S}$. The function $r:\mathcal{S}\times \mathcal{A} \rightarrow \mathbb{R}$ maps state-action pairs to rewards. $D\in\Delta_{\mathcal{S}}$ is the initial state distribution and $\gamma \in [0,1)$ the discount factor. We model the decision agent as a stationary policy $\pi:\mathcal{S}\rightarrow\Delta_{\mathcal{A}}$, where $\pi(a|s)$ is the probability of taking action $a$ in state $s$. The goal is to choose a policy that maximizes the expected total discounted reward $J(\pi)=\mathbb{E}_{\tau}[\sum_{t=0}^{\infty}\gamma^t r(s_t,a_t)]$, where $\tau = (s_0,a_0,s_1,a_1,\dots)$ is a trajectory sampled according to $s_0 \sim D$, $a_t\sim\pi(\cdot|s_t)$ and $s_{t+1}\sim P(\cdot|s_t,a_t)$. A policy $\pi$ induces a normalized discounted state visitation distribution $d^{\pi}$, where $d^{\pi}(s) = (1-\gamma)\sum_{t=0}^{\infty}\gamma^t\mathbb{P}(s_t=s | D,\pi,P)$. We define the corresponding normalized discounted state-action visitation distribution as $\mu^{\pi}(s,a)= d^{\pi}(s)\pi(a|s)$. We denote the state value function of $\pi$ as $V^{\pi}(s) = \mathbb{E}_{\tau \sim \pi}[\sum_{t=0}^{\infty}\gamma^t r(s_t,a_t)|S_0=s]$, the state-action value function $Q^{\pi}(s,a) = \mathbb{E}_{\tau \sim \pi}[\sum_{t=0}^{\infty}\gamma^t r(s_t,a_t)|S_0=s, A_0=a]$ and the advantage function as $A^{\pi}(s,a) = Q^{\pi}(s,a) - V^{\pi}(s)$. Finally, when a function is parameterized with parameters $\bm{\theta} \in \varTheta \subset \mathbb{R}^k$ we write $\pi_{\bm{\theta}}$.

\section{Domain Adaptation and actions-rewards estimation from Videos}

Given $s\in \mathcal{S}$ a single or an aggregation of RGB frames, we can extract from videos in a source MDP the offline dataset $(s_S, s_S')_{1:M}$ made only of observations. Further, we assume an isomorphism exists between the source MDP, from which collect the dataset, and the target MDP of the learner.
\begin{assumption}
Given a source MDP defined by the tuple $(\mathcal{S}_S, \mathcal{A}_S, P, r, D, \gamma)$ in which the source agent (an animal in our case) operates and a target MDP defined by $(\mathcal{S}_T, \mathcal{A}_T, P, r, D, \gamma)$, we assume there exist two invertible functions $\iota_o:\mathcal{S}_S \rightarrow \mathcal{S}_T$ and $\iota_a:\mathcal{A}_S \rightarrow \mathcal{A}_T$ such that $P(s'_T|s_T,a_T) = P(\iota_o(s_S')|\iota_o(s_S), \iota_a(a_S))$.
\label{Ass:isomorphism}
\end{assumption}
Under assumption~\ref{Ass:isomorphism}, we can build an algorithm which collects data from the source MDP and learns in the target. Moreover, we assume the dataset $(s_S, s_S')_{1:M}$ depicts an optimal behavior with respect to the reward function $r$ which is shared by both the source and target MDPs.

\subsection{Domain adaptation}
In the following we focus on $\iota_o:\mathcal{S}_S \rightarrow \mathcal{S}_T$ in assumption~\ref{Ass:isomorphism} which provides a connection between the source state space $\mathcal{S}_S$ and the target $\mathcal{S}_T$. In order to leverage data collected in $\mathcal{S}_S$ for learning in $\mathcal{S}_T$, we have two main alternatives: either directly estimating $\iota_o$ or mapping both $\mathcal{S}_S$ and $\mathcal{S}_T$ in a common lower dimensional subspace $\mathcal{H}$. In the specific case of learning from videos, estimating $\iota_o$ can be accomplished using unpaired image-to-image translation as done in \cite{zhu2017unpaired, jha2018disentangling, park2020contrastive}. However, image-to-image translation is an expensive approach which requires large datasets and a conspicuous amount of iterations. As a result, we opt for mapping both $\mathcal{S}_S$ and $\mathcal{S}_T$ in a common encoding subspace $\mathcal{H}$ using the adversarial domain adaptation framework originally proposed in \cite{tzeng2017adversarial}. We define two different encoders $f^{\text{enc}}_{\bm{\delta}}: \mathcal{S}_T \rightarrow \mathcal{H}$ and $f^{\text{enc}}_{\bm{\psi}}: \mathcal{S}_S \rightarrow \mathcal{H}$ and sample from the source MDP $\{s_S^{(i)}\}_{i=1}^m \sim (s_S, s_S')_{1:M}$ and from the target MDP $\{s_T^{(i)}\}_{i=1}^m \sim d^{\pi_{\bm{\theta}_k}}$. Our adversarial DA step consists in the following $\min$-$\max$ game
\begin{align}
    \min_{\bm{\psi}} \max_{\bm{\phi}} \frac{1}{m}\sum_{i=1}^m[ f^{\text{disc}}_{\bm{\phi}}(f^{\text{enc}}_{\bm{\delta}}(s_T^{(i)})) -  f^{\text{disc}}_{\bm{\phi}}(f^{\text{enc}}_{\bm{\psi}}(s_S^{(i)}))],
    \label{eq:adv_loss}
\end{align}
where $f^{\text{disc}}_{\bm{\phi}}:\mathcal{H}\to \mathbb{R}$ is a discriminator function which classifies from which MDP (either target or source) the embedding $h\in\mathcal{H}$ has been generated and the source encoder $f^{\text{enc}}_{\bm{\psi}}:\mathcal{S}_S \rightarrow \mathcal{H}$ works on making the classification step harder for the discriminator. At optimality, the discriminator is no longer able to distinguish from which MDP the embedding $h$ has been generated meaning that both $\mathcal{S}_S$ and $\mathcal{S}_T$ have been mapped into the same common subspace $\mathcal{H}$.  Note that, eq.~\ref{eq:adv_loss} is the Wasserstein loss as introduced in \cite{arjovsky2017wasserstein}. We prefer this loss to the original binary cross-entropy in \cite{goodfellow2014generative} for stability reasons.

\subsection{Actions and Rewards estimation}
Assume we have the following encoders $f^{\text{enc}}_{\bm{\delta}}: \mathcal{S}_T \rightarrow \mathcal{H}$ and $f^{\text{enc}}_{\bm{\psi}}: \mathcal{S}_S \rightarrow \mathcal{H}$ for the target and source MDPs respectively and a set of observations from the source MDP $(s_S, s_S')_{1:M}$. Source actions $(a_S)_{1:M}$ and rewards $(r_S)_{1:M}$ can be estimated via supervised learning. Specifically, we define two inverse models $f^{ia}_{\bm{\chi}}: \mathcal{H} \times \mathcal{H} \rightarrow \Delta_{\mathcal{A}}$ and $f^{ir}_{\bm{\zeta}}: \mathcal{H} \times \mathcal{H} \rightarrow \mathbb{R}$ for actions and rewards respectively, collect learner interactions $(s_T, a_T, s_T', r_T)_{1:N}$ in the target MDP and train the two inverse models in order to minimize the following losses
\begin{align}
\begin{split}
    \mathcal{L}^{ia}(\bm{\chi}) &= ||a_T - f^{ia}_{\bm{\chi}}(f^{\text{enc}}_{\bm{\delta}}(s_T),f^{\text{enc}}_{\bm{\delta}}(s_T'))||^2, \\
    \mathcal{L}^{ir}(\bm{\zeta}) &= ||r_T - f^{ir}_{\bm{\zeta}}(f^{\text{enc}}_{\bm{\delta}}(s_T),f^{\text{enc}}_{\bm{\delta}}(s_T'))||^2.
\end{split}
\label{eq:actions-rewards-loss}
\end{align}
We then sample $\hat{a}_S \sim f^{ia}_{\bm{\chi}}(f^{\text{enc}}_{\bm{\psi}}(s_S),f^{\text{enc}}_{\bm{\psi}}(s_S'))$ and compute $\hat{r}_S = f^{ir}_{\bm{\zeta}}(f^{\text{enc}}_{\bm{\psi}}(s_S),f^{\text{enc}}_{\bm{\psi}}(s_S'))$ obtaining a set of demonstrations $(f^{\text{enc}}_{\bm{\psi}}(s_S), \hat{a}_S, f^{\text{enc}}_{\bm{\psi}}(s_S'), \hat{r}_S)_{1:M}$ suitable for learning in the target MDP. Note that, the learning agent will end-up using the following set of joint offline data and online interactions
\begin{align}
    (f^{\text{enc}}_{\bm{\psi}}(s_S), \hat{a}_S, f^{\text{enc}}_{\bm{\psi}}(s_S'), \hat{r}_S) \ \cup \ (f^{\text{enc}}_{\bm{\delta}}(s_T), a_T, f^{\text{enc}}_{\bm{\delta}}(s_T'), r_T)
    \label{eq:learning_data}
\end{align}
where both $s_S$ and $s_T$, the RGB frames observations, are encoded in the embedding space $\mathcal{H}$ for domain adaptation purposes. Moreover, recall we assume the demonstrating agent in the source MDP to be a quasi-optimal expert, we add therefore a small constant bonus $r_i$ to $\hat{r}_S$ in order to encourage the learning agent to follow the demonstrating agent and to compensate for reward sparsity. This solution, based on inverse models, learner's interactions and supervised learning, has been quite studied in the literature and we refer to \cite{torabi2018behavioral, schmeckpeper2020reinforcement, pathak2017curiosity} for more details.

\section{Advantage Weighted Policy Optimization}
After addressing domain adaptation and actions-rewards estimation, we obtain the set of data in \eqref{eq:learning_data} made of joint offline data and online interactions. The final challenge in order to successfully leverage videos of animals in RL is determining the proper algorithm able to learn from off-policy data. We start by stating some important results instrumental for our derivations.

\paragraph*{Policy improvement lower bound} Our starting point is the policy improvement lower bound originally developed in \cite{kakade2002approximately} and refined in \cite{achiam2017constrained}.
\begin{Theorem}[\cite{achiam2017constrained}]
Consider the current policy $\pi_k$. For any future policy $\pi$ we have
\begin{align}
    \begin{split}
            J(\pi) - J(\pi_k) \geq& \frac{1}{1-\gamma} \mathbb{E}_{(s,a) \sim \mu^{\pi_k}}\bigg[\frac{\pi(a|s)}{\pi_k(a|s)}A^{\pi_k}(s,a)\bigg] \\
            &- \frac{2\gamma C^{\pi,\pi_k}}{(1-\gamma)^2}\mathbb{E}_{s\sim d^{\pi_k}}[\text{TV}(\pi_k, \pi)(s)],
    \label{eq:achiam}
    \end{split}
\end{align}
where $C^{\pi,\pi_k} = \max_{s\in \mathcal{S}} |\mathbb{E}_{a\sim \pi(\cdot|s)}[A^{\pi_k}(s,a)]|$ and $\text{TV}(\pi_k, \pi) = \frac{1}{2}\int_{a \in \mathcal{A}}|\pi_k(a|s) - \pi(a|s)| da$ is the total variation distance between the distributions $\pi_k(\cdot|s)$ and $\pi(\cdot|s)$.
\label{Theo:Achiam}
\end{Theorem}

Theorem~\ref{Theo:Achiam} has been further generalized for the off-policy setting in \cite{queeney2021generalized} in order to allow for expectations with respect to any behavioral policy.
\begin{Theorem}[\cite{queeney2021generalized}]
Consider the current policy $\pi_k$ and the behavioral policy $\beta$. For any future policy $\pi$ we have
\begin{align}
    \begin{split}
        J(\pi) - J(\pi_k) \geq& \frac{1}{1-\gamma}\mathbb{E}_{(s,a) \sim \mu^{\beta}}\bigg[\frac{\pi(a|s)}{\beta(a|s)}A^{\pi_k}(s,a)\bigg] \\
        &- \frac{2\gamma C^{\pi, \pi_k}}{(1-\gamma)^2}\mathbb{E}_{s \sim d^{\beta}}[\text{TV}(\beta, \pi)(s)]
        \label{eq:lemma_GePPO}
    \end{split}
\end{align}
where $C^{\pi,\pi_k}$ and $\text{TV}(\beta, \pi)$ are defined as in Theorem~\ref{Theo:Achiam}.
\label{Theo:Queeney}
\end{Theorem}

The first term on the right-hand side in both \eqref{eq:achiam} and \eqref{eq:lemma_GePPO} is referred to as the surrogate objective, while the second is a penalty term or a soft constraint. Note that both Theorem~\ref{Theo:Achiam} and Theorem~\ref{Theo:Queeney} show that, enforcing any future policy $\pi$ to stay close to the behavioral policy $\pi_k$ or $\beta$, a monotonic policy improvement is guaranteed. For practical purposes, the following surrogate optimization problem for Theorem~\ref{Theo:Achiam} is derived in \cite{schulman2015trust}
\begin{align}
\begin{split}
    \max_{\pi} \ \ \ & \mathbb{E}_{(s,a) \sim \mu^{\pi_k}}\bigg[\frac{\pi(a|s)}{\pi_k(a|s)}A^{\pi_k}(s,a)\bigg] , \\
    &\mathbb{E}_{s \sim d^{\pi_k}}[\mathbb{D}_{\text{KL}}(\pi_k(\cdot|s) || \pi(\cdot|s))] \leq \epsilon.
    \label{eq:TRPO_surrogate_problem}
\end{split}
\end{align}
Based on \eqref{eq:TRPO_surrogate_problem}, the popular on-policy Proximal Policy Optimization (PPO) \cite{schulman2017proximal} and its off-policy version Generalized PPO (GePPO) in \cite{queeney2021generalized} are formalized. We provide a thorough review on the state-of-art algorithms commonly adopted in deep RL and of their connections with Theorems~\ref{Theo:Achiam} and \ref{Theo:Queeney} in the appendix.

\paragraph*{AWPO derivations} In the following, we focus on deriving from Theorem~\ref{Theo:Queeney} an off-policy algorithm which we call Advantage Weighted Policy Optimization (AWPO). This belongs to the family of weighted regression algorithms and have strong similarities with REPS, MPO, AWR and AWAC to mention some of them. The main source of difference lies in the way the algorithm is derived and in the need for an off-policy Advantage Estimation.
We start from the lower bound in Theorem~\ref{Theo:Queeney} and proceed analogously to \cite{achiam2017constrained}. Note that, the total variation distance is a metric \cite{tsybakov1997nonparametric}, hence $\text{TV}(\beta, \pi)(s) = \text{TV}(\pi, \beta)(s)$. Through Pinsker's inequality \cite{chandak2021universal} we get $\text{TV}(\pi, \beta)(s) \leq \sqrt{\mathbb{D}_{\text{KL}}(\pi|| \beta)(s)/2}$ and combined with Jensen's inequality yields
\begin{align*}
    \begin{split}
        \mathbb{E}_{s \sim d^{\beta}}[\text{TV}(\pi, \beta)(s)] &\leq \mathbb{E}_{s \sim d^{\beta}}\bigg[\sqrt{\frac{1}{2}\mathbb{D}_{\text{KL}}(\pi|| \beta)(s)}\bigg] \\
        &\leq \sqrt{\frac{1}{2}\mathbb{E}_{s \sim d^{\beta}}[\mathbb{D}_{\text{KL}}(\pi|| \beta)(s)]},    \end{split}
\end{align*}
from which we get the surrogate problem in \eqref{eq:AWAC_surrogate} where $\beta$ is any behavioral policy.
\begin{align}
    \begin{split}
            \max_{\pi} \ \ \ & \mathbb{E}_{s \sim d^{\beta}}\bigg[\mathbb{E}_{a \sim \pi(\cdot|s)}\big[A^{\pi_k}(s,a)\big]\bigg], \\
            &\mathbb{E}_{s \sim d^{\beta}}[\mathbb{D}_{\text{KL}}(\pi(\cdot|s) || \beta(\cdot|s))] \leq \epsilon.
    \end{split}
    \label{eq:AWAC_surrogate}
\end{align}
By converting the KL-divergence in \eqref{eq:AWAC_surrogate} as a soft constraint and solving the Lagrangian problem in closed form, we obtain the following optimal solution
\begin{align*}
    \pi^*(a|s) &= \beta(a|s) \exp \bigg( \frac{A^{\pi_k}(s,a)}{\lambda} \bigg)\frac{1}{Z(s)},
\end{align*}
where $Z(s)$ is a normalizing factor ensuring $\pi^*(\cdot|s)$ is a probability distribution.
Further, we project $\pi^*$ on the manifold of function approximations parameterized by $\bm{\theta}$. Given $\pi_{\bm{\theta}_k}$ the policy parameterized by $\bm{\theta}_k$ at step $k$, for $k+1$ we obtain
\begin{align}
    \bm{\theta}_{k+1} = \arg \max_{\bm{\theta}} \mathbb{E}_{(s,a) \sim \mu^{\beta}}\bigg[\log \pi_{\bm{\theta}}(a|s)\exp \bigg( \frac{A^{\pi_{\bm{\theta}_{k}}}(s,a)}{\lambda} \bigg)\bigg].
\end{align}
We defer to the appendix more details on this derivation.
\begin{remark}[Novelty]
    We show direct connection with the generalized policy improvement in Theorem~\ref{Theo:Queeney} which is a missing piece in all the aforementioned literature \cite{peters2010relative,abdolmaleki2018maximum, peng2019advantage, nair2020awac}. Moreover, a novel small nuance of our formulation is the need for off-policy policy evaluation since we make use of samples from $\beta$ to compute $A^{\pi_{\bm{\theta}_{k}}}$ for the current policy $\pi_{\bm{\theta}_{k}}$.
\end{remark}
\begin{remark}[Off-policy policy evaluation]
    Off-policy policy evaluation is defined as leveraging data generated by a behavioral policy $\beta$ to evaluate the performance metric of the evaluation policy $\pi_{\bm{\theta}_{k}}$ (cf. \cite{chandak2021universal}). Well known algorithms for off-policy policy evaluation are Retrace in \cite{munos2016safe}, V-trace in \cite{espeholt2018impala} and UnO in \cite{chandak2021universal}. All of these require however importance sampling (IS) correction, which implies knowledge of $\beta$. When learning from videos of animals, $\beta$ is unknown and cannot be estimated using Imitation Learning as in \cite{kostrikov2021offline, giammarino2022learning} given the absence of actions. A preferable alternative is to avoid IS correction. Thus, we identify and test three main algorithms for off-policy evaluation which do not perform IS correction, these are: Tree-Backup Lambda (TBL) \cite{precup2000eligibility}, Harutyunyan Q-Lambda (HQL) \cite{harutyunyan2016q} and Peng Q-Lambda (PQL) \cite{peng1994incremental}. We refer to \cite{chandak2021universal,,munos2016safe,kozuno2021revisiting} for additional details on policy evaluation with off-policy correction.
\end{remark}
The full algorithm, consisting in: DA step, source actions-rewards estimation and off-policy RL is summarized in Algorithm~\ref{alg:AWPO}.

\begin{algorithm}
\caption{AWPO from Videos w/ Adv Domain Adaptation}
\label{alg:AWPO}
Input: $(s_S, s_S')_{1:M}$, $\pi_{\bm{\theta}}$, $f^{\text{enc}}_{\bm{\psi}}$, $f^{\text{enc}}_{\bm{\delta}}$, $f^{\text{disc}}_{\bm{\phi}}$, $f^{ia}_{\bm{\chi}}$, $f^{ir}_{\bm{\zeta}}$, $r_i$, $\alpha$, $\lambda$, $K$, $N$.  \\
\For{$k=1,K-1$}{
$(s_T, a_T, s_T', r_T)_{1:N} \sim \pi_{\bm{\theta}_k}$ \\
Adv DA in \eqref{eq:adv_loss},\\
$\hat{a}_S$ and $\hat{r}_S$ estimation as in \eqref{eq:actions-rewards-loss},\\
set $\hat{r}_S =  f^{ir}_{\bm{\zeta}}(f^{\text{enc}}_{\bm{\psi}}(s_S),f^{\text{enc}}_{\bm{\psi}}(s_S')) + r_i$ \\
Given the data $\mathcal{D}$ in \eqref{eq:learning_data} where $h \leftarrow f^{\text{enc}}(s)$:\\
Perform $A^{\pi_{\bm{\theta}_{k}}}(h,a)$ estimation,\\
$\bm{\theta}_{k+1} = \bm{\theta}_{k} + \alpha \nabla_{\bm{\theta}} \mathbb{E}_{\mathcal{D}}\bigg[\log \pi_{\bm{\theta}}(a|h)\exp \bigg( \frac{A^{\pi_{\bm{\theta}_{k}}}(h,a)}{\lambda} \bigg)\bigg]$\\
$\bm{\theta}_{k} \leftarrow \bm{\theta}_{k+1}$
}
$\bm{\theta} \leftarrow \bm{\theta}_{K}$\\
\Return $\pi_{\bm{\theta}}$
\end{algorithm}

\section{Experiments}
We conduct a series of experiments in order to gradually evaluate all the single pieces of our algorithm. We focus on 2D navigation tasks with sparse rewards and run tests on offline RL in order to assess the offline performance of AWPO with different Advantage Estimation algorithms. Then, we continue testing Algorithm~\ref{alg:AWPO} in RL with offline observations without DA, with simplified DA (environment in Fig.~\ref{fig:human} as target and Fig.~\ref{fig:human_modified} as source for both these experiments) and finally with hard DA where the rodent video in Fig.~\ref{fig:rodent_video} is used as offline dataset (Fig.~\ref{fig:human} for target and Fig.~\ref{fig:rodent} for source). For all the implementations details for both the \href{https://github.com/VittorioGiammarino/Offline_RL_pytorch.git}{offline}\footnote{https://github.com/VittorioGiammarino/Offline\_RL\_pytorch.git} and \href{https://github.com/VittorioGiammarino/Leveraging-Videos-of-Rodent-Behavior-to-Accelerate-Reinforcement-Learning-from-Pixels-in-Navigation.git}{online}\footnote{https://github.com/VittorioGiammarino/Leveraging-Videos-of-Rodent-Behavior-to-Accelerate-Reinforcement-Learning-from-Pixels-in-Navigation.git} experiments we defer to the code. In a nutshell, the encoders are $3$ layers Convolutional Neural Networks (CNN) and discriminator, policy network, critic network and inverse models are MultiLayer Perceptrons (MLP) with a single hidden layer.
\begin{figure*}[h!]
    \centering
    \begin{subfigure}[t]{0.3\textwidth}
        \centering
        \includegraphics[width=5cm, height=5cm]{Figures/human_expert.pdf}
        \caption{}
        \label{fig:human}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.3\textwidth}
        \centering
        \includegraphics[width=5cm, height=5cm]{Figures/modified_human_expert.pdf}
        \caption{}
        \label{fig:human_modified}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.3\textwidth}
        \centering
        \includegraphics[width=5cm, height=5cm]{Figures/rodent.pdf}
        \caption{}
        \label{fig:rodent}
    \end{subfigure}
    \caption{Environments used for our online RL experiments. Fig.~\ref{fig:human} shows the target MDP. Fig~\ref{fig:human_modified} and Fig.~\ref{fig:rodent} show two different source MDP from which we collect videos used as offline datasets in Algorithm~\ref{alg:AWPO}. For both the environments in \ref{fig:human} and Fig~\ref{fig:human_modified} we rely on the open-source library minigrid \cite{minigrid}}
    \label{fig:isomorphism}
\end{figure*}

\begin{figure*}[h!]
    \centering
    \includegraphics[width=\linewidth]{Figures/frames_2030.pdf}
    \caption{Frames from the video used as offline dataset. The video consists in a top view of a rodent foraging in an empty maze.}
    \label{fig:rodent_video}
\end{figure*}

\subsection{Offline RL}

In the following we test AWPO with different advantage estimation techniques. In addition to the aforementioned HQL \cite{harutyunyan2016q}, PQL \cite{peng1994incremental} and TBL \cite{precup2000eligibility}, which perform IS free off-policy correction to compensate for the use of off-policy data without knowing $\beta$, we also test Generalized Advantage Estimation (GAE) \cite{schulman2015high} which does not perform off-policy correction. Moreover, we compare AWPO against the offline RL baselines Soft-Actor Critic (SAC) \cite{haarnoja2018soft}, Behavioral Cloning (BC) \cite{osa2018algorithmic}, AWAC \cite{nair2020awac} and the state-of-the-art TD3+BC \cite{fujimoto2021minimalist}. We focus on the navigation tasks available in the popular library "Datasets for Deep Data-Driven RL" \cite{fu2020d4rl}. The full results are in Table~\ref{tab:table_offline} and show all the different versions of AWPO matching SOTA performance. Note that, in these experiments, the datasets are made of demonstrations (state, actions and rewards) and do not require any DA step.

\begin{table}
\centering
\caption{Offline RL experiments. The reported reward is averaged over the final 10 evaluations and 10 seeds. Each trial consists of $200$ iterations and during each iteration $30$k offline samples are used for training.}
\label{tab:table_offline}
\footnotesize
\begin{tabular}{c c c c}\toprule
& & \multicolumn{2}{c}{medium} \\
\cmidrule(lr){3-4}
& maze2d-large & antmaze-diverse & antmaze-play \\
\bottomrule
BC & $277.3 \pm 134.8$ & $0.35 \pm 0.33$ & $0.46 \pm 0.28$ \\
AWPO+GAE & $214.7 \pm 126.5$ & $0.43 \pm 0.29$ & $0.48 \pm 0.18$ \\
AWAC & $270.6 \pm 144.6$ & - & - \\
AWPO+HQL & $262.8 \pm 109.8$ & $0.44$ $\pm 0.16$ & $0.43 \pm 0.2$\\
AWPO+PQL & $300.3 \pm 74.2$ & $0.43 \pm 0.27$ & $0.63$ $\pm 0.26$ \\
AWPO+TBL & $204.2 \pm 98.3$ & $0.36 \pm 0.21$ & $0.46 \pm 0.24$ \\
SAC & $186.02 \pm 105.5$ & - & - \\
TD3+BC & $321.1$ $\pm 157.7$ & $0.16 \pm 0.2$ & $0.43 \pm 0.35$ \\
\bottomrule
\end{tabular}
\end{table}

\subsection{Online RL}

\begin{table}
\centering
\caption{Online RL experiments. The reported reward is averaged over the final 10 evaluations and 10 seeds. Each experiment consists of $400$k training steps. In RL only we focus on simple RL where we make no use of off-policy data or offline datasets. In RL+\ref{fig:human_modified} we collect an expert video in Fig.~\ref{fig:human_modified} and leverage this dataset without explicit DA step. Similarly in RL+\ref{fig:human_modified}+DA where however we use a DA step and set $r_i=0.01$. Finally, in RL+\ref{fig:rodent}+DA (1) and RL+\ref{fig:rodent}+DA (2) we use the video in Fig.~\ref{fig:rodent_video} as offline dataset and we set $r_i=0.01$ for the former and $r_i=0.005$ for the latter.}
\label{tab:table_online}
\footnotesize
\begin{tabular}{c c c c c}\toprule
& \multicolumn{4}{c}{AWPO +} \\
\cmidrule(lr){2-5}
& GAE & HQL & PQL & TBL \\
\bottomrule
RL only & $.71 \pm .09$ & $.02 \pm .04$ & $.01 \pm .03$ & $0$ \\
RL+\ref{fig:human_modified} & $.04 \pm .07$ & $.08 \pm .16$ & $.13 \pm .20$ & $.03 \pm .06$ \\
RL+\ref{fig:human_modified}+DA & $.44 \pm .29$ & $.36 \pm .27$ & $.35 \pm .26$ & $.26 \pm .23$ \\
RL+\ref{fig:rodent}+DA (1) & $.27 \pm .26$ & $.28 \pm 0.24$ & $.39 \pm .24$ & $.25 \pm .16$ \\
RL+\ref{fig:rodent}+DA (2) & $.42 \pm .32$ & $.26 \pm .23$ & $.21 \pm .18$ & $.27 \pm .26$ \\
\bottomrule
\end{tabular}
\end{table}

We focus on the 2D navigation task in Fig.~\ref{fig:human} as target MDP, where the state $s_T \in \mathcal{S}_T$ is the $128 \times 128$ RGB image showing the full grid and the goal is to reach the bottom right corner labelled by the green square. We run all the experiments for $400k$ training steps and evaluate the policy every $4k$ steps for $10$ episodes. This procedure is repeated for $10$ different random seeds. Table~\ref{tab:table_online} summarizes the final results where we report the reward obtained at the end of the training process and averaged over episodes and seeds. Note that, almost all the different versions of AWPO benefit from the introduction of videos as offline dataset when associated with DA. In particular, Algorithm~\ref{alg:AWPO} combined with AWPO+HQL, AWPO+PQL and AWPO+TBL always outperforms its RL only counterpart. However, when we turn the DA step off, as in the $2$nd row in Table~\ref{tab:table_online}, we observe limited benefit in leveraging offline videos. The only exception is given by AWPO+GAE which shows great stability in its RL only version (only $.09$ of standard deviation) and higher instability when associated with offline datasets which leads to lower average performance. We include plots of the training curves for each experiment in the appendix.

\section{Conclusions}
This paper analyzes the challenges of using videos in RL as offline datasets in order to increase learning efficiency and performance. The main problems we identify consist in dealing with domain discrepancies between source and target domains, learning from offline observations rather than demonstrations and being able to opportunely accomplish off-policy learning. For domain discrepancies we mean using data collected in conditions which are not necessarily the one encountered by the learning agent. As described in Assumption~\ref{Ass:isomorphism}, we allow for structural differences in $\mathcal{S}$ and $\mathcal{A}$ and our DA step works on smoothing these discrepancies transforming the offline dataset in order to make it exploitable by the learning agent. The importance of this step is evident in Table~\ref{tab:table_online} for the Fig.~\ref{fig:human_modified} experiments. Furthermore, learning from offline observations deals with the problem of estimating actions and rewards from static dataset made of observations only. Our solutions is based on inverse models trained using supervised learning and learner's interactions. Finally, considering recent works on principled off-policy learning, we re-derive the AWPO algorithm and test it with and without offline videos on a series of challenging scenarios.
Future works will focus on improving each individual step and test our solutions not only in navigation but also for learning locomotion and manipulation.

\printbibliography

\section{Appendix}

\subsection{An overview of online RL algorithms}
In this section, we provide an overview of the RL algorithms we considered throughout the paper. We focus on actor-critic structures \cite{konda1999actor} where in addition to a policy network $\pi_{\bm{\theta}}$, parameterized by a function approximation with weights $\bm{\theta}$, we also use a critic network, either $V_{\bm{\theta}_V}(\bm{s})$ or $Q_{\bm{\theta}_V}(\bm{s}, \bm{a})$, parameterized by $\bm{\theta}_V$ or $\bm{\theta}_Q$. All the algorithms presented here update the actor via policy gradient \cite{sutton1999policy} and the critic by minimizing the mean-squared error between the critic network value and a target value as in \cite{mnih2013playing}.

We classify actor-critic algorithms in three main categories: pure off-policy algorithms, which rely on old transitions stored in a replay buffer; pure on-policy algorithms, which use only transitions generated by the current policy; and hybrid methods, which can be easily adapted to both the updates.

\subsection{Off-policy methods}
We consider SAC and Twin Deep Delayed Deterministic Policy Gradient (TD3) \cite{fujimoto2018addressing}. Both algorithms, at each iteration, sample transitions $(s, a, r, s', a')$ from a replay buffer $\mathcal{D}$ and solve a version of the following problem
\begin{align}
    \bm{\theta}_{k+1} = \arg\max_{\bm{\theta}_k} \ \ \mathbb{E}_{\bm{s} \sim \mathcal{D}}\big[\mathbb{E}_{a \sim \pi_{\bm{\theta}_k}}\big[Q_{\bm{\theta}_Q}(s, a)\big]\big],
\end{align}
where $\bm{\theta}_Q$ is computed as follows
\begin{align}
    \begin{split}
        y_k &= \mathbb{E}_{s' \sim \mathcal{D}}\big[r + \gamma\mathbb{E}_{a' \sim \pi_{\bm{\theta}_k}}\big[Q_{\bm{\theta}_Q}(s', a')\big]\big] \\
        \bm{\theta}_Q &= \arg\min_{\bm{\theta}_Q} \ \ \mathbb{E}_{(s, a)\sim\mathcal{D}}\big[(y_k - Q_{\bm{\theta}_Q}(s, a))^2\big].
    \end{split}
    \label{eq:off_policy_critic_update}
\end{align}
Specifically, SAC augments the objective with an entropy term $\mathcal{H} = \mathbb{E}_{a \sim \pi_{\bm{\theta}_k}}\big[-\log\pi_{\bm{\theta}_k}(s, a)\big]$ to incentive wider exploration while TD3 uses a series of mechanisms to address function approximation error. We defer to \cite{haarnoja2018soft} and \cite{fujimoto2018addressing} for more details. These approaches are really successful in practice but are not backed up by some theoretical analysis.

\subsection{On-policy methods}
We consider Proximal Policy Optimization (PPO) in \cite{schulman2017proximal}.

Given a current policy $\pi_k$, as showed in \cite{kakade2002approximately}, for any future policy $\pi$ we have
\begin{equation}
    J(\pi) - J(\pi_k) = \frac{1}{1-\gamma}\mathbb{E}_{s \sim d^{\pi}}\bigg[\mathbb{E}_{a\sim\pi(\cdot|s)}[A^{\pi_k}(s,a)]\bigg],
\end{equation}
From which follows the result in Theorem~\ref{Theo:Achiam}
\begin{align}
\begin{split}
    J(\pi) - J(\pi_k) &\geq \frac{1}{1-\gamma} \mathbb{E}_{(s,a) \sim \mu^{\pi_k}}\bigg[\frac{\pi(a|s)}{\pi_k(a|s)}A^{\pi_k}(s,a)\bigg] \\
    &- \frac{2\gamma C^{\pi,\pi_k}}{(1-\gamma)^2}\mathbb{E}_{s\sim d^{\pi_k}}[\text{TV}(\pi_k, \pi)(s)],
    \label{eq:app_achiam}
\end{split}
\end{align}
where TV is the total variation distance defined as
\begin{equation}
    \text{TV}(\pi_k,\pi)(s) = \frac{1}{2}\int_{a \in \mathcal{A}}|\pi_k(a|s) - \pi(a|s)| da.
    \label{eq:Achiam}
\end{equation}
The inequality in \eqref{eq:app_achiam} shows that, provided a small TV between the current and future policies, the difference in expected return (aka expected improvement) is lower bounded by a positive value, which means we can guarantee a monotonic improvement throughout policy gradient iterations.
By further combining Pinsker's inequality with Jensen's inequality the following holds
\begin{align*}
    \mathbb{E}_{s\sim d^{\pi_k}}[\text{TV}(\pi_k, \pi)(s)] &\leq \mathbb{E}_{s\sim d^{\pi_k}}\bigg[\sqrt{\frac{1}{2}\mathbb{D}_{\text{KL}}(\pi_k || \pi)} \bigg] \\
    &\leq \sqrt{\frac{1}{2}\mathbb{E}_{s\sim d^{\pi_k}}\bigg[\mathbb{D}_{\text{KL}}(\pi_k || \pi)} \bigg].
\end{align*}
The following problem is therefore proposed in \cite{schulman2015trust} to maximize the lower bound in \eqref{eq:app_achiam}:
\begin{align}
\begin{split}
    \max_{\pi} \ \ \ & \mathbb{E}_{(s,a) \sim \mu^{\pi_k}}\bigg[\frac{\pi(a|s)}{\pi_k(a|s)}A^{\pi_k}(s,a)\bigg] , \\
    &\mathbb{E}_{s \sim d^{\pi_k}}[\mathbb{D}_{\text{KL}}(\pi_k(\cdot|s) || \pi(\cdot|s))] \leq \epsilon.
\end{split}
    \label{eq:app_TRPO_surrogate_problem}
\end{align}
Furthermore, converting the hard KL constraint in a soft constraint we get
\begin{align*}
    \pi =& \arg\max_{\pi} \mathbb{E}_{s \sim d^{\pi_k}}\bigg[\mathbb{E}_{a \sim \pi(\cdot|s)}\big[A^{\pi_k}(s,a)\big]\bigg] \\
    &+ \alpha \bigg( \epsilon - \mathbb{E}_{s \sim d^{\pi_k}}[\mathbb{D}_{\text{KL}}(\pi_k(\cdot|s) || \pi(\cdot|s))] \bigg) \\
    =& \arg\max_{\pi}\mathbb{E}_{s \sim d^{\pi_k}}\bigg[\mathbb{E}_{a \sim \pi(\cdot|s)}\big[A^{\pi_k}(s,a)\big]\bigg] \\
    &+ \alpha\epsilon - \alpha\mathbb{E}_{s \sim d^{\pi_k}}\bigg[\int_{a\in\mathcal{A}} \pi_k(a|s)\log\frac{\pi_k(a|s)}{\pi(a|s)} da\bigg] \\
    =& \arg\max_{\pi}\mathbb{E}_{s \sim d^{\pi_k}}\bigg[\mathbb{E}_{a \sim \pi(\cdot|s)}\big[A^{\pi_k}(s,a)\big]\bigg] \\
    &- \alpha \mathbb{E}_{s \sim d^{\pi_k}}\bigg[\int_{a\in\mathcal{A}} \pi_k(a|s)\log\frac{1}{\pi(a|s)} da \bigg] \\
    =& \arg\max_{\pi}\mathbb{E}_{(s,a) \sim \mu^{\pi_k}}\bigg[\frac{\pi(a|s)}{\pi_k(a|s)}A^{\pi_k}(s,a) \\
    &+ \alpha\log\pi(a|s)\bigg] \\
    =& \arg\max_{\pi}\mathbb{E}_{s \sim d^{\pi_k}}\bigg[\mathbb{E}_{a \sim \pi_k}\bigg[\frac{\pi(a|s)}{\pi_k(a|s)}A^{\pi_k}(s,a)\bigg] \\
    &- \alpha\mathcal{H}(\pi_k, \pi)\bigg].
\end{align*}
PPO \cite{schulman2017proximal} converts the previous into the following
\begin{align*}
    \max_{\pi} \mathcal{L}(\pi, \epsilon) =& \max_{\pi}\mathbb{E}_{(s,a) \sim \mu^{\pi_k}}\Bigg[\min\Bigg(\frac{\pi(a|s)}{\pi_k(a|s)}A^{\pi_k}(s,a), \\ &\text{clip}\bigg(\frac{\pi(a|s)}{\pi_k(a|s)}, 1-\epsilon, 1+\epsilon\bigg)A^{\pi_k}(s,a)\Bigg)\Bigg].
\end{align*}
Note that, given the use of importance sampling in the derivation, PPO is an algorithm which necessarily requires knowledge of the policy generating the data $\pi_k$.

\subsection{AWPO derivations}
In the following we derive AWPO. Similar steps lead to AWR \cite{peng2019advantage}, AWAC \cite{nair2020awac} and Generalized Proximal Policy Optimization (GePPO) \cite{queeney2021generalized}.
We start again from the bound in Theorem~\ref{Theo:Achiam}. Note that, The total variation distance is a metric \cite{tsybakov1997nonparametric}, which means
\begin{itemize}
    \item $\text{TV}(\pi,\pi_k)(s) = 0 \text{ iff } \pi = \pi_k \ \ \ \forall a \in \mathcal{A}$,
    \item $\text{TV}(\pi,\pi_k)(s) = \text{TV}(\pi_k,\pi)(s) $,
    \item $\text{TV}(\pi,\pi_k)(s) \leq \text{TV}(\pi,\pi')(s) + \text{TV}(\pi',\pi_k)(s) $.
\end{itemize}
Considering the symmetric property (second bullet in the previous) eq.~\eqref{eq:app_achiam} can be rewritten as
\begin{align*}
    J(\pi) - J(\pi_k) &\geq \frac{1}{1-\gamma} \mathbb{E}_{(s,a) \sim \mu^{\pi_k}}\bigg[\frac{\pi(a|s)}{\pi_k(a|s)}A^{\pi_k}(s,a)\bigg] \\
    &- \frac{2\gamma C^{\pi,\pi_k}}{(1-\gamma)^2}\mathbb{E}_{s\sim d^{\pi_k}}[\text{TV}(\pi,\pi_k)(s)]
\end{align*}
and the following holds as well
\begin{align*}
    \mathbb{E}_{s\sim d^{\pi_k}}[\text{TV}(\pi, \pi_k)(s)] &\leq \mathbb{E}_{s\sim d^{\pi_k}}\bigg[\sqrt{\frac{1}{2}\mathbb{D}_{\text{KL}}(\pi || \pi_k)} \bigg] \\
    &\leq \sqrt{\frac{1}{2}\mathbb{E}_{s\sim d^{\pi_k}}\bigg[\mathbb{D}_{\text{KL}}(\pi || \pi_k)} \bigg].
\end{align*}
These passages yield the following surrogate problem
\begin{align}
\begin{split}
    \max_{\pi} \ \ \ & \mathbb{E}_{s \sim d^{\pi_k}}\bigg[\mathbb{E}_{a \sim \pi(\cdot|s)}\big[A^{\pi_k}(s,a)\big]\bigg], \\
    &\mathbb{E}_{s \sim d^{\pi_k}}[\mathbb{D}_{\text{KL}}\big(\pi(\cdot|s) || \pi_k(\cdot|s)\big)] \leq \epsilon.
\end{split}
\label{eq:app_surrogate2}
\end{align}
Note that $\mathbb{D}_{\text{KL}}\big(\pi(\cdot|s) || \pi_k(\cdot|s)\big) \neq \mathbb{D}_{\text{KL}}\big(\pi_k(\cdot|s) || \pi(\cdot|s)\big)$. Therefore, despite having the same starting point, \eqref{eq:app_TRPO_surrogate_problem} and \eqref{eq:app_surrogate2} lead to different algorithms.

We proceed as illustrated in \cite{peters2010relative}, the hard KL constraint can be converted in a soft constraint, yielding
\begin{align*}
    \mathcal{L}(\pi, \lambda) &= \mathbb{E}_{s \sim d^{\pi_k}}\bigg[\mathbb{E}_{a \sim \pi(\cdot|s)}\big[A^{\pi_k}(s,a)\big]\bigg] \\
    &+ \lambda \bigg( \epsilon - \mathbb{E}_{s \sim d^{\pi_k}}[\mathbb{D}_{\text{KL}}(\pi(\cdot|s) || \pi_k(\cdot|s))] \bigg).
\end{align*}
We can rewrite this problem as
\begin{align*}
    \pi =& \arg\max_{\pi}\mathbb{E}_{s \sim d^{\pi_k}}\bigg[\mathbb{E}_{a \sim \pi(\cdot|s)}\big[A^{\pi_k}(s,a)\big]\bigg] \\
    &- \lambda \bigg(\mathbb{E}_{s \sim d^{\pi_k}}\bigg[\int_{a \in \mathcal{A}}\pi(a|s)\log \frac{\pi(a|s)}{\pi_k(a|s)}\bigg]\bigg), \\
    =& \arg\max_{\pi}\mathbb{E}_{s \sim d^{\pi_k}}\bigg[\mathbb{E}_{a \sim \pi(\cdot|s)}\big[A^{\pi_k}(s,a) \\
    &- \lambda\log\pi(a|s) +\lambda \log\pi_k(a|s)) \big]\bigg], \\
    =& \arg\max_{\pi}\mathbb{E}_{s \sim d^{\pi_k}}\bigg[\mathbb{E}_{a \sim \pi(\cdot|s)}\big[A^{\pi_k}(s,a)\big] \\
    &+ \lambda\mathcal{H}(\pi) -\lambda\mathcal{H}(\pi, \pi_k)\big]\bigg].
\end{align*}
Given all the terms expressed as an expectation of the functional $\pi(\cdot|s)$, we can solve the problem in closed form differentiating wrt $\pi$
\begin{align*}
    \frac{\partial  \mathcal{L}(\pi, \lambda)}{\partial \pi} =& d^{\pi_k}(s)A^{\pi_k}(s,a) - \lambda d^{\pi_k}(s) \log \pi(a|s) \\
    &+ \lambda d^{\pi_k}(s) \log \pi_k(a|s) - \lambda d^{\pi_k}(s).
\end{align*}
We look for a stationary point by setting the differentiation equal to zero and we obtain
\begin{align*}
    \log \pi(a|s) &= \frac{A^{\pi_k}(s,a)}{\lambda} +\log \pi_k(a|s) - Z(s),  \\
    \pi^*(a|s) &= \pi_k(a|s) \exp \bigg( \frac{A^{\pi_k}(s,a)}{\lambda} \bigg)\frac{1}{Z(s)}.
\end{align*}
Using a function approximation we need to project $\pi^*$ on the function manifold as follows
\begin{align*}
    \theta^* =& \arg \min_{\theta} \mathbb{E}_{s \sim d^{\pi_k}}[\mathbb{D}_{\text{KL}}(\pi^* || \pi_{\theta})] \\
    =& \arg \min_{\theta} \mathbb{E}_{s \sim d^{\pi_k}}\bigg[\mathbb{D}_{\text{KL}}\bigg(\pi_k(a|s) \\
    &\exp \bigg( \frac{A^{\pi_k}(s,a)}{\lambda} \bigg)\bigg|\bigg|\pi_{\theta}(a|s)\bigg)\bigg] \\
    =& \arg \min_{\theta} \mathbb{E}_{s \sim d^{\pi_k}}\bigg[- \int_a \pi_k(a|s) \\
    &\exp \bigg( \frac{A^{\pi_k}(s,a)}{\lambda} \bigg)\log \pi_{\theta}(a|s)da \bigg]\\
    =& \arg \max_{\theta} \mathbb{E}_{s \sim d^{\pi_k}, a \sim \pi_k}\bigg[\log \pi_{\theta}(a|s)\exp \bigg( \frac{A^{\pi_k}(s,a)}{\lambda} \bigg)\bigg].
\end{align*}
This last formulation is interesting since the only dependence on $\pi_k$ is through the samples in the expectation.

According to how we estimate $A^{\pi_k}$ either AWR, AWAC or AWPO are recovered which can allow for both on-policy and off-policy data. AWAC leverages the same structure as the presented off-policy algorithms where transitions are stored in a replay buffer, the critic $Q_{\bm{\theta}_Q}$ is updated as in \eqref{eq:off_policy_critic_update} and the actor updated via
\begin{align}
\begin{split}
    \theta_{k+1} =& \arg \max_{\theta} \mathbb{E}_{(s,a) \sim \mathcal{D}}\bigg[\log \pi_{\theta}(a|s)\\
    &\exp\bigg( \frac{Q_{\bm{\theta}_Q}(s,a) - \mathbb{E}_{\bm{a} \sim \pi_{\theta}}[Q_{\bm{\theta}_Q}(s,a)]}{\lambda} \bigg)\bigg].
    \label{eq:AWAC}
\end{split}
\end{align}

On the other hand, AWR starts from the observation that the replay buffer $\mathcal{D}$ contains data collected from $N$ previous policies $\{\pi_1, \dots, \pi_N\}$. Sampling from $\mathcal{D}$ is then analagous to modeling the sampling policy as $\mu(\bm{a}|\bm{s}) = \frac{\sum_{i=1}^N w_id_{\pi_i}(s)\pi_i(a|s)}{\sum_{j=1}^N w_jd_{\pi_j}(s)}$ where the weights $\sum_i w_i = 1$ specify the probabilities of selecting each policy $\pi_i$. The following surrogate problem is then proposed
\begin{align}
    \begin{split}
        \max_{\pi} \ \ \ & \mathbb{E}_{s \sim d^{\mu}}\bigg[\mathbb{E}_{a \sim \pi(\cdot|s)}\big[A^{\pi_i}(s,a)\big]\bigg], \\
        &\mathbb{E}_{s \sim d^{\mu}}[\mathbb{D}_{\text{KL}}\big(\pi(\cdot|s) || \mu(\cdot|s)\big)] \leq \epsilon,
    \end{split}
\end{align}
where $d^{\mu} = \sum_{i=1}^N w_id_{\pi_i}$. The resolution of this problem via the Lagrangian turns similar to \eqref{eq:AWAC}. However, in \cite{peng2019advantage} the advantage is computed as
\begin{align}
\begin{split}
    A^{\pi_i}(s,a) &= \mathcal{R}^{\pi_i}_{s,a} - \bar{V}, \\
    \bar{V} &= \arg\min_{V} \sum_i w_i \mathbb{E}_{s \sim d^{\pi_i}, a \sim \pi_i}[||\mathcal{R}^{\pi_i}_{s,a} - V(s)||^2],
\end{split}
\label{eq:Adv_AWR}
\end{align}
where $\mathcal{R}^{\pi_i}_{s,a}$ is a k-steps Monte Carlo estimate of the discounted expected reward following $\pi_i$. We notice that by averaging the previous k-steps advantage estimators we recover GAE in \cite{schulman2015high}.

Finally, GePPO in \cite{queeney2021generalized} formalizes a principled way to use a replay buffer containing $N$ previous policies $\{\pi_1, \dots, \pi_N\}$. Theorem~\ref{Theo:Queeney} is proved in \cite{queeney2021generalized} and states that
\begin{align*}
    J(\pi) - J(\pi_k) \geq& \frac{1}{1-\gamma}\mathbb{E}_{(s,a) \sim \mu^{\pi_{\text{ref}}}}\bigg[\frac{\pi(a|s)}{\pi_{\text{ref}}(a|s)}A^{\pi_k}(s,a)\bigg] \\
    &- \frac{2\gamma C^{\pi, \pi_k}}{(1-\gamma)^2}\mathbb{E}_{s \sim d^{\pi_{\text{ref}}}}[\text{TV}(\pi, \pi_{\text{ref}})(s)],
\end{align*}
In the previous we distinguish $3$ policies: $\pi$ which is any future policy, $\pi_k$ which is the current policy and $\pi_{\text{ref}}$ which is the actual behavioral policy. This result can be generalized for the case in which $\pi_{\text{ref}}$ is a mixture of policies from previous iterations yielding
\begin{align}
    \begin{split}
        J(\pi) - J(\pi_k) \geq& \frac{1}{1-\gamma}\sum_{i=1}^N w_i \mathbb{E}_{s \sim d^{\pi_i}, a \sim \pi_i}\bigg[\frac{\pi(a|s)}{\pi_i(a|s)}A^{\pi_k}(s,a)\bigg] \\
        &- \frac{2\gamma C^{\pi, \pi_k}}{(1-\gamma)^2}\sum_{i=1}^N w_i \mathbb{E}_{s \sim d^{\pi_i}, a \sim \pi_i}[\text{TV}(\pi, \pi_i)(s)].
    \end{split}
    \label{eq:theo_GePPO}
\end{align}
Theorem~\ref{Theo:Queeney} justifies the use of a replay buffer storing previous experience and shows a difference in the formulation of the advantage function with respect to \eqref{eq:Adv_AWR}. The GePPO formulation shows $A^{\pi_k}$ rather than $A^{\pi_i}$ and uses an Advantage version of V-Trace \cite{espeholt2018impala} for off-policy policy evaluation. Note that both V-Trace and GePPO use importance sampling and therefore require knowledge of the behavioral policies $\{\pi_1, \dots, \pi_N\}$ used to generate the data stored in the buffer.

\subsection{Experiments training curves}

In the following we show the training curves for the experiments in Table~\ref{tab:table_online}.
\begin{figure}[h!]
    \centering
    \includegraphics[width=5cm, height=5cm]{Figures/MiniGrid-Empty-16x16-v0_AWR.pdf}
    \caption{Table~\ref{tab:table_online}: RL only}
\end{figure}

\begin{figure}[h!]
    \begin{subfigure}{\linewidth}
        \centering
        \includegraphics[width=5cm, height=5cm]{Figures/on_off_RL_from_observations_MiniGrid-Empty-16x16-v0_dataset_modified_human_expert_domain_adaptation_False.pdf}
        \caption{RL+\ref{fig:human_modified}}
    \end{subfigure}
    ~
    \begin{subfigure}{\linewidth}
        \centering
        \includegraphics[width=5cm, height=5cm]{Figures/on_off_RL_from_observations_MiniGrid-Empty-16x16-v0_dataset_modified_human_expert_domain_adaptation_True.pdf}
        \caption{RL+\ref{fig:human_modified}+DA}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{\linewidth}
        \centering
        \includegraphics[width=5cm, height=5cm]{Figures/on_off_RL_from_observations_MiniGrid-Empty-16x16-v0_dataset_rodent_intrinsic_reward_0.005.pdf}
        \caption{RL+\ref{fig:rodent}+DA (1)}
    \end{subfigure}
        ~
    \begin{subfigure}[t]{\linewidth}
        \centering
        \includegraphics[width=5cm, height=5cm]{Figures/on_off_RL_from_observations_MiniGrid-Empty-16x16-v0_dataset_rodent_intrinsic_reward_0.01.pdf}
        \caption{RL+\ref{fig:rodent}+DA (2)}
    \end{subfigure}
    \caption{Table~\ref{tab:table_online}}
\end{figure}

\end{document}


\section{Introducción}

El cómputo de la cerradura convexa, o \textit{convex hull} en ingles, es una tarea fundamental y de suma importancia para diversas aplicaciones de la ciencia y tecnología, como por ejemplo en geometría computacional, gráfica computacional y modelamiento de formas, entre otras. El convex hull es requerido frecuentemente para detección de colisiones, cálculo de interferencias, análisis de formas, reconocimiento de patrones, generación de estadísticas, sistemas de información geográfica, entre otras muchas disciplinas~\cite{barber_1996, Rourke_1985, Lin_1996}.

La convexidad es una propiedad relevante para variadas ramas de la ciencia que requieren de análisis geométrico. Una cerradura, que contiene un conjunto de puntos, es convexa si para todo par de puntos pertenecientes a esta, el segmento de línea recta que los une está completamente dentro de ella~\cite{cormen_2001}. Uno de los primeros problemas identificados en el campo de la geometría computacional es el de calcular la forma convexa más pequeña, llamada convex hull, que encierra un conjunto de puntos~\cite{preparata_1985}. Así, dado un conjunto de puntos en un espacio en $d$ dimensiones, el convex hull es su mínimo subconjunto convexo que contiene todos los puntos, tal como se muestra en la Figura \ref{fig1}.
\begin{figure}[ht!]
\centerline{\includegraphics[scale=0.6]{images/ch1.pdf}}
\caption{Cálculo del Convex Hull} \label{fig1}
\end{figure}
En otras palabras, dado un conjunto de puntos $S$, definimos que $CH(S)$ es la envolvente o cierre convexo de $S$ si i) $CH(S)$ es convexo, ii) $S$ está contenido en $CH(S)$ y iii) $CH(S)$ es el menor conjunto convexo que contiene a $S$. Para el caso de 2 dimensiones ($2D$), el convex hull corresponde al polígono convexo más pequeño que contiene a todos los puntos del conjunto. Intuitivamente podemos definir al convex hull como ``si $S$ es un conjunto finito de puntos en el plano e imaginamos rodear el conjunto con una banda elástica; cuando se suelte el elástico este asumirá la forma del convex hull''~\cite{preparata_1985}.

El cómputo del convex hull suele realizarse directamente sobre los puntos de entrada. Sin embargo, en muchos casos --- dependiendo de la instancia de entrada --- es conveniente primero realizar un filtrado de puntos, descartando elementos que son evidentemente innecesarios para el proceso final. Lo anterior da lugar a las siguientes dos etapas:
\begin{enumerate}
    \item \textbf{Filtrado}. Este paso consiste en descartar el máximo número de puntos que evidentemente no pertenecen al convex hull, por medio de una validación matemática de bajo costo computacional; esto incluye a puntos que están sobre alguno de los lados de este polígono y no son vértices del segmento de recta.
    \item \textbf{Convex hull}. Recibe como entrada los puntos que han pasado el filtro de la etapa anterior y calcula el convex hull con algún algoritmo tradicional. De esta forma, el posterior cálculo del convex hull es realizado con un número reducido de puntos, siempre y cuando el paso anterior haya logrado descartar algunos puntos, acelerando y simplificando la selección de puntos del convex hull. El cálculo en sí es realizado por algoritmos conocidos como Graham Scan~\cite{graham_1972}, Jarvis March~\cite{jarvis_1973} u otro.
\end{enumerate}
\subsection{Algoritmos para el cómputo del convex hull}
El problema de calcular el convex hull de un conjunto de puntos ha sido ampliamente estudiado en la geometría computacional. Se han propuesto muchos algoritmos teóricamente óptimos para conjuntos de puntos de baja y alta dimensión. Este fue uno de los primeros problemas en el campo a ser estudiado desde el punto de vista de la complejidad computacional, determinando que el tiempo asintótico óptimo para todo algoritmo basado en comparaciones dentro del modelo RAM\footnote{Bajo el modelo de costos RAM, es posible realizar cualquier cálculo matemático simple entre dos palabras de memoria, de tamaño $O(\log n)$ bits, en tiempo constante.}, es de es $O(n \log n)$, dónde $n$ es la cantidad de puntos de entrada~\cite{cormen_2001,akl_1984}. Además, mediante un proceso lineal, es posible reducir el problema de ordenar $n$ puntos a encontrar el convex hull en este mismo conjunto; confirmando que la envolvente convexa se puede calcular en un tiempo óptimo de $O(n \log n)$, ya que esta es la cota inferior al problema de ordenamiento.
La Tabla \ref{tabla:tiempos_algoritmos} muestra los tiempos asintóticos de los algoritmos de convex hull históricamente más conocidos para dos dimensiones. Consideremos $n$ el número de puntos en el conjunto de entrada, $h$ el número de vértices del polígono de salida. Nótese que $h \le n$, por lo que $nh \le n^2$.

\begin{table}[htbp]
\caption{Tiempos de algoritmos para el cálculo del convex hull en $2D$.}
\begin{center}
\begin{tabular}{|l|l|l|}
\hline
\textbf{Algoritmo} & \textbf{Tiempos} & \textbf{Autor} \\
\hline
Fuerza Bruta & $O( n^4 )$ & Anónimo \\ \hline
Gift Wrapping & $O ( nh ) $ & Chand \& Kapur, 1970 \\ \hline
Graham Scan & $O ( n \log n ) $ & Graham, 1972 \\ \hline
Jarvis March & $O ( nh ) $ & Jarvis, 1973 \\ \hline
QuickHull & $O ( nh ) $ & Eddy, 1977, Bykat, 1978 \\ \hline
Divide-and-Conquer & $O ( n \log n ) $ & Preparata \& Hong, 1977 \\ \hline
Monotone Chain & $O ( n \log n ) $ & Andrew, 1979 \\ \hline
Incremental & $O ( n \log n ) $ & Kallay, 1984 \\ \hline
Marriage-before-Conquest & $O ( n \log n ) $ & Kirkpatrick \& Seidel, 1986 \\ \hline
\end{tabular}
\label{tabla:tiempos_algoritmos}
\end{center}
\end{table}

\subsection{Algoritmo Heaphull}
El trabajo de Ferrada et al.~\cite{ferrada_2019} está basado en el desarrollo de una técnica de optimización que reduce el costo computacional para construir el convex hull. Este método preprocesa el conjunto de puntos de entrada, filtrando todos los puntos dentro de un polígono de ocho vértices en tiempo $O(n)$ y devuelve un conjunto reducido de puntos candidatos para encontrar el convex hull, los cuales se encuentran semi ordenados y distribuidos en cuatro colas de prioridad (una para cada cuadrante).

Los resultados experimentales muestran que en el peor de los casos (cuando todos los puntos se encuentran sobre una circunferencia) un pequeño desplazamiento radial aleatorio de los puntos hace de este método el más rápido. Además, al aumentar la magnitud de este desplazamiento, el rendimiento del método propuesto escala a un ritmo más rápido que las otras implementaciones testeadas. En términos de eficiencia de memoria, esta implementación logra usar de $3\times$ a $6\times$ menos memoria que otros métodos.

Este algoritmo es esencialmente una técnica de filtrado que se conforma de cuatro etapas relevantes, (1) construcción de un octágono con los puntos extremos, (2) filtrado de puntos y agrupación en colas de prioridad, (3) ordenamiento parcial y (4) cálculo del convex hull. Las primeras tres etapas toman tiempo $O(n)$, mientras que la etapa de cálculo del convex hull toma $O (n' \log n' )$, donde $n'$ es el tamaño del conjunto filtrado $P' \subseteq P$. El Algoritmo~\ref{alg_heaphull} muestra la idea de heaphull.

\begin{center}
\begin{algorithm}[H]
\small
\caption{\textbf{heaphull}, Ferrada et al.~\cite{ferrada_2019}. To compute the convex hull $CH[1..h]$ in $2D$}\label{euclid}
\begin{algorithmic}[1]
\Require $n$ points in $2D$, stored in an array $P[1..n]$.
\Ensure an array $CH[1..h]$ with the $h$ points of the hull in counter-clockwise order.
\Procedure{heaphull}{$P,n$}
\State $E\gets\textsc{findExtremes}($P$, $n$)${\scriptsize\Comment{{\color{blue}find extreme points ($O(n)$ time)}}}
\State $n_1=n_2=n_3=n_4=0${\scriptsize\Comment{{\color{blue}counter for the 4 priority queues $Q_i$}}}
\For{\texttt{$j=1$ to $n$}}
	\If{$P[j]$ is outside the convex octagon $CP(E)$}
		\State $i\gets\textsc{findQueue}(P,n,j$) {\scriptsize\Comment{{\color{blue}find $Q_i$ for $P[j]$}}}
		\State store $j$ in the priority queue $Q_i$
		\State $n_i\gets n_i+1$
	\EndIf
\EndFor
\State $CH\gets \emptyset$
\For{\texttt{$i=1$ to $4$}}
    \scriptsize\Comment{{\color{blue}add the CH for points of $Q_i$ ($O(n_i\log n_i)$ time)}}
	\State $CH\gets CH\cup\textsc{hull}(Q_i, n_i)$
\EndFor
\State \textbf{return} $CH$
\EndProcedure
\end{algorithmic}
\label{alg_heaphull}
\end{algorithm}
\end{center}

\subsection{Paralelismo GPU en el cómputo del Convex Hull}

Las GPU (Graphic Process Unit) se han utilizado para acelerar el rendimiento de diversas aplicaciones, como la simulación de partículas, el modelado molecular y el procesamiento de imágenes. La arquitectura GPU esta basada en el paralelismo masivo, proporcionando cientos de cores (elementos de procesamiento), donde cada core es un procesador pipeline multi-etapa. Los cores se agrupan para generar un multiprocesador (SM) de caracter SIMD al nivel de grupos de 32 threads. Cada SM tiene su propio conjunto de registros y memoria local, y los cores dentro del mismo SM tienen memoria compartida limitada.

Este trabajo propone una adaptacion del algoritmo heaphull al modelo de programacion GPU, para acelerar su rendimiento y aprovechar el paralelismo de datos disponible en la etapa de filtrado. En particular se busca responder las siguientes preguntas:
\begin{itemize}
    \item ¿Cuánto mejora el desempeño de la etapa de filtrado al paralelizar en GPU el algoritmo heaphull?
    \item Al paralelizar heaphull ¿Se logra que el algoritmo sea más eficiente que otros que ya usan paralelismo?
\end{itemize}
El resto de este artículo se divide en 5 secciones. En la primera sección se describe el trabajo relacionado, considerando especialmente los algoritmos que realizan un esfuerzo por paralelizar la etapa de filtrado. Para esto se clasificaron en algoritmos secuenciales tradicionales, algoritmos secuenciales que explotan el filtrado previo, algoritmos paralelos de filtrado en CPU y algoritmos paralelos de filtrado en GPU. En la segunda sección se define en detalle la metodología de investigación implementada. En la sección 3 se presenta el diseño del algoritmo paralelo en CPU y en GPU. En la sección 4 se presentan los resultados de rendimiento (speedup, eficiencia y comparación entre CPU y GPU). Finalmente en la sección 5 se desarrollan las conclusiones.

\section{Trabajo Relacionado}
Existen varios métodos prácticos conocidos e implementaciones de software robustas que abordan el cómputo del convex hull de manera secuencial. Desde 1970 se han desarrollado varios algoritmos clásicos~\cite{jiayin_2015}, como Graham scan~\cite{graham_1972}, Gift wrapping~\cite{jarvis_1973}, Incremental method~\cite{kallay_1984}, Divide-and-Conquer~\cite{preparata_1977}, Monotone chain~\cite{andrew_1979}, and Quick-Hull~\cite{barber_1996}.

\subsection{Convex hull con etapa de filtrado}
Realizar un filtrado eficiente del conjunto de puntos puede tener un gran impacto en el tiempo final de cálculo. Vyšniauskaitė~\cite{laura_2006} en 2006 propuso una idea de la filtración a priori de puntos, presentado dos nuevos algoritmos basados en la búsqueda de puntos extremos y en la subdivisión del conjunto de puntos en arreglos más pequeños.

También Sharif~\cite{sharif_2011} el 2011 propone un método híbrido para calcular el convex hull, basado en dos algoritmos ya existentes, es decir, Quickhull y GrahamScan, intentando eliminar las deficiencias en las dos técnicas mencionadas anteriormente usando el primero para realizar un filtrado y realizando el cálculo final con el segundo de estos.

\subsection{Convex hull paralelo}
Para mejorar la eficiencia computacional del convex hull se han hecho algunas contribuciones valiosas rediseñando e implementando algoritmos de convex hull secuenciales en paralelo, explotando el potencial de cálculo masivo de las GPU. La mayoría de estas implementaciones se han desarrollado teniendo como base el algoritmo QuickHull~\cite{jiayin_2015}.
\subsubsection{Paralelismo en CPU}
Como se muestra la publicación de Chen y otros~\cite{chen_1992} el convex hull puede ser calculado de forma paralela en el modelo PRAM (Máquina de acceso aleatorio paralelo), siendo el algoritmo directamente implementable en CPUs multicore.

En 1987 Goodrich~\cite{goodrich_1987} diseñó un algoritmo paralelo para encontrar el convex hull que se ejecuta en tiempo $O(log ~n)$ utilizando  $O (n / log ~n)$ procesadores en el modelo computacional CREW (concurrent read and exclusive write) PRAM, que es óptimo. Una de las técnicas que utilizan para alcanzar estos límites óptimos es el uso de una estructura de datos paralela que llaman ``hull tree''.

Miller~\cite{miller_1988} en 1988 presentó algoritmos paralelos para identificar los puntos extremos del convex hull, concentrándose en el desarrollo de estos en tiempo polilogarítmico para una variedad de máquinas paralelas y analizándolos usando notación O. Asimismo, Liu~\cite{liu_2015} en 2015 hace una revisión completa de los algoritmos paralelos existentes a la fecha, categorizándolos según su arquitectura y considerando sus tiempos de ejecución en notación O y el número de procesadores necesarios para esto.

Berkman~\cite{berkman_1996} en 1996 desarrolló un algoritmo paralelo para encontrar el convex hull en tiempo $O (log ~log ~n)$ usando $n/log ~log ~n$ procesadores en una arquitectura PRAM CRCW (Concurrent read concurrent write) común. Para romper la barrera de tiempo $\Omega (log ~n / log ~log ~n)$ requerida para generar el convex hull, introduce una estructura de datos, un árbol balanceado de altura doblemente logarítmica y usa este para representarlo. El algoritmo demuestra el poder del ``paradigma doblemente logarítmico de divide y vencerás''.

Nakagawa~\cite{nakagawa_2009} presentó una implementación de algoritmo paralelo simple para calcular el convex hull y evaluar el rendimiento en los procesadores de cuatro núcleos duales, logrando un factor de aceleración de aproximadamente 7 utilizando 8 procesadores. Como el factor de aceleración de más de 8 no es posible, su implementación paralela para calcular el convex hull es casi óptima.

Waghmare~\cite{waghmare_2010} presenta un algoritmo de convex hull usando clustering K-means, en el que los puntos en 2D se agrupan en clusters diferentes y luego se calculan los convex hull para cada uno de estos. El algoritmo se implementa en modo MPI, OpenMP e híbrido. Los resultados indican que el enfoque híbrido supera al enfoque MPI y OpenMP por separado.

\subsubsection{Paralelismo en GPU}
El trabajo relacionado en computación paralela para el cálculo del convex hull en GPU muestra que se puede alcanzar mejor velocidad y rendimiento, como se evidencia en el trabajo de Tang et al.~\cite{tang_2012}, y Jiayin et al~\cite{jiayin_2015}. Estos algoritmos basados en GPU trabajan sobre la optimización del filtrado de puntos, realizando la etapa posterior de cálculo en forma secuencial.

Srungarapu et al.~\cite{srung_2011} propuso en 2011 una implementación optimizada en GPU para encontrar el convex hull para conjuntos de puntos de dos dimensiones. Su implementación trata de minimizar el impacto de los patrones de acceso irregular a los datos, logrando un speedup de hasta 14 veces sobre la implementación estándar secuencial en CPU.

En 2011 Jurkiewicz y Danilewski~\cite{jurki_2011} presentaron una propuesta de implementación de algoritmo de cálculo del convex hull implementada en su totalidad en GPU, sin filtrado y basado en 3 algoritmos secuenciales distintos, con pruebas experimentales con conjuntos de entrada de $10^7$ puntos, obteniendo resultados en reducción de tiempo de hasta $100$ veces para el máximo de puntos al compararse con una librería que utiliza algoritmos similares, pero que su desempeño se encuentra muy por debajo de la media.

Mei G.~\cite{mei_2016} en 2016 propuso una versión paralela en CUDA (GPU) para el cálculo de convex hull para dos~\cite{mei_2017} y tres~\cite {mei_2015} dimensiones, alcanzando aceleraciones de aproximadamente 4x en promedio y $5x \sim 6x$ en los mejores casos, la cual llama CudaChain. Realizando solo una mejora de filtrado, afirma que el 95\% de puntos de entrada se pueden descartar en la mayoría de las pruebas experimentales. Anterior a esto, usó GPUs para mejorar la fase de filtrado en la búsqueda de puntos extremos sobre datos sintéticos y con el uso de la biblioteca qhull (implementación de Quickhull en C++), la idea básica era descartar los puntos que se ubican dentro de un polígono convexo formado por 16 extremos~\cite{mei_2014}.

Por otro lado, Stein~\cite{stein_2012} en 2012 hizo una versión paralela para el convex hull en 3 dimensiones que se basa en el enfoque QuickHull y comienza construyendo un tetraedro inicial utilizando cuatro puntos extremos, descarta los puntos internos y distribuye los puntos externos a las cuatro caras y continúa luego iterativamente. En su trabajo, Stein afirma que su implementación superó al Qhull basado en CPU en 30 veces para 10 millones de puntos y 40 veces para 20 millones de puntos.

En 2020 Masnadi~\cite{masnadi_2020} presenta ConcurrentHull, una nueva técnica basada en el filtrado de puntos para conjuntos de datos 2D y 3D. Su implementación, que es una combinación de filtrado, divide y vencerás, y computación paralela, permite ser empleada en un entorno de computación distribuida. Su algoritmo tiene una versión para CPU y otra para GPU (CUDA). Los resultados muestran que tiene una ganancia de rendimiento con grandes tamaños de datos de entrada y tiene la ventaja de poder manejar grandes conjuntos de datos.

Si bien las soluciones propuestas informan un speedup importante para las condiciones establecidas para las pruebas, nuestro trabajo busca no solo mejorar los tiempos de ejecución si no ponemos especial atención en el porcentaje de filtrado a obtener en la fase paralela, lo cual tiene el mayor impacto en el tiempo de cálculo posterior. Nuestro trabajo compara las dos soluciones más actuales (CudaChain y ConcurrentHull) con la nuestra (GPU Heaphull), comparando la solución en el peor caso y en el caso promedio.

\section{Paralelización del Algoritmo Heaphull}
El primer paso de la versión paralela consiste en encontrar los cuatro puntos extremos en las direcciones norte, sur, este y oeste del conjunto de puntos y formar el cuadrilátero que definen. Tras este paso, se tienen puntos distribuidos en cuatro regiones, formadas por los puntos fuera del cuadrilátero. En cada uno de ellas existe un punto que es el más cercano a la esquina correspondiente. Tras la obtención de estos cuatro puntos adicionales, se obtiene un octágono que divide los puntos restantes externos en ocho regiones que se pueden tratar individualmente siguiendo la misma regla.

La paralelización en GPU se realiza acelerando la reducción de mínimos y máximos para formar el polígono de filtrado, usando el enfoque de Shuffle Warp Reduce, el cual permite que conjuntos de 32 threads, llamados \textit{warps}, puedan comunicar sus registros de forma eficiente.

A diferencia de la versión tradicional de heaphull, que implementa la búsqueda de los $8$ puntos en cada paso del bucle, en la versión GPU se ejecuta un \textit{kernel} de reducción para encontrar los primeros $4$ puntos extremos, seguido de un segundo \textit{kernel} de reducción para encontrar los $4$ puntos de las esquinas. La razón de por qué se deben ejecutar dos \textit{kernels} en sucesión es debido a que la búsqueda de los puntos secundarios requiere de los primeros $4$ puntos.

En el segundo paso de la búsqueda de los puntos de esquina, se utilizan distancias de Manhattan, es decir, para cada esquina $c_i, 1 \leq i \leq 4$, se encuentra el punto en $P$ que minimiza la suma de las distancias vertical y horizontal a la esquina correspondiente. Si bien la distancia de Manhattan es una métrica que puede entregar puntos mas distantes a la esquina que la distancia Euclidiana, esto tiende a ocurrir con baja frecuencia, y cuando ocurre, se puede ver que el punto elegido todavía se encuentra muy cerca del más cercano por distancia Euclidiana~\cite{ferrada_2019}. La ventaja de usar la distancia Manhattan es que no requiere del cómputo de raíces cuadradas, lo cual lo convierte en una métrica mas rápida de calcular que la distancia euclidiana.

En el paso siguiente se realiza el descarte de los puntos que deben ser filtrados y se crean las colas. Para la creación de las 4 colas se usa un arreglo donde se asigna el número de cola ($1$, $2$, $3$ o $4$) a la posición donde se encuentra un punto que si pertenece a la cola respectiva y un $0$ cuando ha sido descartado. Luego, en la CPU se crean las 4 colas con los respectivos índices guardados para ser entregadas a la sección correspondiente del algoritmo secuencial heaphull.

El Algoritmo~\ref{alg_filterGPU} muestra los pasos computacionales del proceso.
\begin{center}
\begin{algorithm}
\small
\caption{\textbf{GPUfilter}. Para calcular el octágono de filtrado de $P[1..n]$ en $2D$}
\begin{algorithmic}[1]
\Require $n$ valores de punto flotante, en $2D$, almacenados en un arreglo $P[1..n]$.
\Ensure un arreglo $Q[1..n']$ con los índices de la cola a la que pertenece cada uno de los puntos.
\Procedure{GPUfilter}{$P,n$}
\State$E\gets\textsc{findExtremes}($P$, $n$)${\scriptsize\Comment{{\color{blue} extremos en tiempo $O(n)$}}}
\State Sea $Q[]$ un arreglo de n enteros
\For{\texttt{$j=1$ to $n$}}
	\If{$P[j]$ está fuera del octágono $CP(E)$}
	    \State{\scriptsize\Comment{{\color{blue} encuentra $Q_i$ para $P[j]$ (tiempo $O(1)$)}}}
		\State $i\gets\textsc{findQueue}(P,n,j$)
		\State $Q[j] = i$

	\EndIf
\EndFor
\State \textbf{return} $Q$
\EndProcedure
\end{algorithmic}
\label{alg_filterGPU}
\end{algorithm}
\end{center}

Como referencia, también se implementó una versión paralela de heaphull en CPU multi-core. Se probaron dos versiones; una con el patrón \textit{reduce} de OpenMP para mínimos y máximos, y otra usando la técnica de \textit{divide and conquer}. Comparaciones de rendimiento mostraron que esta última logra mayor aceleración, de hasta $3x$ para el filtrado, siendo la escogida.

\section{Evaluación Experimental}
A fin de comparar experimentalmente el desempeño de nuestra solución paralela en GPU, incluimos como baseline las implementaciones de: heaphull tradicional~\cite{ferrada_2019}, CudaChain~\cite{mei_2016} y ConcurrentHull~\cite{masnadi_2020}; este último permite el manejo de una mayor cantidad de puntos y consigue mejores resultados en esas condiciones.
Para las pruebas se utilizó un computador con Procesador Intel\textregistered Core\texttrademark i5-8300H CPU (2.30GHz × 8), 8 GB de memoria RAM, tarjeta gráfica NVIDIA GeForce GTX $1050$ Ti, CUDA Version 11.6 en Ubuntu 20.04 LTS.

\subsection{Tiempos de Ejecución}
Para cada implementación se crean tests puntos de tamaños desde $10^4$ hasta $10^8$. Cada test se ejecuta 100 veces y se calcula el tiempo promedio de ejecución en cada algoritmo.

Los tiempos promedio de búsqueda de puntos extremos para las pruebas experimentales en las versiones de heaphull secuencial, paralelo en CPU y paralelo en GPU se pueden ver en la Tabla~\ref{tabla:tiempos_extremos} y en el gráfico de la Figura~\ref{fig2}. Los resultados de heaphull paralelo CPU solo son una referencia; los gráficos futuros solo usan la versión secuencial (Heaphull CPU).

\begin{table}[htbp]
\caption{Tiempo promedio de búsqueda de puntos extremos en ms.}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
\textbf{Puntos} & \textbf{Tiempo CPU} & \textbf{Tiempo CPU Paralelo} & \textbf{Tiempo GPU} \\
\hline
$10^4$ & $0.0460$ & $0.1000$ & $0.0433$ \\ \hline
$10^5$ & $0.4600$ & $0.4000$ & $0.0814$ \\ \hline
$10^6$ & $5.0120$ & $2.4000$ & $0.8505$ \\ \hline
$10^7$ & $56.8020$ & $18.0000$ & $6.9785$ \\ \hline
$10^8$ & $793.7420$ & $289.1300$ & $67.0598$ \\ \hline
\end{tabular}
\label{tabla:tiempos_extremos}
\end{center}
\end{table}

\begin{figure}[htbp]
\centerline{\includegraphics[width=\columnwidth]{images/filtrado.pdf}}
\caption{Tiempos promedio de filtrado para las versiones secuencial, paralela en CPU y paralela en GPU.} \label{fig2}
\end{figure}

\subsubsection{Desempeño caso Promedio (Dist. Normal)}

En la Tabla~\ref{tabla:tiempos_ch_gpu} se pueden ver los tiempos de cálculo de la versión secuencial de heaphull, CudaChain, ConcurrentHull y nuestra implementación en GPU, en milisegundos para $10^4$ a $10^8$ puntos. Nuestra solución es la que obtiene mejores tiempos en los casos con $10^5$ puntos o más. En la Figura~\ref{fig3} podemos ver el gráfico de los tiempos en milisegundos para el cálculo de convex hull en las cuatro implementaciones anteriores.

\begin{table}[htbp]
\begin{center}
\caption{Tiempo promedio cómputo del convex hull en el caso promedio.}
\begin{tabular}{|c|c|c|c|c|}
\hline
\textbf{Puntos} & \textbf{Heaphull} & \textbf{CudaChain} & \textbf{ConcurrentHull} & \textbf{GPU HH} \\
\hline
$10^4$ & $0.1201$ & $6.4803$ & $882.0300$ & $0.1575$ \\ \hline
$10^5$ & $1.2014$ & $12.0930$ & $902.0700$ & $0.9940$ \\ \hline
$10^6$ & $12.0824$ & $24.1600$ & $956.3200$ & $7.0058$ \\ \hline
$10^7$ & $139.9730$ & $224.2051$ & $1211.6300$ & $35.1825$ \\ \hline
$10^8$ & $2054.5416$ & $1517.1612$ & $2201.5500$ & $464.8182$ \\ \hline
\end{tabular}
\label{tabla:tiempos_ch_gpu}
\end{center}
\end{table}

\begin{figure}[htbp]
\centerline{\includegraphics[width=\columnwidth]{images/tiempos_normal.pdf}}
\caption{Tiempos para el cómputo del convex hull en el caso promedio} \label{fig3}
\end{figure}

El Speedup obtenido se puede ver en la Tabla~\ref{tabla:speedup} y en la Figura~\ref{fig5}, pudiendo identificarse que  al trabajar con $10^8$ puntos se logra una disminución de tiempo de hasta $4.4\times$ respecto al algoritmo secuencial de heaphull y de $3.2\times$ respecto a CudaChain, el cual ya logra un speedup de hasta $6\times$ comparado con la librería de Quickhull (Qhull en C++). No se considera en el gráfico el speedup frente a ConcurrentHull de menos de $10^7$ puntos para poder mantener una escala apreciable.

\begin{table}[htbp]
\begin{center}
\caption{Speedup solución desarrollada sobre heaphull, CudaChain y ConcurrentHull, caso promedio.}
\begin{tabular}{|c|c|c|c|}
\hline
\textbf{Puntos} & \textbf{Heaphull} & \textbf{CudaChain} & \textbf{ConcurrentHull} \\
\hline
$10^4$ & $0.7625$ & $41.1448$ & $5600,1905$ \\ \hline
$10^5$ & $1.2087$ & $12.1660$ & $907.5151$ \\ \hline
$10^6$ & $1.7246$ & $3.4486$  & $136.5040$ \\ \hline
$10^7$ & $3.9785$ & $3.5303$  & $34.4384$ \\ \hline
$10^8$ & $4.4201$ & $3.2640$  & $4.7364$ \\ \hline
\end{tabular}
\label{tabla:speedup}
\end{center}
\end{table}

\begin{figure}[htbp]
\centerline{\includegraphics[width=\columnwidth]{images/speedup_GPU_HH.pdf}}
\caption{Speedup implementación desarrollada sobre heaphull, CudaChain y ConcurrentHull, caso promedio.} \label{fig5}
\end{figure}

En relación al porcentaje de filtrado de puntos, heaphull secuencial y la versión en GPU de este logran porcentajes para puntos provenientes de una distribución normal sobre el 99,99\% en promedio (igual porcentaje para mismos puntos de entrada). En el caso de CudaChain logra un filtrado de 98,86\% en promedio en el mejor caso ($10^8$ puntos de entrada) y en el peor caso para este ($10^4$ puntos de entrada) llega a un 91,6\%, mientras que heaphull en sus dos versiones logra un 99,87\% para esa misma cantidad de puntos.

\subsubsection{Desempeño Peor Caso (Circunferencia)}

Para el peor caso (los puntos sobre la circunferencia) con $10^7$ puntos nuestra implementación paralela en GPU demora en encontrar el convex hull $6104.58$ ms, sin embargo con 2\% de distorsión el tiempo disminueye a $4973.3$ ms; el primero representa un speedup de $0.96$ respecto del algoritmo secuencial de heaphull, sin embargo en el segundo caso se obtiene un speedup de $1.41$. Para $10^8$ puntos el speedup es de $1,04$ en el peor caso. En la Tabla~\ref{tabla:tiempos_circunsferencia} se encuentran los tiempos de ejecución en el peor caso. En la figura~\ref{fig6} se puede ver el gráfico con los tiempos de para cada versión.

\begin{table}[htbp]
\begin{center}
\caption{Tiempos para el cómputo del convex hull en el peor caso.}
\begin{tabular}{|c|c|c|c|c|}
\hline
\textbf{Puntos} & \textbf{Heaphull} & \textbf{CudaChain} & \textbf{ConcurrentHull} & \textbf{GPU HH} \\
\hline
$10^4$ & $1.0153$ & $86.7300$ & $90.0223$ & $1.6400$ \\ \hline
$10^5$ & $15.0671$ & $112.5680$ & $96.3250$ & $23.6950$ \\ \hline
$10^6$ & $241.7010$ & $445.2650$  & $1256.3600$ & $427.0750$ \\ \hline
$10^7$ & $5861.8900$ & $6459.5251$  & $6201.5522$ & $6104.5800$ \\ \hline
$10^8$ & $95278.10$ & $121426.00$  & $123045.20$ & $91559.80$ \\ \hline
\end{tabular}
\label{tabla:tiempos_circunsferencia}
\end{center}
\end{table}

\begin{figure}[htbp]
\centerline{\includegraphics[width=\columnwidth]{images/tiempos_circunferencia.pdf}}
\caption{Tiempos de cómputo del convex hull en el peor caso} \label{fig6}
\end{figure}

Con una distorsión de sólo 2\%, los porcentajes de filtrado para $10^8$ puntos se mantienen en 10.50\% en promedio, con tiempos de ejecución del algoritmo completo similares al caso promedio. La versión codificada del algoritmo tiene un comportamiento con un speedup no relevante cuando se trata de cantidades de puntos inferiores a los $10^7$ puntos. Además, cuando se trata del peor caso no hay una mejora significativa respecto a la versión secuencial.

En general, los mejores tiempos se logran al trabajar con $10^7$ y $10^8$ puntos, en estos casos el speedup es mayor a 4x.

\section{Discusión y Conclusiones}
Nuestro trabajo presenta un implementación paralela en GPU basada en el algoritmo heaphull que mejora el tiempo de filtrado y el tiempo total de cálculo del convex hull, manteniendo las mismas garantías asintóticas óptimas del algoritmo tradicional para $2D$. El algoritmo se basa en realizar un filtrado con un octágono que permite descartar todos los puntos que se encuentran en su interior y lograr en el caso promedio filtrar el $99,99\%$ de los puntos del conjunto de entrada, quedando una cantidad muy reducida de puntos para el cálculo del convex hull; conservando así la eficiente del filtrado del algoritmo heaphull original.

El mayor impacto en desempeño se logra al trabajar con grandes volúmenes de puntos de entrada (sobre $10^7$), siendo un aporte importante en esta arista para el cálculo del convex hull, lo que es característico en la presente era del Big Data.

En comparación con otras implementaciones en GPU se obtienen mejores tiempos totales, considerando que se está trabajando sobre un algoritmo que ya resulta ser el más rápido secuencial y en $2D$, esto es un importante aporte al desarrollo futuro de implementaciones que puedan requerir el cálculo en tiempo real de grandes volúmenes de puntos.

Dentro de las limitaciones de la solución propuesta se encuentra el uso de memoria de GPU, el cual se ve aumentado por la necesidad de crear arreglos para almacenar los índices que son resultado de las reducciones para poder conservar el arreglo original de puntos de entrada.

Como trabajo futuro tenemos previsto ofrecer un nuevo diseño e implementación del convex hull con algoritmos que sean paralelizables (Jarvis March, Quickhull, Kirkpatrick–Seidel, Chan's algorithm) para comprobar si alguno de ellos, a pesar de aumentar la complejidad del algoritmo, puede lograr mejores tiempos aún al aprovechar el paralelismo de la GPU. Por otra parte, la implementación de la solución para más dimensiones puede ser un importante aporte, por sobre todo en $3D$, pero se debe considerar que el tamaño requerido de memoria crece en potencias según la cantidad de dimensiones, lo cual puede dificultar el cálculo de grandes volúmenes de datos en la GPU.

También en el trabajo futuro está, siguiendo el modelo de heaphull, el uso eficiente de memoria, lo cual podría ser un importante aporte de la mano de grandes volúmenes de datos, logrando data más compacta y facilitando el cálculo de mayores números de datos.