\section{Introduction}
Diffusion equations with fractional-order derivatives in time (which is called as time-fractional diffusion equations) have been introduced in Physics by Nigmatullin \cite{Nigmatullin_86} to describe super slow diffusion process in a porous medium with the structure type of fractal geometry (Koch's tree). Then, by the probabilistic point of view, in the paper \cite{Metzler_00}, Metzler and Klafter have pointed out that a time-fractional diffusion equation generates a non-Markovian diffusion process with a long memory. After that, Roman and Alemany \cite{Roman_94} have considered continuous-time random walks on fractals and observed that the average probability density of random walks on fractals obeys a diffusion equation with a fractional time derivative asymptotically. Another context where such systems appear is the modelling of evolution processes in materials with memory, see e.g., \cite{Pruss, Caputo_99}.

The existence of solutions to time-fractional partial differential equations has been studied by many authors. In \cite{Kochubei_04}, using Fourier transform, the authors have obtained a fundamental solution for time-fractional elliptic equations with smooth coefficient. Combining the Galerkin method and the Yoshida approximation sequence, in \cite{Zacher_09},  the author has proposed a way to prove existence of certain weak solutions to abstract evolutionary integro-differential equations in Hilbert spaces. By virtue of the operator theory in functional analysis and  the eigenfunction expansion method for symmetry elliptic operators, in \cite{Yamamoto_11}, Sakamoto and Yamamoto have proved  the existence and uniqueness of the weak solution for a fractional diffusion-wave equation. Based on a definition of the Caputo derivative on a finite interval in fractional Sobolev spaces, Gorenflo, Luchko and Yamamoto  \cite{Gorenflo_15} have investigated solutions (in the distribution sense) to time-fractional diffusion equations from the operator theoretic viewpoint. In the recent works, by combining a special approximate solution with the Galerkin approximation method, Kubica and Yamamoto \cite{Adam} have proved the existence of weak solutions to a large class of fractional diffusion equations with coefficients are dependent on spatial and time variables, by a classical variational approach,  K.V. Bockstal \cite{Bockstal, Bockstal_1} has established the existence of a unique weak solution to a certain class of fractional diffusion equations with Caputo derivative.

However, to the best of our knowledge, the development of this theory is still in its infancy and requires further researches. The main difficulty which one have to face is the meaning of the initial condition of solutions and the correctness of the formulation of weak solutions. These are due to the fact that the solutions of time-fractional partial differential equations with with time-dependent coefficients are not in standard Sobolev spaces and so classical embedding theorems in functional analysis do not apply. In this paper, we only focus on an initial-boundary value problem with the zero initial condition for time-fractional elliptic equations having time dependent variable coefficients which can not apply Fourier transformation to solve them. To overcome the aforementioned obstacles, our strategy is to use the Galerkin method, a basic inequality for the fractional derivative of convex Lyapunov candidate functions, the Yoshida approximation sequence and the weak compactness argument.

The paper is organized as follows. In section 2, we recall some preliminary results on fractional calculus. Then, we give the setting of the problem and propose a clear definition of a weak solution to a time fractional elliptic equation with time dependent variable coefficients. The main result of the paper is Theorem 3.4 on the existence and uniqueness of weak solutions introduced in section 3.

To conclude the introduction, we will introduce some notations used throughout the rest of the article. Denote by $\mathbb{N}$ the set of natural numbers and by $\mathbb{R}$ the set of real numbers. For any natural $d\in \mathbb{N}$, let $\mathbb{R}^d$ be the $d$-dimensional Euclidean space. For a open subset $\Omega$ of $\mathbb{R}^d$, let $C^\infty_c(\Omega)$ denote the space of infinitely differentiable functions $f: \Omega\rightarrow \mathbb{R}$ with the compact support in $\Omega$, $L^p(\Omega)$, $p\in \mathbb{N}$, be the set of measurable functions $f:\Omega\to\mathbb{R}$ such that
\[
\int_\Omega |f(x)|^pdx<\infty,
\]
$H^1(\Omega)$ be the Sobolev space containing all locally integrable functions $f:\Omega\to \mathbb{R}$ such that $f$ and its weak derivatives belong to $L^2(\Omega)$, $H^1_0(\Omega)$ be the closure of $C^\infty_c(\Omega)$ in $H^1(\Omega)$, and $H^{-1}(\Omega)$ be the dual space of $H^1_0(\Omega)$. Fix $T>0$, denote by $W^1_1([0,T];\mathbb{R})$ the space of functions $f:[0,T];\mathbb{R}$ such that $f$ and its weak derivative belong to the space $L^1([0,T];\mathbb{R})$ and $W^1_2([0,T];H^{-1}(\Omega))$ by the space of functions $f:[0,T]\rightarrow H^{-1}(\Omega)$ such that $f$ and the weak derivative belong to the space $L^2([0,T];H^{-1}(\Omega))$.
\section{Fractional calculus}
We briefly recall an abstract framework of fractional calculus.

Let $\alpha\in (0,1]$, $[0,T]\subset \mathbb{R}$ and $x:[0,T]\rightarrow \mathbb{R}$ be a measurable function such that $\int_0^T|x(\tau)|\;d\tau<\infty$. The right-handed \emph{Riemann--Liouville integral operator of order $\alpha$} is defined by
\[
(I_{0+}^{\alpha}x)(t):=\frac{1}{\Gamma(\alpha)}\int_0^t(t-\tau)^{\alpha-1}x(\tau)\;d\tau,
\]
where $\Gamma(\cdot)$ is the Gamma function. The left-handed \emph{Riemann--Liouville integral operator of order $\alpha$} is defined by
\[
(I_{T^-}^{\alpha}x)(t):=\frac{1}{\Gamma(\alpha)}\int_t^T(\tau-t)^{\alpha-1}x(\tau)\;d\tau.
\]
We have the following result on the relation between the right-handed and left-handed Riemann--Liouville integral operators.
\begin{theorem}[The Hardy--Littlewood form of fractional integration by parts]
If $p>1$, $q>1$, $0<\alpha<1$, $\frac{1}{p}+\frac{1}{q}-1\leq \alpha$, and \[ f\in L^p([0,T];\mathbb{R}),\;g\in L^q([0,T];\mathbb{R}),
\]
then
\[
\int_0^T (I^\alpha_{0+} f(t))g(t)dt=\int_0^T f(t)I^\alpha_{T^-}g(t)dt.
\]
\end{theorem}
\begin{proof}
See \cite{LY_37}.
\end{proof}
The right-handed \emph{Riemann--Liouville fractional derivative} $^{RL\!}D_{a+}^\alpha x$ of $x$ on $[0,T]$ is defined by
\[
^{RL\!}D_{0+}^\alpha x(t):=(DI_{0+}^{1-\alpha}x)(t)\;\text{for almost}\;t\in [0,T],
\]
where $D=\frac{d}{dt}$ is the usual derivative. The left-handed \emph{Riemann--Liouville fractional derivative} $^{RL\!}D_{T^-}^\alpha x$ of $x$ on $[0,T]$ is defined by
\[
^{RL\!}D_{T^-}^\alpha x(t):=-(DI_{T^-}^{1-\alpha}x)(t)\;\text{for almost}\;t\in [0,T].
\]
\begin{theorem}\label{integration_by_parts}
The formula
\[
\int_0^T f(t)^{RL}D^\alpha_{0+}g(t)dt=\int_0^T g(t)^{RL}D^\alpha_{T^-}f(t)dt
\]
is valid for $0<\alpha<1$, $f\in I^\alpha_{T^-}(L^p)$, $g\in I^\alpha_{0+}(L^q)$ and $\frac{1}{p}+\frac{1}{q}-1\leq \alpha$.
\end{theorem}
\begin{proof}
See \cite[Corollary 2, pp. 46]{Samko}.
\end{proof}
The right-handed \emph{Caputo fractional derivative} of $x$ on $[0,T]$ is defined by
\[
^{\!C}D^\alpha_{0+}x(t)=^{RL\!}D_{a+}^\alpha (x(t)-x(0))\quad\text{for almost}\;t\in [0,T]
\]
and the left-handed \emph{Caputo fractional derivative} of $x$ on $[0,T]$ is defined by
\[
^{\!C}D^\alpha_{T^-}x(t)=^{RL\!}D_{T^-}^\alpha (x(t)-x(T))\quad\text{for almost}\;t\in [0,T].
\]
We have a sufficient condition for the existence of fractional derivative.
\begin{theorem}
Let $f\in AC([0,T];\mathbb{R})$, $\alpha\in (0,1)$, then $^{RL}D^\alpha_{0+}f$ and $^{RL}D^\alpha_{T^-}f$ exist almost everywhere. Moreover, $^{RL}D^\alpha_{0+}f, ^{RL}D^\alpha_{T^-}f\in L^r([0,T];\mathbb{R})$, $1\leq r\leq \frac{1}{\alpha}$, and
\[
^{RL}D^\alpha_{0+}f(t)=\frac{1}{\Gamma(1-\alpha)}\Big[\frac{f(0)}{t^\alpha}+\int_0^t \frac{f'(\tau)}{(t-\tau)^\alpha}d\tau\Big],
\]
\[
^{RL}D^\alpha_{T^-}f(t)=\frac{1}{\Gamma(1-\alpha)}\Big[\frac{f(T)}{(T-t)^\alpha}-\int_t^T \frac{f'(\tau)}{(\tau-t)^\alpha}d\tau\Big].
\]
\end{theorem}
\begin{proof}
See \cite[Lemma 2.2, pp. 35--36]{Samko}.
\end{proof}
\begin{definition}\label{w_sol}
Let $u\in L^1([0,T];H^1_0(\Omega))$. We define the weak Riemann--Liouville fractional derivative of the order $\alpha$ of $u$, $^{RL}D^\alpha_{0+}u(t)$, as below
\[
\int_0^T \varphi(t){^{RL}D^\alpha_{0+}}u(t)dt=\int_0^T{^{RL}D^\alpha_{T^-}}\varphi(t)u(t)dt
\]
for all $\varphi\in C^\infty_c((0,T);\mathbb{R})$.
\end{definition}
\section{Weak solutions to time-fractional elliptic equations}
Let $\Omega$ be a bounded domain in $\mathbb{R}^d$ with the boundary $\partial \Omega\in C^1$, $T>0$ and $\alpha\in (0,1)$. Denote $\Omega_T=(0,T]\times\Omega$. We consider the equation of the order $\alpha$
\begin{align}\label{main_eq}
\frac{\partial^{\alpha} u(t,x)}{\partial t^\alpha}-\sum_{i,j=1}^d \partial_{x_i}(a_{ij}(t,x)\partial_{x_j} u(t,x))+\sum_{j=1}^d b_j(t,x)\partial_{x_j}u(t,x)+c(t,x)u(t,x)=f(t)
\end{align}
for $(t,x)\in \Omega_T,$ where $\frac{\partial^{\alpha} u(t,\cdot)}{\partial t^\alpha}$ is the weak Riemann--Liouville fractional derivative of the order $\alpha$ of $u$ with respect to the time variable $t$ and
\begin{itemize}
\item[(a1)] $a_{ij}, b_j, c\in L^\infty (\Omega_T;\mathbb{R})$ for all $1\leq i,j\leq d$;
\item[(a2)] $a_{ij}=a_{ji}$ for all $1\leq i,j\leq d$;
\item[(a3)] there exists $\theta>0$ such that $\sum_{i,j=1}^d a_{ij}(t,x)\xi_i\xi_j\geq \theta \|\xi\|^2$ for a.e. $t\in (0,T)$, $x\in \Omega$ and for all $\xi\in \mathbb{R}^d$;
\item[(a4)] $f\in L^\infty([0,T];H^{-1}(\Omega))$.
\end{itemize}
Assume that
\begin{equation}\label{boundary_eq}
u(t,x)=0\; \text{on}\;[0,T]\time\partial\times\partial \Omega.
\end{equation}
Denote
\[
a(u,v;t):=\int_{\Omega}\sum_{1\leq i,j\leq d}a_{ij}\partial_{x_j}u\partial_{x_i}v+\sum_{1\leq j\leq d}b_j\partial_{x_j}u v+cuv
\]
for all $u,v\in H^1_0(\Omega)$ and $\bold{u}(t):=[u(t)](x)$ for all $t\in [0,T]$, $x\in \Omega$.
\begin{definition}[Weak solution]
A function $\bold{u}:[0,T]\rightarrow H^1_0(\Omega)$ is a weak solution to the problem \eqref{main_eq} with the condition \eqref{boundary_eq} if
\begin{itemize}
\item[(i)] $\bold{u}\in L^2([0,T];H^1_0(\Omega))$ and $^{RL}D^\alpha_{0+}\bold{u}\in L^2([0,T];H^{-1}(\Omega))$;
\item[(ii)] for all $v\in H^1_0(\Omega)$
\[
\langle ^{RL}D^\alpha_{0+}\bold{u}(t),v\rangle_{H^{-1}\times H^1_0}+a(\bold{u}(t),v;t)=\langle f(t),v\rangle_{H^{-1}\times H^1_0}
\]
for a.e. $t\in (0,T)$, where $\langle \cdot,\cdot\rangle_{H^{-1}\times H^1_0}$ is the duality pairing between $H^{-1}(\Omega)$ and $H^1_0(\Omega)$.
\end{itemize}
\end{definition}
\subsection{Galerkin approximation solution}
Let $\{e_j\}_{j=1}^\infty$ be smooth functions which constitutes an orthonormal basic of $L^2(\Omega)$ and a basis of $H^1_0(\Omega)$ such that
\begin{itemize}
\item[(i)] $-\Delta e_k=\lambda_k e_k$, $k\in \mathbb{N}$;
\item[(ii)] \[
\int_\Omega e_j e_j=\begin{cases}
1,\quad i=j\\
0,\quad i\ne j
\end{cases}
\text{and}\;
\int_{\Omega}De_i De_j =\begin{cases}
\lambda_i,\quad i=j\\
0,\quad i\ne j
\end{cases}
\]
\end{itemize}
For the existence of these functions, see \cite[p. 335]{Evans_98}. Fix $N\in \mathbb{N}$. Let $E_N=\text{span}\{e_1,\dots,e_N\}$ and
$P_N$ be the project map from $L^2(\Omega)$ to $E_N$ defined by
$P_Nu=\sum_{i=1}^Nc^ie_i$ with $u\in L^2(\Omega)$ has the form $u=\sum_{i=1}^\infty c^i e_i$ in which $c^i=\int_{\Omega}ue_idx$, $i\in\mathbb{N}$. In this section, using Galerkin method, we will construct approximation solutions to the problem \eqref{main_eq}--\eqref{boundary_eq} in $E_N$.

Let the function $\bold{u}_N:[0,T]\rightarrow E_N$ having the form
\begin{equation}\label{Galerkin_sol}
\bold{u}_N(t)=\sum_{i=1}^N c^i(t)e_i,\quad t\in [0,T],
\end{equation}
here $c^i(\cdot)$, $1\leq i\leq N$, is continuous and has the Riemann--Liouville fractional derivative of the order $\alpha$ on $[0,T]$. Assume that
\begin{equation}\label{gas}
(^{RL}D^\alpha_{0+}\bold{u}_{N},v)_{L^2}+a(\bold{u}_N(t),v;t)=\langle f(t),v\rangle_{H^{-1}\times H^1_0}
\end{equation}
for a.e. $t\in (0,T)$ and
\begin{equation}\label{ini_cond}
\bold{u}_N(0)=0.
\end{equation}
We obtain the following result.
\begin{proposition}
For any $N\in \mathbb{N}$, there exists a unique solution to the problem \eqref{gas}--\eqref{ini_cond} having the form \eqref{Galerkin_sol}.
\end{proposition}
\begin{proof}
Consider the function $\bold{u}_N:[0,T]\rightarrow E_N$ having the form
\[
\bold{u}_N(t)=\sum_{i=1}^N c^i(t)e_i,\quad t\in [0,T],
\]
where $c^i(\cdot)$, $1\leq i\leq N$, is continuous and has the Riemann--Liouville fractional derivative of the order $\alpha$ on $[0,T]$. To $\bold{u}_N(\cdot)$ is a solution to \eqref{gas}--\eqref{ini_cond} then $c^i(\cdot)$, $1\leq i\leq N$, have to satisfy the following condition
\begin{align}
&(\sum_{i=1}^N {^{RL}D^\alpha_{0+}c^{i}(t)}e_i,e_j)_{L^2}+a(\sum_{i=1}^N c^i(t)e_i,e_j;t)\label{e1}\\
\notag&\hspace{1cm}={^{RL}D^\alpha_{0+}}c^{j}(t)+\sum_{i=1}^N c^i(t)a(e_i,e_j;t)=\langle f(t),e_j\rangle_{H^{-1}\times H^1_0}=f^j(t)
\end{align}
for $1\leq j\leq N$ and almost every $t\in (0,T]$. Moreover,
\begin{equation}\label{e2}
c^j(0)=0,\quad 1\leq j\leq N.
\end{equation}
Put $\overrightarrow{c}(t)=(c^1(t),\dots,c^N(t))^{\rm T}$, $A(t)=(a(e_i,e_j;t))_{1\leq i,j\leq N}$, and $$\overrightarrow{f}(t)=(f^1(t),\dots,f^N(t))^{\rm T}.$$ Then the system \eqref{e1}--\eqref{e2} is rewritten in the form
\begin{align}
\label{e3}^{RL}D^\alpha_{0+}\overrightarrow{c}(t)+A(t)\overrightarrow{c}(t)={^{C}D^\alpha_{0+}}\overrightarrow{c}(t)+A(t)\overrightarrow{c}(t)&=\overrightarrow{f}(t),\quad t\in (0,T],\\
\label{e4}\overrightarrow{c}(0)&=0.
\end{align}
Hence, the system \eqref{e3}--\eqref{e4} has a solution on $[0,T]$ if and only if the following integral equation has a continuous solution
\begin{equation}\label{equi}
\overrightarrow{c}(t)=-\frac{1}{\Gamma(\alpha)}\int_0^t (t-\tau)^{\alpha-1}A(\tau)\overrightarrow{c}(\tau)d\tau+\frac{1}{\Gamma(\alpha)}\int_0^t (t-\tau)^{\alpha-1}\overrightarrow{f}(\tau)d\tau,\;t\in [0,T].
\end{equation}
On the space $C([0,T];\mathbb{R}^N)$, we establish an operator as
\[
\mathcal{T}_0\varphi(t)=-\frac{1}{\Gamma(\alpha)}\int_0^t (t-\tau)^{\alpha-1}A(\tau)\varphi(\tau)d\tau+\frac{1}{\Gamma(\alpha)}\int_0^t (t-\tau)^{\alpha-1}\overrightarrow{f}(\tau)d\tau,\;t\in [0,T].
\]
For any $\gamma>0$, define a norm $\|\cdot\|_\gamma$ on $C([0,T];\mathbb{R}^N)$ by
\[
\|\varphi\|_\gamma:=\max_{t\in [0,T]}\frac{\|\varphi(t)\|}{\exp{(\gamma t)}}.
\]
It is obvious that $(C([0,T];\mathbb{R}^N), \|\cdot\|_\gamma)$ is a Banach space. On the other hand, for any $\varphi,\tilde{\varphi}\in C([0,T];\mathbb{R}^N)$ and $t\in [0,T]$, we have
\begin{align*}
\frac{\|\mathcal{T}_0\varphi(t)-\mathcal{T}_0\tilde{\varphi}(t)\|}{\exp{(\gamma t)}}&\leq \frac{\text{ess sup}_{t\in [0,T]}\|A(t)\|}{\Gamma(\alpha)}\int_0^t (t-\tau)^{\alpha-1}\exp{(-\gamma(t-\tau))}\frac{\|\varphi(\tau)-\tilde{\varphi}(\tau)\|}{\exp{(\gamma \tau)}}d\tau\\
&\leq \frac{\text{ess sup}_{t\in [0,T]}\|A(t)\|}{\gamma^\alpha\Gamma(\alpha)}\int_0^{\gamma t} u^{\alpha-1}\exp{(-u)}du\|\varphi-\tilde{\varphi}\|_\gamma\\
&\leq\frac{\text{ess sup}_{t\in [0,T]}\|A(t)\|}{\gamma^\alpha}\|\varphi-\tilde{\varphi}\|_\gamma.
\end{align*}
Hence, for $\gamma>0$ large enough, the operator $\mathcal{T}_0$ is contractive in $(C([0,T];\mathbb{R}^N), \|\cdot\|_\gamma)$ and has a unique fixed point which is also the solution to the system \eqref{equi}. The proof is complete.
\end{proof}
We now give some estimate of the Galerkin approximation solution.
\begin{proposition}
There exists a positive constant $C$, depending on $\Omega$, $T$ and the coefficients of the equation \eqref{main_eq}, such that for all $N\in \mathbb{N}$ it holds that
\begin{align*}
&\|\bold{u}_N\|^2_{L^2([0,T];L^2(\Omega))}+\|\bold{u}_N\|^2_{L^2([0,T];H^1_0(\Omega))}+\|^{RL}D^\alpha_{0+}\bold{u}_N\|^2_{L^2([0,T];H^{-1}(\Omega))}\\
&\hspace{2cm}\leq C\|f\|^2_{L^\infty([0,T];H^{-1}(\Omega)}.
\end{align*}
\end{proposition}
\begin{proof}
First, using \cite[Theorem 2]{Tuan_18}, we have
\begin{align}
\notag \langle^{RL}D^\alpha_{0+}\bold{u}_N(t),\bold{u}_N(t)\rangle_{H^{-1}\times H^1_0}&=\langle^{C}D^\alpha_{0+}\bold{u}_N(t),\bold{u}_N(t)\rangle_{H^{-1}\times H^1_0}\\
\notag&=\Big(\sum_{i=1}^N{^{C}D^\alpha_{0+}}c^i(t)e_i,\sum_{i=1}^N c^i(t)e_i\Big)_{L^2}\\
\notag &=\sum_{i=1}^N {c^i(t)^{C}D^\alpha_{0+}}c^i(t)\\
\notag&\geq \frac{1}{2}\sum_{i=1}^N {^{C}D^\alpha_{0+}}(c^i(t))^2\\
\notag&=\frac{1}{2}{^{C}D^\alpha_{0+}}\sum_{i=1}^N (c^i(t))^2\\
\notag&=\frac{1}{2}{^{C}D^\alpha_{0+}}\|\bold{u}_N(t)\|^2_{L^2(\Omega)}.
\end{align}
On the other hand, by \cite[Theorem 3, p. 300]{Evans_98}, there exist $\beta>0$ and $\nu\geq 0$ such that
\begin{equation*}
\beta \|\bold{u}_N(t)\|^2_{H^1_0}\leq a(\bold{u}_N(t),\bold{u}_N(t);t)+\nu \|\bold{u}_N(t)\|^2_{L^2}\; \text{for almost every}\; t\in (0,T).
\end{equation*}
Moreover, from the assumption of $f$ and the Cauchy inequality
\begin{align}
\notag\langle f(t),\bold{u}_N(t)\rangle_{H^{-1}\times H^1_0}&\leq \|f(t)\|_{H^{-1}(\Omega)}\|\bold{u}_N(t)\|_{H^1_0(\Omega)}\\
\notag&\leq \frac{1}{4\beta}\|f(t)\|^2_{H^{-1}(\Omega)}+\beta \|\bold{u}_N(t)\|^2_{H^1_0(\Omega)}.
\end{align}
Thus, for almost every $t\in (0,T)$,
\begin{equation}
\notag ^{C}D^\alpha_{0+}\|\bold{u}_N(t)\|^2_{L^2(\Omega)}\leq 2\nu \|\bold{u}_N(t)\|^2_{L^2(\Omega)}+\frac{1}{2\beta}\|f(t)\|^2_{H^{-1}(\Omega)}.
\end{equation}
Put $v(t):=\|\bold{u}_N(t)\|^2_{L^2}$, $h(t):=\frac{1}{2\beta}\|f(t)\|^2_{H^{-1}}$ and use the comparison principle for solutions to fractional differential equation and the variation of constants formula for solutions to the equations (see \cite[Lemma 3.1]{Tuan_17}), we obtain the estimate
\begin{equation}
\notag v(t)\leq\int_0^t (t-\tau)^{\alpha-1}E_{\alpha,\alpha}(2\nu (t-\tau)^\alpha)h(\tau)d\tau \;\text{for a.a.}\; t\in [0,T].
\end{equation}
This implies that
\begin{equation}\label{e5}
\|\bold{u}_N(t)\|^2_{L^2(\Omega)}\leq \frac{T^\alpha E_{\alpha,\alpha+1}(2\nu T^\alpha)\|f\|^2_{L^\infty(0,T;H^{-1}(\Omega))}}{2\beta}\;\text{for a.a.}\; t\in [0,T].
\end{equation}
Next, by the similar arguments as above, we see that
\begin{align}
\notag ^{C}D^\alpha_{0+}\|\bold{u}_N(t)\|^2_{L^2(\Omega)}+\beta \|\bold{u}_N(t)\|^2_{H^1_0(\Omega)}&\leq \frac{1}{\beta}\|f(t)\|^2_{H^{-1}(\Omega)}+2\nu \|\bold{u}_N(t)\|^2_{L^2(\Omega)}.
\end{align}
Furthermore,
\begin{align*}
^{C}D^\alpha_{0+}\|\bold{u}_N(t)\|^2_{L^2(\Omega)}={^{RL}D^\alpha_{0+}}\|\bold{u}_N(t)\|^2_{L^2(\Omega)}=\frac{d}{dt}I^{1-\alpha}_{0+}(\|\bold{u}_N(t)\|^2_{L^2(\Omega)}),
\end{align*}
\begin{align}
\notag I^{1-\alpha}_{0+}(\|\bold{u}_N(t)\|^2_{L^2(\Omega)})|_0^T&+\beta\|\bold{u}_N\|^2_{L^2(0,T;H^1_0(\Omega))}\\
\notag&\leq \frac{T\|f\|^2_{L^\infty(0,T;H^{-1}(\Omega))}}{\beta}+2\nu \int_0^T \|\bold{u}_N(t)\|^2_{L^2(\Omega)}dt,
\end{align}
Thus,
\[
\beta\|\bold{u}_N\|^2_{L^2([0,T];H^1_0(\Omega))}\leq \frac{T\|f\|^2_{L^\infty([0,T];H^{-1}(\Omega))}}{\beta}+2\nu \int_0^T \|\bold{u}_N(t)\|^2_{L^2(\Omega)}dt.
\]
This combines with \eqref{e5} implies that there is a constant $C_1>0$ such that
\begin{equation}\label{e6}
\|\bold{u}_N\|^2_{L^2([0,T];H^1_0(\Omega))}\leq C_1 \|f\|^2_{L^\infty ([0,T];H^{-1}(\Omega))}.
\end{equation}
Finally, fix any $v\in H^1_0(\Omega)$ with $\|v\|_{H^1_0}\leq 1$, and write $v=v_0+v_1$, where $v_1\in \text{span}\{e_1,\dots,e_N\}$ and $(e_i,v_0)=0$, $1\leq i\leq N$. Using the estimate concerning the bilinear operator $a(\cdot,\cdot)$ and the facts that $$(^{RL}D^\alpha_{0+}u_N(t),v_1)_{L^2}=\langle ^{RL}D^\alpha_{0+}u_N(t),v\rangle_{H^{-1}\times H^1_0},$$ $(^{RL}D^\alpha_{0+}u_N(t),v_1)_{L^2}+a(u_N(t),v_1;t)=\langle f(t),v_1\rangle_{H^{-1}\times H^1_0}$ and $\|v_1\|_{H^1_0(\Omega)}\leq \|v\|_{H^1_0(\Omega)}\leq 1$, we obtain
\begin{align}
\notag|\langle ^{RL}D^\alpha_{0+}\bold{u}_N(t), v\rangle_{H^{-1}\times H^1_0}|&\leq |a(\bold{u}_N(t),v_1;t)|+|\langle f(t),v_1\rangle_{H^{-1}\times H^1_0}|\\
\notag & \leq C_2\|\bold{u}_N(t)\|_{H^1_0(\Omega)}\|v\|_{H^1_0(\Omega)}+\|f(t)\|_{H^{-1}(\Omega)}\|v\|_{H^1_0(\Omega)}.
\end{align}
Hence,
\begin{equation}\label{e7}
\notag \|^{RL}D^\alpha_{0+}\bold{u}_N(t)\|_{H^{-1}(\Omega)}\leq C_2\|\bold{u}_N(t)\|_{H^1_0(\Omega)}+\|f(t)\|_{H^{-1}(\Omega)}
\end{equation}
which together with \eqref{e7} and the estimate \eqref{e6} implies
\begin{equation}\label{e8}
\|^{RL}D^\alpha_{0+}\bold{u}_N\|^2_{L^2([0,T];H^{-1}(\Omega))}\leq C_3\|f\|^2_{L^\infty([0,T];H^{-1}(\Omega))}.
\end{equation}
From \eqref{e5}, \eqref{e6} and \eqref{e8}, the proof is complete.
\end{proof}
\subsection{Existence and uniqueness of weak solutions}
We now in a position to state the main result of the paper.
\begin{theorem}
Consider the problem \eqref{main_eq}--\eqref{boundary_eq}. Suppose that assumptions (a1)--(a4) hold. Then, this problem has a unique weak solution.
\end{theorem}
\begin{proof}
First, we prove the system \eqref{main_eq}--\eqref{boundary_eq} has at least one weak solution. From \eqref{e6} and \eqref{e8}, by the Banach--Aloaglu theorem, there exist a sequence $\{n_k\}_{k=1}^\infty$ such that
\begin{align}
\bold{u}_{n_k}&\rightharpoonup u\;\text{in}\; L^2([0,T];H^1_0(\Omega)),\label{e9}\\
^{RL}D^\alpha_{0+}\bold{u}_{n_k}&\rightharpoonup v\;\text{in}\;L^2([0,T];H^{-1}(\Omega)).\label{e10}
\end{align}
Let $\varphi\in C^\infty_c((0,T);\mathbb{R})$ and $\psi\in H^1_0(\Omega)$ be arbitrary. Then,
\begin{align*}
\notag\int_0^T \langle v(t),\phi(t)\psi\rangle_{H^{-1}\times H^1_0}dt&=\lim_{k\to\infty}\int_0^T \langle ^{RL}D^\alpha_{0+}\bold{u}_{n_k}(t),\phi(t)\psi\rangle_{H^{-1}\times H^1_0}dt\\
\notag&=\lim_{k\to\infty}\int_0^T \phi(t)\langle ^{RL}D^\alpha_{0+}\bold{u}_{n_k}(t),\psi\rangle_{H^{-1}\times H^1_0}dt\\
\notag&=\lim_{k\to\infty}\int_0^T {^{RL}D^\alpha_{T-}}\phi(t)\langle \bold{u}_{n_k}(t),\psi\rangle_{H^{-1}\times H^1_0}dt\\
\notag&=\int_0^T {^{RL}D^\alpha_{T-}}\phi(t)\langle \bold{u}(t),\psi\rangle_{H^{-1}\times H^1_0}dt\\
\notag&=\int_0^T\langle {^{RL}D^\alpha_{0+}}\bold{u}(t),\varphi(t)\psi\rangle_{H^{-1}\times H^1_0}dt\\
&=\int_0^T\varphi(t)\langle {^{RL}D^\alpha_{0+}}\bold{u}(t),\psi\rangle_{H^{-1}\times H^1_0}dt,
\end{align*}
which implies
\begin{equation}\label{e_11}
v(t)={^{RL}D^\alpha_{0+}}\bold{u}(t).
\end{equation}
Fix $N,M\in \mathbb{N}$ and $N>M$. For any $\varphi\in C^\infty_0([0,T];\mathbb{R})$ and $w\in E_M$, we see that
\begin{align*}
\notag\int_0^T \langle {^{RL}D^\alpha_{0+}\bold{u}_N(t)},\varphi(t)w\rangle_{H^{-1}\times H^1_0}dt&=\int_0^T \varphi(t)\langle {^{RL}D^\alpha_{0+}\bold{u}_N(t)},w\rangle_{H^{-1}\times H^1_0}\\
&\rightarrow \int_0^T \varphi(t)\langle {^{RL}D^\alpha_{0+}\bold{u}(t)},w\rangle_{H^{-1}\times H^1_0},
\end{align*}
\begin{align*}
\notag\int_0^Ta(\bold{u}_N(t),\varphi(t)w;t)dt&=\int_0^T\varphi(t)a(\bold{u}_N(t),w;t)dt\\
&\rightarrow\int_0^T\varphi(t)a(\bold{u}(t),w;t)dt,
\end{align*}
and
\begin{equation*}
\int_0^T \langle f(t),\varphi(t)w\rangle_{H^{-1}\times H^1_0}dt=\int_0^T \varphi(t)\langle f(t),w\rangle_{H^{-1}\times H^1_0}dt.
\end{equation*}
Thus,
\begin{equation}\label{e13}
\langle {^{RL}D^\alpha_{0+}\bold{u}(t)},w\rangle_{H^{-1}\times H^1_0}+a(\bold{u}(t),w;t)=\langle f(t),w\rangle_{H^{-1}\times H^1_0}
\end{equation}
for all $w\in E_M$. Moreover, $\bigcup_{M\in \mathbb{N}}E_M$ is dense in $H^1_0(\Omega)$ which implies that \eqref{e13} is also true for all $w\in H^1_0(\Omega)$ and thus the problem \eqref{main_eq}--\eqref{boundary_eq} has a weak solution.

Next, we give a proof for the uniqueness of the weak solution to \eqref{main_eq}--\eqref{boundary_eq}. Suppose that $u_1, u_2$ are two weak solution to this system. Define $u:=u_1-u_2$. Then, $u$ satisfies
\begin{equation*}
\langle (k\ast u)'(t),v\rangle_{H^{-1}\times H^1_0}+a(u(t),v;t)=0,\;v\in H^1_0(\Omega),\;a.a.\;t\in (0,T).
\end{equation*}
Let $v=u(t)$, then for almost every $t\in (0,T)$
\begin{align}
\label{eq_ex_1}\langle^{RL}D^\alpha_{0+}u(t),u(t)\rangle_{H^{-1}\times H^1_0}+a(u(t),u(t);t)&=\langle (k\ast u)'(t),u(t)\rangle_{H^{-1}\times H^1_0}\\
\notag&\hspace{-1cm}+a(u(t),u(t);t)=0,
\end{align}
where $k(t)=\frac{1}{\Gamma(1-\alpha)t^\alpha}$ for $t>0$. Motivated by Rico Zacher \cite[pp. 291--292]{Zacher_08}, we will approximate the operator $\frac{d}{dt}(k\ast u)$ by the sequences $\{\frac{d}{dt}(k_n\ast u)\}_n$, where $k_n(t):=ns(t)=nE_\alpha(-nt^\alpha)$ (note that by using Laplace transform, we see that $s(\cdot)$ is the unique solution to the equation $s(t)+n(l\ast s)(t)=1$, $t>0$, where $l(t)=\frac{t^{\alpha-1}}{\Gamma(\alpha)}$, $t>0$, see \cite[p. 3]{Ke_2020}). We rewrite the equation \eqref{eq_ex_1} as below
\begin{equation}\label{eq_ex_2}
\langle (k_n\ast u)'(t),u(t)\rangle_{H^{-1}\times H^1_0}+a(u(t),u(t);t)=h_n(t),\;a.a.\;t\in (0,T)
\end{equation}
with
\[
h_n(t):=\langle (k_n\ast u)'(t)-(k\ast u)'(t),u(t)\rangle_{H^{-1}\times H^1_0},\;a.a.\;t\in (0,T).
\]
By virtue of \cite[Lemma 2.1]{Zacher_09},
\[
\frac{1}{2}\frac{d}{dt}(k_n\ast\|u(\cdot)\|^2_{L^2(\Omega)})(t)\leq (\frac{d}{dt}(k_n\ast u)(t),u(t))_{L^2(\Omega)},\;a.a.\;t\in (0,T).
\]
On the other hand, there exists $\nu\geq 0$ such that
\[
a(u(t),u(t);t)\geq -\nu\|u(t)\|^2_{L^2(\Omega)}.
\]
Hence, from \eqref{eq_ex_2}, we have
\[
\frac{d}{dt}(k_n\ast\|u(\cdot)\|^2_{L^2(\Omega)})(t)\leq 2\nu\|u(t)\|^2_{L^2(\Omega)}+2h_n(t),\;a.a.\;t\in (0,T).
\]
This together with the positivity of $l$ implies that
\[
l\ast\frac{d}{dt}(k_n\ast \|u(\cdot)\|^2_{L^2(\Omega)})\leq 2\nu(l\ast \|u(\cdot)\|^2_{L^2(\Omega)})(t)+2l\ast h_n(t),\; a.a.\; t\in (0,T).
\]
We will show that
\begin{equation}\label{key_est}
\lim_{n\to\infty}h_n=0\;\text{in}\; L^1([0,T];\mathbb{R}).
\end{equation}
From the facts that $u\in L^2([0,T];H^{-1}(\Omega))$, $^{RL}D^\alpha_{0+}u=\frac{d}{dt}(k\ast u)\in L^2([0,T];H^{-1}(\Omega))$, \cite[Proposition 2.1, p. 293]{Zacher_08} and  \cite[Example 2.1, p. 294]{Zacher_08}, we have
\[
k\ast \|u(\cdot)\|^2_{H^{-1}(\Omega)}\in W^1_1([0,T];\mathbb{R}).
\]
This implies that
\[
u\in D(B_2),\; \|u(\cdot)\|^2_{H^{-1}(\Omega)}\in D(B_1),
\]
where
\[
B_1(u)=\frac{d}{dt}k\ast u,\;D(B_1)=\{u\in L^1(0,T):k\ast u\in W^1_1([0,T];\mathbb{R})\},
\]
and
\[
B_2(u)=\frac{d}{dt}k\ast u,\;D(B_2)=\{u\in L^2([0,T];H^{-1}(\Omega)):k\ast u\in W^1_2([0,T];H^{-1}(\Omega))\}.
\]
Hence, by using \cite[Estimate (18), p. 292]{Zacher_08}, we obtain
\begin{equation}\label{uniq_e}
\lim_{n\to \infty}\int_0^T \|\frac{d}{dt}[(k-k_n)\ast u](t)\|^2_{H^{-1}(\Omega)}dt=\lim_{n\to \infty}\|\frac{d}{dt}[(k-k_n)\ast u](\cdot)\|^2_{L^2([0,T];H^{-1}(\Omega))}=0.
\end{equation}
By \eqref{uniq_e} and the Holder inequality, the following estimates hold
\begin{align*}
\lim_{n\to \infty}\int_0^T|h_n(t)|dt&=\lim_{n\to \infty}\int_0^T|\langle (k_n\ast u)'(t)-(k\ast u)'(t),u(t)\rangle_{H^{-1}\times H^1_0}|dt\\
&\leq \lim_{n\to \infty}\int_0^T \|\frac{d}{dt}[(k-k_n)\ast u](t)\|_{H^{-1}(\Omega)}\|u(t)\|_{H^1_0(\Omega)}\\
&\leq \lim_{n\to \infty} \|\frac{d}{dt}[(k-k_n)\ast u](\cdot)\|_{L^2([0,T];H^{-1}(\Omega))}\|u\|_{L^2([0,T];H^1_0(\Omega))}\\
&=0,
\end{align*}
which shows that $\lim_{n\to\infty}h_n=0$ in $L^1([0,T])$.
Notice that
\begin{align*}
(l\ast k)(t)&=\frac{1}{\Gamma(\alpha)\Gamma(1-\alpha)}\int_0^t s^{\alpha-1}(t-s)^{-\alpha}ds\\
&=\frac{1}{\Gamma(\alpha)\Gamma(1-\alpha)}\int_0^1 u^{-\alpha}(1-u)^{\alpha-1}du\\
&=\frac{1}{\Gamma(\alpha)\Gamma(1-\alpha)}B(\alpha,1-\alpha)\\
&=1,\;\forall t>0,
\end{align*}
\[
l\ast\frac{d}{dt}(k_n\ast \|u(\cdot)\|^2_{L^2(\Omega)})=\frac{d}{dt}(k_n\ast l\ast\|u(\cdot)\|^2_{L^2(\Omega)})\to \frac{d}{dt}(k\ast l\ast\|u(\cdot)\|^2_{L^2(\Omega)})=\|u(\cdot)\|^2_{L^2(\Omega)}
\]
in $L^1([0,T];\mathbb{R})$ as $n\to \infty$ (see \cite[Estimate (19), p. 292]{Zacher_08}) and $l\ast h_n\to 0$ in $L^1([0,T];\mathbb{R})$ as $n\to \infty$, we obtain
\begin{equation}\label{f_est}
\|u(t)\|^2_{L^2(\Omega)}\leq 2 \nu(l\ast \|u(\cdot)\|^2_{L^2(\Omega)})(t),\;a.a.\; t\in (0,T).
\end{equation}
From \eqref{f_est}, using a Gronwall type inequality as in \cite[Lemma 7.1.1, p. 188]{Henry}, then $\|u(t)\|^2_{L^2(\Omega)}=0$ a.e. in $(0,T)$, that is, $u_1=u_2$. The proof is complete.
\end{proof}
\begin{remark}
The key point in the proof of the uniqueness of weak solution to the problem \eqref{main_eq}--\eqref{boundary_eq} is to show that  $h_n\to 0$ in $L^1([0,T])$ as $n\to\infty$. The approach proposed in \cite{Zacher_09} can not apply directly to this situation because the operator $B=\frac{d}{dt}(k\ast u)$ with domain $D(B)=\{u\in L^2([0,T];H^1_0(\Omega)):\frac
{d}{dt}k\ast u\in L^2([0,T];H^{-1}(\Omega))\}$ is not m-accretive.
\end{remark}
\begin{remark}
In \cite[Theorem 3.1]{Bockstal}, the author has proved the existence and uniqueness of a weak solution to a class of fractional diffusion equations with Caputo derivative. The right hand side of these equations has the same form as in the equation \eqref{main_eq} in our paper. However, to obtain this result, they need additional following assumptions
\begin{itemize}
\item $\|\partial_t f(t)\|_{H^{-1}(\Omega)}\leq Ct^{-\alpha}$;
\item $b_j=0$ ($1\leq j\leq d$), $\partial_t a_{i,j}\in L^\infty(\overline{\Omega_T})$ ($1\leq i,j\leq d$), and $\partial_t c\in L^\infty (\overline{\Omega_T})$,
\end{itemize}
see \cite[lemmas 3.3, 3.4]{Bockstal}. In this paper, by another approach, we have studied the existence and uniqueness of a weak solution to the fractional-order equation \eqref{main_eq} (with Riemann--Liouville fractional derivative) without the assumptions as mentioned above.
\end{remark}
\begin{remark}
An important step in proving the existence of the weak solution is to show that the function $v$ in \eqref{e10} is the fractional derivative $^{RL}D^\alpha_{0+}u$ (in the weak sense) of the solution $u$. In this paper, we do that by using Definition \ref{w_sol}. This definition is inspired by the integration by the parts formula for Riemann--Liouville fractional derivatives (Theorem \ref{integration_by_parts}). Unfortunately, this property is generally not true for Caputo fractional derivatives (except for the case where these two fractional derivatives coincide). Therefore, the approach as in the present paper cannot be applied to study the existence of weak solutions to partial differential equations with Caputo fractional derivatives in time.
\end{remark}
\section*{Acknowledgement}
This research is supported by The International Center for Research and Postgraduate Training in  Mathematics--Institute of Mathematics--Vietnam Academy of Science and Technology under the Grant ICRTM01-2020.09. A part of this paper was completed while the author was a postdoc at the Vietnam Institute for Advanced Study in Mathematics (VIASM). He would like to thank VIASM for support and hospitality. The author also would like to thank Nguyen Anh Tu, Ha Duc Thai, and the anonymous referee for the valuable time, constructive suggestions, and interesting comments that have helped him improve the quality and presentation of the paper.


\subsection{Our Contributions}
\label{subsec:contrib}
\textit{Lower bound on AoI regret:} We show a lower bound of $\Omega(\log T)$ on AoI regret for instances that have at least one \textit{strictly competitive} arm (formally defined in Section \ref{sec:setting}).

\textit{Performance of variants of UCB and Thompson Sampling:} We show that the AoI regret for Correlated-UCB (CUCB) and Correlated-Thompson Sampling (CTS) (proposed in \cite{gupta2019multiarmed}) is $\text{O}(C \log T) + \text{O}(1)$. Here $C$ is the number of sub-optimal \textit{competitive arms} (formally defined in Section \ref{sec:setting}).

\textit{Empirical validation:} Through simulations, we compare the performance of UCB, Thompson Sampling, CUCB, CTS, and their \textit{AoI-Aware} variants proposed in \cite{fatale2020regret}.

\subsection{Related Work}
Multi-Armed Bandits (MABs) are a sequential decision making framework where at every time step $t$, a choice has to be made between $K$ possible \textit{bandit arms} with unknown statistics.
UCB \cite{lai1985asymptotically} and Thompson Sampling \cite{thompson1933likelihood} are two widely studied algorithms for the MAB problem. In this work, we study variants of these policies more suited for our setting.

Recently, variants of the traditional MAB framework that are capable of incorporating additional structure into the decision making problem have been introduced. In addition to observing rewards, the \textit{contexual bandit} framework \cite{zhou2015survey}, learns a mapping between a context vector $\theta$ and the best arm $k^{*}$. Another example is the structured bandit framework \cite{structured}, in which the mean rewards for all arms as a function of the context $\theta$ are known but $\theta$ itself is hidden. The Correlated Multi-Armed Bandit framework of \cite{gupta2019multiarmed} is a variant of the MAB problem that presents the scheduler with arms whose rewards are not independent of each another. That is, sampling arm $k$ can reveal information about the rewards we can expect from another arm $\ell$. In addition to observing the rewards obtained from sampling an arm $k$, scheduling algorithms cognizant of correlation can track additional side information to identify some arms as sub-optimal without sampling them often, thereby reducing the accumulation of regret.

AoI or Age-of-Information is measure of the freshness of data available at the central monitoring station. AoI has been the focus of a variety of work, and \cite{kosta2017age} can be referred to for a comprehensive survey. Previously, a large focus of the work on AoI has been on problems where channel statistics are known \cite{tripathi2017age, jhunjhunwala2018age, kadota2018optimizing, Sombabu:2018:AAS:3241539.3267734}. In our work, we take channel statistics to be unknown and apply Multi-Armed Bandits to the scheduling problem, which is the approach taken by \cite{fatale2020regret} for minimizing AoI regret. Building upon the work of \cite{fatale2020regret}, we drop the assumption of independence between channels and instead take that they follow the Correlated Multi-Armed Bandit model.

\color{black}
\subsection{Lower Bound for Correlated Bandit Instances}
First, we define the class of $\alpha$-consistent policies and then in Theorem \ref{thm:lower_bound} provide a lower bound on the AoI regret achievable by any policy $\rho$ belonging to this class.
\begin{definition}[$\alpha$-consistent policies \cite{lai1985asymptotically}]
	Let $k_s$ denote the index of the channel scheduled in time-slot $s$. The index $k^{*}$ denotes the index of the optimal channel. A scheduling policy is called $\alpha$-consistent, for a constant $\alpha \in (0,1)$, if there exists an instance dependent constant $M$ such that,
	\begin{align}
		\mathbb{E} \Big[ \sum_{s=1}^{t} \mathds{1}\{k_s = k\} \Big] \leq M t^{\alpha}, \, \, \forall \, k \neq k^{*}.
	\end{align}
	\label{def:alpha_consistent}
\end{definition}

\begin{theorem}[Lower bound on AoI regret]
	If a bandit instance $\mathrm{I}$ has at least one competitive arm $k$ with $\tilde{\Delta}_{k, k^{*}} < 0$, then for any $\alpha$-consistent policy $\rho$, we have,
	\begin{align*}
	R_{\rho}(T) &\geq
	\max_{k \in \mathcal{C'}} \frac{\Delta_k}{\mathrm{D}(P_k, P'_k)} \frac{ (1 - \alpha)\log{T} - \log{(4M)} }{\mu^{*}}.
	\end{align*}
	Otherwise, if $\tilde{\Delta}_{k, k^{*}} \geq 0 \, \, \forall \, k \in [K]$,
$
	R_{\rho}(T) \geq  0.
$
	\label{thm:lower_bound}
\end{theorem}
Here $\mathrm{D}(P_k, P'_k)$ is the KL divergence between the reward distribution of arm $k$ and a suitably chosen perturbed reward distribution. The set $\mathcal{C'}$ is the set of strictly competitive arms and $M$ is an instance dependent constant as in Definition \ref{def:alpha_consistent}.

\subsection{Upper Bound for Correlated Bandit Instances}
We now characterize upper bounds on AoI regret for the policies CUCB and CTS proposed in \cite{gupta2019multiarmed}. The key idea behind these policies is to exploit the correlation in the rewards of various arms to obtain a set of competitive arms in each step. The algorithms then play an arm from this set.
For the sake of completeness, we provide details of these policies in the appendix (Algorithms \ref{algo:CUCB} and \ref{algo:CTS} respectively).
\begin{theorem}[Upper bound on AoI regret under CUCB]
For a Correlated Bandit instance, let the number of competitive sub-optimal arms be equal to $C$. Further, let $\mathcal{C}$ denote the set of competitive arms inclusive of the optimal arm $k^{*}$. Let $\mu_{\min} = \min_k \mu_k$,
\begin{align*}
	t_o &= \inf \Big \{ \tau \geq 2: \Delta_{\min}, \tilde{\Delta}_{k, k^{*}} \geq 4\sqrt{\frac{2K \log \tau }{\tau}} \Big \},\\
    U^{(nc)}_{k, \scalebox{.6}{$\mathrm{CUCB}$}} &= Kt_0 + K^{3}\sum_{t = Kt_0}^{T} 2\Big( \frac{t}{K} \Big)^{-2} + \sum_{t = 1}^{T}3t^{-3},\\
    U^{(c)}_{k, \scalebox{.6}{$\mathrm{CUCB}$}} &= 8 \frac{\log (T)}{\Delta_k^{2}} + \Big (1 + \frac{\pi^2}{3}\Big)

	+ \sum_{t=1}^{T} 2 K t \exp \Big( -\frac{t \Delta_{\min}^{2}}{2K} \Big).
\end{align*}
Then, for $T > t_0$,
\begin{align*}
\mathbb{E} [R_{\scalebox{.6}{$\mathrm{CUCB}$}}(T)] &\leq \frac{1 - \mu^{*}}{\mu^* \mu_{\min}} + \Big ( \frac{1}{\mu_{\min}} - \frac{1}{\mu^{*}}\Big ) \times\\
&\, \,  \bigg( \sum_{k^{'} \in [K] \backslash \mathcal{C}} \Delta_{k^{'}} U^{(nc)}_{k, \scalebox{.6}{$\mathrm{CUCB}$}} + \sum_{k \in \mathcal{C} \backslash \{k^{*}\}} \Delta_k U^{(c)}_{k, \scalebox{.6}{$\mathrm{CUCB}$}} \bigg)\\
&= \text{O}(1) + \text{O}(C \log T),
\end{align*}
and for $T \leq t_0$,
$
\mathbb{E} [R_{\scalebox{.6}{$\mathrm{CUCB}$}}(T)] \leq \big( \frac{1}{\mu_{\min}} - \frac{1}{\mu^{*}} \big) T.
$
\label{thm:upper_bound_cucb}
\end{theorem}

\begin{theorem}[Upper bound on AoI regret under CTS]
For a Correlated Bandit instance, let the number of competitive sub-optimal arms be equal to $C$ and let $\mathcal{C}$ denote the set of competitive arms inclusive of the optimal arm $k^{*}$. Further, let $\mu_{\min} = \min_k \mu_k$,
\begin{align*}
	t_b &= \inf \Big \{ \tau \geq \exp{(11 \beta)}: \Delta_{\min}, \tilde{\Delta}_{k, k^{*}} \geq 6\sqrt{\frac{2K \beta \log \tau }{\tau}} \Big \},\\
    U^{(nc)}_{k, \scalebox{.6}{$\mathrm{CTS}$}} &= K t_b + \sum_{t=1}^{T} 3 t^{-3} \\

	& \ \ + K^{2}\sum_{t = K t_b}^{T} \bigg( (2K + 3) \Big(\frac{t}{K}\Big)^{-2} + \Big( \frac{t}{K}\Big)^{1 - 2\beta}\bigg),\\
    U^{(c)}_{k, \scalebox{.6}{$\mathrm{CTS}$}} &= 18 \frac{\log (T \Delta_{k}^{2})}{\Delta_k^{2}} + \exp{(11\beta)} + \frac{9}{\Delta_{k}^{2}}\\
	& \ \ + \sum_{t=1}^{T} 2 K t \exp \Big( -\frac{t \Delta_{\min}^{2}}{2K} \Big).
\end{align*}
Then, for any choice of $\beta > 1$ and for $T > t_b$,
\begin{align*}
\mathbb{E} [R_{\scalebox{.6}{$\mathrm{CTS}$}}(T)] &\leq \frac{1 - \mu^{*}}{\mu^* \mu_{\min}} + \Big ( \frac{1}{\mu_{\min}} - \frac{1}{\mu^{*}}\Big ) \times\\
\, \, & \bigg( \sum_{k^{'} \in [K] \backslash \mathcal{C}} \Delta_{k^{'}} U^{(nc)}_{k, \scalebox{.6}{$\mathrm{CTS}$}} + \sum_{k \in \mathcal{C} \backslash \{k^{*}\}} \Delta_k U^{(c)}_{k, \scalebox{.6}{$\mathrm{CTS}$}} \bigg)\\
&= \text{O}(1) + \text{O}(C \log T),
\end{align*}
and for $T \leq t_b$,
$
\mathbb{E} [R_{\scalebox{.6}{$\mathrm{CTS}$}}(T)] \leq \big( \frac{1}{\mu_{\min}} - \frac{1}{\mu^{*}} \big) T.
$
\label{thm:upper_bound_cts}
\end{theorem}

\subsection{Discussion of Implications}
From Theorems \ref{thm:upper_bound_cucb} and \ref{thm:upper_bound_cts} it is clear that both CUCB and CTS are $\alpha$-consistent and therefore AoI regret under these policies will satisfy Theorem \ref{thm:lower_bound}. Since, the bounds on AoI regret depend on the number of sub-optimal competitive arms $C$, we consider different possibilities for $C$ to understand the regret bounds. If $C > 0$ for a Correlated Bandit instance, and if at least one arm is strictly competitive, then the lower bound and upper bound on AoI regret are both $\text{O}(\log T)$. However, when there are no competitive arms, that is when $C=0$, then there is no meaningful lower bound on the expected AoI regret. The $C = 0$ case agrees with the fact that the set $\mathcal{C}\backslash\{k^{*}\}$ being empty results in a constant $\text{O}(1)$ upper bound on AoI regret. Hence for both these cases of Correlated Bandit instances, the bounds are order-optimal.

The Correlated Bandit model used in this work assumes the knowledge of deterministic reward functions. In practice, it may be challenging to determine these functions, and while conducting such an exercise we may even be able to learn the distribution of the underlying state $X$ itself, thereby doing away with the need for any Bandit algorithm. However, the strength of this model is that once these functions are determined, they can be utilized in other communication systems with a similar configuration but a different and unknown distribution of $X$. Occlusions of different nature can repeatedly disrupt multiple channels at the same time engendering correlation in their performances. Communication systems similar to the one considered in this work might be placed in very different environments while following a standardized configuration. Due to differences in the environment, occlusion events analogous to the abstract entries in $\{x_1, x_2, \ldots, x_n\}$ used in this work, would occur with different relative frequencies at every installation. This variety in environments would make it difficult to scale solutions customized for every location. The distribution agnostic model and algorithms considered in this work would be highly beneficial in such scenarios.

\section{Introduction}
\input{introduction.tex}

\section{Setting}
\input{setting.tex}

\section{Main Results and Discussion}
\input{main_results.tex}

\section{Simulations}
\input{simulation_results.tex}

\section{Proofs}
\input{proofs.tex}

\section{Conclusions and Future Work}
With next generation communication technologies relying on higher carrier frequencies, the problem of line-of-sight occlusion becomes more significant. In systems of the kind considered in this work, occlusions can effect multiple channels simultaneously resulting in their performances being correlated. This correlation can be exploited to identify certain channels as sub-optimal within a few time steps. We have shown theoretical bounds on AoI regret to corroborate this fact. Moreover, through simulations we observed that policies capable of exploiting correlation perform significantly better than those that disregard the presence of correlation.

The Correlated Bandit framework used in this work does not put constraints on the nature of rewards received. It would likely be beneficial to exploit the fact that in our system, rewards are always either $0$ or $1$. Additionally, the theoretical upper bound on AoI regret for CTS was for \textit{Gaussian} priors. An analogous result remains to be shown for CTS with \textit{Beta} priors. These topics will be the focus of future work.

\bibliographystyle{IEEEtran}

\subsection{Proof of Theorem \ref{thm:lower_bound}}
\label{subsec:lower_lemmas}
\begin{lemma}[Lower bound on AoI regret for any policy]
	Let $n_k(t)$ be the number of times arm $k$ was scheduled in the time slots $1$ to $t-1$. Then, AoI regret $R_{\rho}(T)$ under any policy $\rho$ is lower bounded as,
	\begin{align*}
	R_{\rho}(T) \geq \frac{1}{\mu^{*}} \sum_{i \neq k^{*}} \Delta_{i} \mathbb{E}[n_i (T + 1)],
	\end{align*}
	where $\Delta_{i} = \mu^{*} - \mu_{i}$.
	\label{lemma:aoi_lb_general}
\end{lemma}
\begin{proof}
	Let $S_{\rho}(t)$ and $S^{*}(t)$ be indicator random variables denoting a successful update in time slot $t$ by policy $\rho$ and the oracle's policy respectively. Then, from Lemma 1 of \cite{fatale2020regret} we have,
	\begin{align}
		R_{\rho}(T) \geq \frac{1}{\mu^{*}} \sum_{t=1}^{T} \mathbb{E} [S^{*}(t) - S_{\rho}(t)].
		\label{eq:lb_statement}
	\end{align}
	Further, from \cite{fatale2020regret}, we also have,
	\begin{align}
		\mathbb{E} [S^{*}(t) - S_{\rho}(t)] &\geq \sum_{i \neq k^{*}} (\mu^{*} - \mu_k) \mathbb{P}(\mathds{1}\{ k_t = i \} = 1)\\
		&= \sum_{i \neq k^{*}} \Delta_{i} \mathbb{P} (\mathds{1}\{k_t = i\} = 1).
	\end{align}
	Therefore from inequality (\ref{eq:lb_statement}) we have,
	\begin{align}
	R_{\rho}(T) &\geq \frac{1}{\mu^{*}}\sum_{t=1}^{T} \sum_{i \neq k^{*}} \Delta_i \mathbb{P}(\mathds{1}\{ k_t = i \} = 1)\\
	&= \frac{1}{\mu^{*}} \sum_{i \neq k^{*}} \Delta_{i} \mathbb{E} [n_i (T + 1)].
	\end{align}
\end{proof}

\begin{lemma}[Bretagnolle-Huber Inequality \cite{DBLP:Tsybakov}]\label{lemma:Bretagnolle_Huber}
	Consider two probability measures $P$ and $Q$, both absolutely continuous with respect to a given measure. Then for any event $\mathcal{A}$,
	\begin{align*}
	P(\mathcal{A}) + Q(\mathcal{A}^{c}) \geq \frac{1}{2} \exp\big( -\mathrm{KL}(P||Q) \big).
	\end{align*}
\end{lemma}

\begin{lemma}[Divergence Decomposition Lemma \cite{lattimore}]
	Let $\nu = (P_1 , \ldots, P_K)$ be the reward distributions associated with one $K$-armed bandit, and let $\nu' = ( P'_1 , \ldots , P'_K)$ be the reward distributions associated with another $K$-armed bandit. Let $\mathbb{P}_{\nu}^{t} = \mathbb{P}_{\nu \rho}^{t}$ and $\mathbb{P}_{\nu'}^{t} = \mathbb{P}_{\nu' \rho}^{t}$ be the joint distributions corresponding to the schedule of bandit arms chosen under policy $\rho$ and the rewards received. Then,
\begin{align*}
	\mathrm{KL}(\mathbb{P}_{\nu}^{t}||\mathbb{P}_{\nu'}^{t}) = \mathrm{D} (\mathbb{P}_{\nu}^{t}, \mathbb{P}_{\nu'}^{t}) = \sum_{i=1}^{K} \mathbb{E}_{\nu} [n_i (t + 1)] \mathrm{D}(P_{i}, P'_{i}),
\end{align*}
	where $\mathbb{E}_{\nu} [n_i (t + 1)]$ is the expected number of pulls of arm $i$ in $t$ rounds of play for the bandit instance described by $\nu$.
	\label{lemma:divDecomposition}
\end{lemma}

\begin{proof}[Proof of Theorem \ref{thm:lower_bound}]
Consider a Correlated Bandit instance $\mathrm{I}$ having just two arms, the optimal arm with index $1$ and a lone sub-optimal arm with index $2$. If the sub-optimal arm is strictly competitive, that is if $\tilde{\Delta}_{2, 1} < 0$, then from Theorem 3 in \cite{gupta2019multiarmed} we can construct a perturbed bandit instance $\mathrm{I'}$ such that $\mathbb{E}_{X'}[\tilde{Y}_{2}(X)] > \mu_{1}$. Let $\mathbb{P}^{t}_{\mathrm{I}}$ and $\mathbb{P}^{t}_{\mathrm{I'}}$ be the distributions corresponding to rewards and scheduled arms in the first $t$ time steps for instances $\mathrm{I}$ and $\mathrm{I'}$ respectively. By construction, arm $2$ will be the optimal arm for the bandit instance $\mathrm{I'}$. Now, using Definition \ref{def:alpha_consistent}, for $\alpha$-consistent policies $\rho$, there exists a constant $M$ such that,
\begin{align}
\mathbb{E}^{t}_{\mathrm{I}} \Big[ \sum_{\tau = 1}^{t} \mathds{1} \{k_{\tau} = 2\} \Big] &\leq M t^{\alpha} \label{eq:markov1}\\
\mathbb{E}^{t}_{\mathrm{I'}} \Big[ \sum_{\tau = 1}^{t} \mathds{1} \{k_{\tau} = 1\} \Big] &\leq M t^{\alpha}.
\label{eq:markov2}
\end{align}
Defining the event $\mathcal{A} = \{n_2(t + 1) > t/2\}$ and using inequalities (\ref{eq:markov1}) and (\ref{eq:markov2}), the following Markov inequalities hold,
\begin{align}
\mathbb{P}^{t}_{\mathrm{I}} (\mathcal{A}) &\leq \frac{2M}{t^{1 - \alpha}}\\
\mathbb{P}^{t}_{\mathrm{I'}} (\mathcal{A}^{c}) &\leq \frac{2M}{t^{1 - \alpha}}.
\end{align}
Now, using Lemma \ref{lemma:Bretagnolle_Huber} we can write,
\begin{align}
\mathrm{D} (\mathbb{P}_{\mathrm{I}}^{t}, \mathbb{P}_{\mathrm{I'}}^{t}) \geq (1 - \alpha)\log{t} - \log{(4M)}.
\label{eq:div_lower_bound}
\end{align}
Next, using Lemma \ref{lemma:divDecomposition} we can expand $\mathrm{D} (\mathbb{P}_{\mathrm{I}}^{t}, \mathbb{P}_{\mathrm{I'}}^{t})$ as,
\begin{align}
\mathrm{D} (\mathbb{P}_{\mathrm{I}}^{t}, \mathbb{P}_{\mathrm{I'}}^{t}) &=  \mathbb{E}_{\mathrm{I}} [n_2(t+1)] \mathrm{D} (P_2, P'_2).
\label{eq:lb_on_pulls}
\end{align}
Combining inequality (\ref{eq:lb_on_pulls}) and Lemma \ref{lemma:aoi_lb_general}, we get the following lower bound on AoI regret for instance $\mathrm{I}$,
\begin{align}
R_{\rho}(T) &\geq \frac{\Delta_2 \Big( (1 - \alpha)\log{T} - \log{(4M)} \Big)}{\mu^{*}{\mathrm{D} (P_2, P'_2)} }. \label{eq:lower_bound_I}
\end{align}
If a Correlated Bandit instance has more than one strictly competitive sub-optimal arm, then the expected number of sub-optimal pulls, and by extension, the lower bound on AoI regret, can only be higher due to the greater exploration required. Hence, in general, if the collection of strictly competitive arms is given by $\mathcal{C'}$, then,
\begin{align}
R_{\rho}(T) &\geq \max_{k \in \mathcal{C'}} \frac{\Delta_k}{\mathrm{D}(P_k, P'_k)} \frac{ (1 - \alpha)\log{T} - \log{(4M)} }{\mu^{*}}. \label{eq:lower_bound}
\end{align}
If among the sub-optimal arms, there is no strictly competitive arm, then the lower bound on AoI Regret is simply $0$.
\end{proof}

\subsection{Proof of Theorems \ref{thm:upper_bound_cucb} and \ref{thm:upper_bound_cts}}
As in ordinary MAB instances, the reward $Y_k(X)$ of each arm $k$ in Correlated Bandit instances is distributed according to the Bernoulli distribution. Hence, under the technical Assumption 1 of \cite{fatale2020regret} the following Lemma would apply to Correlated Bandit instances.
\begin{lemma}[Lemma 4 in \cite{fatale2020regret}]\label{lemma:gen_upper_bound}
	Given the expected number of pulls of sub-optimal arms, the AoI regret over a run of $T$ rounds is upper bounded by,
	\begin{align*}
	\sum_{t=1}^{T} \mathbb{E} [a(t)] - \frac{T}{\mu^{*}} \leq
	\frac{1 - \mu^{*}}{\mu^* \mu_{\min}} + \Big ( \frac{1}{\mu_{\min}} - \frac{1}{\mu^{*}}\Big ) \mathbb{E} \Big[\sum_{t=1}^{T} \mathds{1}_{k_t \neq k^{*}}\Big].
	\label{UB_general}
	\end{align*}
	\label{lemma:aoi_upper_bound}
\end{lemma}

\begin{lemma}[Expected number of pulls of a non-competitive arm, Theorem 1 in \cite{gupta2019multiarmed}]
	The expected number of pulls of a non-competitive sub-optimal arm under CUCB is upper bounded by,
	\begin{align*}
	\mathbb{E} [n_{k}^{(nc)}(T)]  &\leq Kt_0 + K^{3}\sum_{t = Kt_0}^{T} 2\Big( \frac{t}{K} \Big)^{-2} + \sum_{t = 1}^{T}3t^{-3}\\
	&= U^{(nc)}_{k, \scalebox{.6}{$\mathrm{CUCB}$}} = O(1),

	\intertext{and under CTS with Gaussian priors by,}

	\begin{split}
	\mathbb{E} [n_{k}^{(nc)}(T)] &\leq K t_b + \sum_{t=1}^{T} 3 t^{-3} \, \\
	& \quad \, + K^{2}\sum_{t = K t_b}^{T} \bigg( (2K + 3) \Big(\frac{t}{K}\Big)^{-2} + \Big( \frac{t}{K}\Big)^{1 - 2\beta}\bigg)
	\end{split}\\
	&= U^{(nc)}_{k, \scalebox{.6}{$\mathrm{CTS}$}} = O(1),
	\end{align*}
	where $t_0, t_b > 0$ are constants defined as,
	\begin{align*}
	t_o &= \inf \Big \{ \tau \geq 2: \Delta_{\min}, \tilde{\Delta}_{k, k^{*}} \geq 4\sqrt{\frac{2K \log \tau }{\tau}} \Big \},\\
	t_b &= \inf \Big \{ \tau \geq \exp{(11 \beta)}: \Delta_{\min}, \tilde{\Delta}_{k, k^{*}} \geq 6\sqrt{\frac{2K \beta \log \tau }{\tau}} \Big \},
	\end{align*}
	where $\beta > 1$ is a parameter in CTS (Gaussian priors).
	\label{lemma:non_comp}
\end{lemma}

\begin{lemma}[Expected number of pulls of a competitive arm, Theorem 2 in \cite{gupta2019multiarmed}]
	The expected number of pulls of a competitive sub-optimal arm under CUCB is upper bounded by,
	\begin{align*}
	\begin{split}
	\mathbb{E} [n_{k}^{(c)}(T)] &\leq 8 \frac{\log (T)}{\Delta_k^{2}} + \Big (1 + \frac{\pi^2}{3}\Big)\\
	&\qquad \, \,  + \sum_{t=1}^{T} 2 K t \exp \Big( -\frac{t \Delta_{\min}^{2}}{2K} \Big)
	\end{split}\\
	&= U^{(c)}_{k, \scalebox{.6}{$\mathrm{CUCB}$}} =  O(\log T),

	\intertext{and under CTS with Gaussian priors by,}

	\begin{split}
	\mathbb{E} [n_{k}^{(c)}(T)] &\leq 18 \frac{\log (T \Delta_{k}^{2})}{\Delta_k^{2}} + \exp{(11\beta)} + \frac{9}{\Delta_{k}^{2}}\\
	& \qquad \, \, + \sum_{t=1}^{T} 2 K t \exp \Big( -\frac{t \Delta_{\min}^{2}}{2K} \Big)
	\end{split}\\
	&= U^{(c)}_{k, \scalebox{.6}{$\mathrm{CTS}$}} = O(\log T),
	\end{align*}
	where $\beta > 1$ is a parameter in CTS (Gaussian priors).
	\label{lemma:comp}
\end{lemma}

\begin{proof}[Proof of Theorems \ref{thm:upper_bound_cucb} and \ref{thm:upper_bound_cts}]
	The results follow by substituting the appropriate expression for the total number of sub-optimal pulls from Lemmas \ref{lemma:non_comp} and \ref{lemma:comp} into Lemma \ref{lemma:gen_upper_bound}.
\end{proof}
\subsection{Correlated Bandit Model}
\begin{figure}[t]
	\centerline{\includegraphics[scale=0.46]{figures/final_plots/i-1/all_arms_plot_AoI.pdf}}
	\caption{Bandit instance $\mathrm{I}_1$}
	\label{fig:instance_i1}
\end{figure}
\begin{figure}[t]
	\centerline{\includegraphics[scale=0.46]{figures/final_plots/i-6/all_arms_plot_AoI.pdf}}
	\caption{Bandit instance $\mathrm{I}_2$}
	\label{fig:instance_i2}
\end{figure}
For a system with $K$ communication channels we construct a Correlated Bandit instance with the same number of arms. The random variable $X$ captures the underlying state that determines the rewards for the correlated arms. The reward obtained on playing arm $k$ is denoted by a Bernoulli random variable $Y_k(X)$, where $Y_k$ is a known deterministic function.
We define $\mu_k = \mathbb{E}_X [Y_k (X)]$. The optimal arm is denoted by $k^{*}$, and has mean $\mu^{*}$. The difference between the largest mean and the mean of a sub-optimal arm is called the sub-optimality gap and is given by $\Delta_k = \mu^{*} - \mu_{k}$.

Since the distribution of $X$ is unknown, the means $\mu_1, \mu_2, \ldots, \mu_K$ are also not known. MAB algorithms such as UCB empirically estimate the means $\hat{\mu}_k$ for use in decision making. Under the Correlated Bandit model, rewards for all the arms are functions of the same random variable $X$. Hence, we can infer the possible rewards that an arm $\ell$ could have returned if we had chosen arm $\ell$ instead of arm $k$. To exploit this, \cite{gupta2019multiarmed} introduces the notion of \textit{pseudo rewards}.
\begin{definition}[Expected pseudo reward and pseudo gap]
	Pseudo reward for arm $\ell$ with respect to arm $k$ is given by,
	\begin{align}
		s_{\ell, k} (r) &= \sup_{x: Y_k(x) = r} Y_{\ell} (x).
	\end{align}
	Expected pseudo reward in turn is defined as,
	\begin{align}
	{\phi}_{\ell, k} &= \mathbb{E}_X [\,s_{\ell, k} (Y_{k} (X))\,].
	\end{align}
	The pseudo gap is defined as $\tilde{\Delta}_{\ell, k^{*}} = \mu^{*} - \phi_{\ell, k^{*}}$.
	\label{def:pseudo}
\end{definition}
We say arm $k$ is \textit{competitive} if $\tilde{\Delta}_{\ell, k^{*}} \leq 0$ and \textit{strictly competitive} if the inequality is strict. If $\tilde{\Delta}_{\ell, k^{*}} > 0$, we call arm $k$ non-competitive. We use $C$ to denote the number of competitive arms excluding arm $k^{*}$.

\begin{example}
Two examples of Correlated Bandit instances $\mathrm{I}_1$ and $\mathrm{I}_2$ are shown in Figures \ref{fig:instance_i1} and \ref{fig:instance_i2} respectively. Both instances have $K=4$ arms and have $X$ as a discrete random variable taking values in the abstract alphabet $\{x_1, x_2, x_3, x_4\}$. In the figures, the length of the interval corresponding to $x_i$ is equal to $\mathbb{P}\{X = x_i\}$. Consider $\mathrm{I}_1$, in it, the mean reward for the bandit arms is given by, $\mu_1 = 1 \times 0.2$, $\mu_2 = 1\times0.3 + 1\times0.3$, $\mu_3 = 1\times0.3 + 1\times0.2$ and $\mu_4 = 1\times0.3$. Hence arm 2 is optimal with $\mu^{*} = 0.6$. The expected pseudo rewards can be computed using Definition \ref{def:pseudo} to obtain ${\phi}_{1, 2} = 0.4$, ${\phi}_{3, 2} = 1.0$ and ${\phi}_{4, 2} = 0.4$. Hence, for instance $\mathrm{I}_1$ arm 2 is optimal and arm 3 is the only competitive sub-optimal arm. Therefore $C=1$ for instance $\mathrm{I}_1$. Performing a similar analysis on $\mathrm{I}_2$, we find that arm 4 is optimal and no other arm is competitive. In other words $C=0$ for $\mathrm{I}_2$.
\end{example}

\subsection{The AoI Regret Metric}
Here, the notions of Age-of-Information (AoI) and AoI regret, that were earlier described informally, are made precise.
\begin{definition}[Age-of-Information (AoI)]
At the start of time slot $t$, let $a(t)$ denote the AoI at the central monitoring station and let $u(t)$ denote the time index at which the recent most successful update was received by the monitoring station. Then,
$
	a(t) = t- u(t).
$
Alternatively,
\begin{align*}
a(t) =
\begin{cases}
1 & \text{if the update in $t-1$ succeeds} \\
a(t-1) + 1 & \text{otherwise.}
\end{cases}
\end{align*}
\end{definition}
Under a given scheduling policy $\rho$, let $a_{\rho}(t)$ denote the AoI in time slot $t$. Further, let $a^{*}(t)$ be the AoI under an oracle's policy that uses the optimal arm $k^{*}$ in all time slots. The AoI regret at time $T$ is the cumulative difference in expected AoI for the two policies from time-slots $1$ to $T$.
\begin{definition}[Age-of-Information Regret (AoI Regret)]
AoI regret for a policy $\rho$, over $T$ slots is given by,
\begin{align}
R_{\rho}(T) =
\sum_{t=1}^{T} \mathbb{E}[a_{\rho}(t) - a^{*}(t)] = \sum_{t=1}^{T} \mathbb{E}[a_{\rho}(t)] - \frac{T}{\mu^{*}},
\label{eq:aoi_regret_def}
\end{align}
where \eqref{eq:aoi_regret_def} follows from the expectation of a geometric random variable with parameter $\mu^{*}$.
\end{definition}
\color{black}


\subsection{Graph neural networks}
Graph neural networks (GNNs) attempt to learn embeddings (i.e, vector representations) of nodes and/or graphs,
utilizing the mechanisms of neural networks adapted to the topology of graphs.
The core idea of GNNs is to allow messages to pass between neighbors so that
the
representation of each node can incorporate the information
from its neighborhood and thus to enable the GNNs to indirectly learn the graph
structures. Numerous novel architectures for GNNs have been proposed and
tested,
which differ by the information-passing mechanisms. Among the most recent
architectures are graph convolutions~\cite{kipf2016semi},  attention
mechanisms~\cite{velickovic2018graph}, and those inspired by convolutional
neural networks~\cite{gao2019graph,niepert2016learning,defferrard2016convolutional}. The final
embeddings obtained from GNNs can be utilized
for various graph mining tasks, such as node
classification~\cite{kipf2016semi}, link
prediction~\cite{schlichtkrull2018modeling}, and graph
classification~\cite{zhang2018end,dai2016discriminative}.

\subsection{Graph classification by GNN}
In graph classification, GNNs are tasked with predicting the label of an unseen
graph. While node embeddings can be updated within a graph, the elusive step
here is how to combine them into a vector representation of the entire graph
that can distinguish among different labels. Two of the most common approaches
are global pooling~\cite{duvenaud2015convolutional} and hierarchical
pooling~\cite{zhitao2018hierarchical,yao2019graph,lee2019self}. The simplest
ways for
global pooling are global mean and global max of the final node embeddings. In
contrast, hierarchical pooling iteratively reduces the number of nodes either by merging
similar nodes into supernodes~\cite{zhitao2018hierarchical,yao2019graph} or
selecting most significant nodes~\cite{lee2019self} until reaching a final
supernode whose embedding is used to represent the whole
graph.

\subsection{Transfer learning for graphs}
While most existing methods attempt to train GNNs as an end-to-end
classification system, some studies considered transfer learning in which the GNN
is trained on a large dataset before being applied to the task of interest, often in
a much smaller dataset. \cite{hlugzlplstrategies} succeeded in
improving 3 (out of 4 attempted) existing GNNs by transfer learning from other
tasks. Rather than training a GNN to classify a dataset right away, the authors
pre-trained that GNN on another massive dataset (up to 456K
graphs);
then they added a classifier and trained the whole architecture on the graph
classification task. However, transfer learning for graph remains a major
challenge, as \cite{ching2018opportunities,wang2019data} pointed out that
considerable domain knowledge is required to design the appropriate
pre-training procedure.

\subsection{Metric learning}
Metric learning aims to approximate a real-valued distance between two objects, and the most common objectives are contrastive loss and triplet loss.
Some work has focused on metric learning on
graphs~\cite{ktena2018metric,liu2019community,ling2020hierarchical}. Some of
them~\cite{ktena2018metric,liu2019community}
employ a Siamese network structure, in which a twin network sharing the same
weights is applied on a pair of graphs, and the two output vectors acting as
representation of the two graphs are passed through a distance measure.

In computer vision,~\cite{schroff2015facenet} learns metric on triplet
of images, where two (anchor and positive) share the same label and
one (negative) has a different label. The model aims to
minimize the distance between the anchor and the positive, while
maximizing the distance between the anchor and the negative. This inspired our
interest in learning graph metrics with a triplet loss.

\subsubsection{Problem definition}
We tackle the task of graph classification. Given
$\mathcal{D}=\{(G_1,y_1),(G_2,y_2),...\}$ where $y_i \in \mathcal{Y}$ is the
class label of the graph $G_i \in \mathcal{G}$, the goal of graph
classification is to learn a mapping $f : \mathcal{G} \to \mathcal{Y}$ that
maps graphs to the set of class labels and predicts the class labels of unknown
graphs.

\subsubsection{Outline of our method}
Our method combines the advantages of both GNN and metric learning.
Specifically, to facilitate a better accuracy of the classifier, our method
first maps input graphs into vectors in the
Euclidean space such that their corresponding vectors are well-separated
based on classes. Below, we first briefly introduce GNNs and a learning
scheme on triplet loss.
We then describe the two stages of our method: pre-training a GNN and
training a classifier.

\subsection{Graph neural networks}
Various GNNs have been proposed and proven effective in
approximating such a function $f$. Starting from a graph $G=(V,E)$ with node
features $X_G = \{x_v | v \in V \}$, GNNs obtain final embeddings
of nodes $Z_G = \{z_v | v \in V \}$ and a final embedding $h_G$ of the
graph after $K$ layers of aggregation. Specifically, at each $(k+1)$-th layer,
the embedding $h_v^{(k+1)}$ of each node $v$ incorporates the embeddings of itself and its neighboring nodes $N(v)$ at the $k$-th layer as follows:
\begin{equation*}
	h_v^{(k+1)} = MERGE\Big(h_v^{(k)}, AGGR\Big(\{h_u^{(k)} | u \in
	N(v)\}\Big) \Big)
\end{equation*}
The embedding $h_G$ of the graph $G$ is then obtained by pooling all node
embeddings into a single vector as follows:
\begin{equation*}
	h_G = POOL\Big( \{h_v^{(i)} | v \in V; i=1,..,K\}\Big)
\end{equation*}
GNNs differ by how the incorporating function $MERGE$, the
aggregating function $AGGR$, and the final pooling function $POOL$ are
implemented.

\begin{figure}[t]
	\centering
	\vspace{-2.5mm}
	\includegraphics[width=0.96\linewidth]{FIG/stage1.pdf} \\
	\vspace{-1.5mm}
	\caption{\label{fig:stage1}
		First training stage for a GNN. $\hat{f}(.)$ is trained to differentiate graphs of different classes via the triplet loss.
	}
\end{figure}

\begin{figure}[t]
	\centering
	\vspace{-4mm}
	\includegraphics[width=0.85\linewidth]{FIG/stage2.pdf} \\
	\vspace{-2mm}
	\caption{\label{fig:stage2}
		Second training stage for a classifier. After $\hat{f}(.)$ maps $g_i$
		to $\hat{f}(g_i)$, a classifier is added to map $\hat{f}(g_i)$ to the
		class probability prediction. At this step, either the classifier is
		trained independently (\textsc{2STG}\xspace) or the whole architecture is
		trained together (\textsc{2STG+}\xspace).
	}
\end{figure}

\subsection{Metric learning based on triplet loss}
Triplet loss was first introduced in ~\cite{schroff2015facenet}. The core idea
is to enforce a
margin between classes of samples. This results in embeddings of the same
class mapped to a cluster distant apart from that of other classes.
Specifically, given a mapping $\bar{f}$, we wish for a graph $g_a$ (anchor)
to be closer to another graph
$g_p$ (positive) of the same class than to a graph $g_n$ (negative) of another
class by at least a margin $\alpha$, which is a hyperparameter, i.e.,
$$\|\bar{f}(g_a)-\bar{f}(g_p)\|_2^2 +
\alpha <
\|\bar{f}(g_a)-\bar{f}(g_n)\|_2^2.$$

The triplet loss for the whole dataset becomes:
\begin{equation*}
\sum_{(a,p,n)}\max\left(\|\bar{f}(g_a)-\bar{f}(g_p)\|_2^2 -
\|\bar{f}(g_a)-\bar{f}(g_n)\|_2^2 +
\alpha, 0\right)
\end{equation*}
with the summation over all considered triplets.

Our two-stage method combines the power of both GNNs and the metric
learning method, as described below.
\subsection{First training stage (pre-training a GNN)}
\label{subsec:firststage}
In the first training stage (depicted in Fig.~\ref{fig:stage1}), given a GNN architecture
$\hat{f}(.)$, its weights are shared
among a triplet network $T$, which consists of three identical GNN architectures having the same weights as $\hat{f}(.)$. The parameters of $T$ are trained on each triplet
of graphs $(g_a,g_p,g_n)$ (anchor, positive, negative), in which the anchor and
the positive graphs are of the same class while the negative graph is of
another class. $T$ maps a triplet of graphs to a triplet of real-valued vectors in the Euclidean space: $T(g_a, g_p, g_n) =
(\hat{f}(g_a), \hat{f}(g_p), \hat{f}(g_n))$. Ideally, $\hat{f}(g_a)$ and
$\hat{f}(g_p)$ should be close while $\hat{f}(g_n)$ is far from them
both. The triplet loss for $(g_a,g_p,g_n)$ is defined as:
\begin{equation*}
	\max\left(\|\hat{f}(g_a)-\hat{f}(g_p)\|_2^2 -
	\|\hat{f}(g_a)-\hat{f}(g_n)\|_2^2 +
	\alpha,0\right)
\end{equation*}

\subsection{Second training stage (training a classifier)}
\label{subsec:secondstage}
In the second stage, a classifier is either trained
independently, or added on top of the trained GNN and trained together on the
graph classification task (see Fig.~\ref{fig:stage2}).

In summary, we propose two training methods for GNNs: \textsc{2STG}\xspace and
\textsc{2STG+}\xspace, both consist of two stages.
\begin{itemize}
	\item \textbf{\textsc{2STG}\xspace} (Pre-training Setting): In the first stage, the
	GNN maps each triplet of graphs to a triplet of
	Euclidean-space vectors, and the GNN is trained on triplet loss. In the
	second stage, a classifier is trained independently to classify the graph
	embeddings.
	\item \textbf{\textsc{2STG+}\xspace} (Fine-tuning setting): It has the same structure as \textsc{2STG}\xspace
	except that in the second stage, the classifier is plugged on top of
	the trained GNN, and then the whole architecture is trained together in an end-to-end
	manner.
\end{itemize}

Note that our methods are compatible to any GNN model that maps each
graph to a representation vector. As shown in the next section, when applied to
this method, each GNN model outperformed itself in the original setting.

\subsection{GNN architectures}

In order to demonstrate that our two-stage method helps realize better
accuracy of GNNs, for each of the following GNN model, we compare the accuracy obtained in
the original setting versus that from our method:
\begin{itemize}
	\item \textsc{GraphSage}~\cite{hamiltoninductive}: Global
	pooling is used to pool node embeddings into one graph embedding.
	\item \textsc{GAT}~\cite{velickovic2018graph}: Global pooling is also applied on the node embeddings.
	\item \textsc{DiffPool}~\cite{zhitao2018hierarchical}: A
	hierarchical approach for pooling the node embeddings.
	\item \textsc{EigenGCN}~\cite{yao2019graph}: A different design for hierarchical
	pooling.
	\item \textsc{SAGPool}~\cite{lee2019self}: A hierarchical graph pooling with
	self-attention mechanisms.
\end{itemize}
In previous studies, these models were trained end-to-end, mapping each graph
to a prediction of class probabilities.
To further illustrate the competency of our method, we also compare it with a
transfer-learning method \cite{hlugzlplstrategies}.

\input{tab_graph_statistics}
\subsection{Datasets}

To validate the claims, we apply our methods on 12
datasets. Some statistics of the datasets are summarized in
Table.~\ref{tab:graph_stats}.
\subsubsection{Benchmark datasets}
They include some commonly tested binary-class datasets
~\cite{morris2020tudataset} for graph classification: \textsc{DD}, \textsc{MUTAG},
\textsc{MUTAGENICITY}, \textsc{PTC-FM}, \textsc{PROTEINS}, and \textsc{IMDB-B}.
\subsubsection{New York Taxi datasets}
In addition, we also test our
method on NYC Taxi
datasets.\footnote{https://www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page}. We extracted the taxi ridership data in 2019 from New York City (NYC) Taxi
Commission. The areas in New York are represented as
nodes, and each taxi trip is an edge connecting the source and destination
nodes. All taxi trips in an 1-hour interval form a graph, and each dataset
spans
a month of taxi operations. We augmented the binary label for each graph
as taxi trips in weekdays (Mon-Thu) vs. weekend (Fri-Sun). We considered two
taxi
operators (Yellow and Green) and processed data in January, February and March
of 2019, making 6 datasets in total: \textsc{Jan. G.}, \textsc{Feb. G.}, \textsc{Mar. G.}, \textsc{Jan. Y.}, \textsc{Feb. Y.}, and \textsc{Mar. Y.}.

\subsection{Experimental procedure}
We test the ability of each GNN architecture to classify graphs in the
following three
settings:
\begin{itemize}
	\item Original setting: The GNN with a final classifier outputs the
	estimated class
	probabilities, and the weights are updated by the cross-entropy loss
	with respect to the ground truth. We use the implementation
	provided by the authors. To enhance the capacity of the final classifier,
	we tune it by using up to three fully-connected layers and
	select the model based on validation sets.
	\item \textsc{2STG}\xspace and \textsc{2STG+}\xspace: See Section~\ref{sec:method}.

\end{itemize}

We also compare our two-stage method with the transfer-learning
method in~\cite{hlugzlplstrategies}, which also claims the effectiveness of a
pre-training strategy. Out of the 5 GNN models investigated in our work,
\textsc{GraphSage} and \textsc{GAT} are provided with trained weights by
\cite{hlugzlplstrategies}, and they are compared with \textsc{GraphSage}/\textsc{GAT} trained in
\textsc{2STG+}\xspace.

Each dataset is randomly split into three sets: train (80\%), validation
(10\%) and test (10\%).
Details about hyperparameter search are in Appendix~\ref{appendix:hyperparameter}.
The reported results are mean and standard deviation
of test accuracy of five splits.

We initialize node features as learnable features that are
optimized with GNN parameters during training.
While input features are provided in some datasets, we empirically find that using learnable features often leads to better accuracy.

In \textsc{2STG}\xspace and \textsc{2STG+}\xspace, each graph can be anchor once, while the positive and negative graphs are chosen
randomly. The classifier is a multi-layer perceptron (MLP).

\subsection{Results and discussion}

\subsubsection{Improvement by our methods}
We draw some observations from the results of comparing \textsc{2STG}\xspace and \textsc{2STG+}\xspace with the original
setting (Tables~\ref{tab:triplet_ppi}
and~\ref{tab:triplet_taxi}):
\begin{itemize}
	\item Pre-training using triplet loss (i.e., the first stage of \textsc{2STG}\xspace and \textsc{2STG+}\xspace)
	consistently enhances the graph classification accuracy
	of each GNN model by $0.9-5.4\%$
	points, compared to its original setting.
	\item Fine-tuning the weights of GNNs (i.e, the second stage of \textsc{2STG+}\xspace) further improves the
	accuracy from \textsc{2STG}\xspace in some cases by up to 1.3\% points.
\end{itemize}

\input{tab_triplet_ppi}
\input{tab_triplet_taxi}
\input{tab_triplet_pretrain}
\input{tab_triplet_pretrain2}

The improvement in accuracy by \textsc{2STG}\xspace and \textsc{2STG+}\xspace suggests two
possible explanations:
\begin{itemize}
	\item The end-to-end training methods fail to realize the full potential
	of the GNN models. Even if the final classifier of is
	upgraded from a fully-connected layer to an MLP, the accuracy is not as high
	as in \textsc{2STG}\xspace and \textsc{2STG+}\xspace.
	\item Learning meaningful embeddings in between that are fairly separated
	based on classes (see Figure~\ref{fig:dd}), for example through metric
	learning as in our methods,
	 facilitates a better accuracy of the final classifier.
\end{itemize}

\begin{figure}[t]
	\addtolength{\tabcolsep}{-3pt}
	\begin{center}
		\vspace{-2mm}

				\includegraphics[width=1.04in]{FIG/dd_original.pdf}
				\includegraphics[width=1.08in]{FIG/dd_triplet.pdf}
				\includegraphics[width=1.08in]{FIG/dd_finetune.pdf}

		\caption{\label{fig:dd}
			Visualization of the final embeddings in Original (left),
			\textsc{2STG}\xspace (middle) and \textsc{2STG+}\xspace (right) settings (\textsc{DD}
			dataset).
			Instances of the two classes are separated better in \textsc{2STG}\xspace
			and \textsc{2STG+}\xspace
			than in the
			original setting.
		}
	\end{center}
\end{figure}

\subsubsection{Comparison with a transfer learning method}
We compare the pre-trained models of \textsc{GraphSage}/\textsc{GAT}
in~\cite{hlugzlplstrategies} with \textsc{GraphSage}/\textsc{GAT} of the same architecture (5
layers, 300 dimensional hidden units and global mean pooling) trained in
\textsc{2STG+}\xspace.
The pre-trained models are fine-tuned on the datasets.
The mean and standard deviation of test accuracy of 5 splits are reported in Tables~\ref{tab:triplet_pretrain} and~\ref{tab:triplet_pretrain2}.
Despite being pre-trained on a much smaller dataset, \textsc{2STG+}\xspace achieves better accuracy in $83\%$ of the considered cases: up to $2\%$
points in  the benchmark
datasets and up to $3\%$ points in the NYC Taxi
datasets. This validates our claims in Table~\ref{tab:compare}.

\subsubsection{Running time}
For each hyperparameter setting in a dataset, the original setting takes up to
half an hour to train a model.  Due to having two stages, \textsc{2STG}\xspace and
\textsc{2STG+}\xspace take up to an hour for both stages.
In~\cite{hlugzlplstrategies}, the transfer-learning method was reported to
take up to one day to pre-train on a rich dataset.

\subsection{Further analysis}
We study the final embeddings, generated by three training methods: original, \textsc{2STG}\xspace, \textsc{2STG+}\xspace. In Fig.~\ref{fig:dd}, despite having 2 dimensions, the embeddings of the original setting can be reduced to 1 dimension. We suspect that given the potential capacity, specifically the dimension of the final embeddings, \textsc{2STG}\xspace and \textsc{2STG+}\xspace utilize such expressive power better than the original setting by resulting in higher intrinsic dimension, weaker correlation between dimensions, and larger distance between embeddings~\cite{chen2020measuring}. While measuring the distance did not provide clear evidence, the results of the other two measurements validate our hypothesis.

\subsubsection{Intrinsic dimension}
We measure the \textit{intrinsic dimension}, which we define as the dimension needed to retain 99$\%$ of the variance in the final embeddings. In the final embeddings of dimension $d$ of the validation set, we apply Principal Component Analysis (PCA) and measure the cumulative explained variance $V(i)$ achieved by setting the final dimension as $i$ for $i=1,2,..,d$. Assume that $V(j) \leq 0.99 \leq V(j+1)$, we linearly interpolate between $j$ and $(j+1)$ to obtain the dimension with $99\%$ explained variance. Such dimensions are reported in Table.~\ref{tab:eff_dim}. Given $64$ dimensions as the capacity, \textsc{2STG}\xspace and \textsc{2STG+}\xspace retain a higher intrinsic dimension, implying that they utilize this capacity better than the original setting.
\input{tab_eff_dim}
\input{tab_pair_corr}

\begin{figure*}[htbp!]
	\vspace{-3mm}
	\begin{subfigure}[b]{0.196\textwidth}
		\includegraphics[width=1.35in]{FIG/Features/Pairwise/green1_sage_5_std.pdf}
		\caption{\textsc{GraphSage}}
	\end{subfigure}
	\begin{subfigure}[b]{0.196\textwidth}
		\includegraphics[width=1.35in]{FIG/Features/Pairwise/green1_gat_5_std.pdf}
		\caption{\textsc{GAT}}
	\end{subfigure}
	\begin{subfigure}[b]{0.196\textwidth}
	\includegraphics[width=1.35in]{FIG/Features/Pairwise/green1_diffpool_5_std.pdf}
	\caption{\textsc{DiffPool}}
	\end{subfigure}
	\begin{subfigure}[b]{0.196\textwidth}
	\includegraphics[width=1.35in]{FIG/Features/Pairwise/green1_eigen_5_std.pdf}
	\caption{\textsc{EigenGCN}}
	\end{subfigure}
	\begin{subfigure}[b]{0.196\textwidth}
	\includegraphics[width=1.35in]{FIG/Features/Pairwise/green1_sag_5_std.pdf}
	\caption{\textsc{SAGPool}}
	\end{subfigure}
	\caption{\label{fig:green1}
		Average absolute correlation coefficients over epochs on the Jan. G. dataset. Means and standard deviations of 15 runs are plotted. Embedding dimensions quickly become strongly correlated in the orginal setting. Such correlation is much weaker in \textsc{2STG}\xspace and \textsc{2STG+}\xspace.
	}
\end{figure*}

\subsubsection{Correlation between final embedding dimensions}
We measure the average absolute correlation coefficients between final embedding dimensions. Specifically, given final embeddings of dimension five, the
average of absolute values of ${5 \choose 2}$ pairwise correlation coefficients is
calculated. As seeen in Table.~\ref{tab:pair_corr}, in most of the cases, the original setting leads to stronger
correlation than \textsc{2STG}\xspace and \textsc{2STG+}\xspace, and several obtained dimensions
are too strongly correlated, making them redundant. This implies that
\textsc{2STG}\xspace and \textsc{2STG+}\xspace utilize the provided capacity (i.e., dimension of final embeddings) better than the original setting by giving fewer redundant dimensions.
As seen in Fig.~\ref{fig:green1}, the average absolute correlation coefficients in the original setting quickly become close to 1. In \textsc{2STG}\xspace and \textsc{2STG+}\xspace, such correlation is much weaker.

\subsubsection*{Source code}
All source code used in the paper can be found at: \href{https://github.com/manhtuando97/two-stage-gnn}{https://github.com/manhtuando97/two-stage-gnn.}
\section{ Supplementary information}

\section{Datasets}
\label{appendix:datasets}

We tested our training method using 12 datasets.
Some statistics of the datasets are summarized in
Table.~\ref{tab:graph_stats}

\subsection{Benchmark datasets}
These are the commonly tested binary-class datasets
~\cite{morris2020tudataset} for the
graph classification task: \textsc{DD}, \textsc{MUTAG},
\textsc{MUTAGENICITY}, \textsc{PTC-FM}, \textsc{PROTEINS}, and \textsc{IMDB-B}.

\subsection{New York City Taxi datasets}
We extracted the taxi ridership data in 2019 from New York City (NYC) Taxi
Commission. The areas in New York are represented as
nodes, and each taxi trip is an edge connecting the source and destination
nodes. All taxi trips in an 1-hour interval form a graph, and each dataset
spans
a month of taxi operations. We augmented the binary label for each graph
as taxi trips in weekdays (Mon-Thu) vs. weekend (Fri-Sun). We considered two
taxi
operators (Yellow and Green) and processed data in January, February and March
of 2019, making 6 datasets in total: \textsc{Jan. G.}, \textsc{Feb. G.}, \textsc{Mar. G.}, \textsc{Jan. Y.}, \textsc{Feb. Y.}, and \textsc{Mar. Y.}.

\section{GNN architectures}
\label{appendix:gnn}

In order to demonstrate that our two-stage training method helps realize a
better
performance of GNNs, for each GNN architecture, we compared the accuracy
obtained in
the original setting versus that from our method. The GNN
architectures we considered in this work are:
\begin{itemize}
	\item \textsc{GraphSage}~\cite{hamiltoninductive}: This is often used as a
	strong baseline in
	graph classification. After obtaining node embeddings, global mean/max
	pooling is applied to combine all node embeddings into one graph embedding.
	\item \textsc{GAT}~\cite{velickovic2018graph}: Instead of uniformly passing
	neighbor information into
	a node embedding, \cite{velickovic2018graph} employs an attention
	mechanism for the
	importance of each neighboring node.
	\item \textsc{DiffPool}~\cite{zhitao2018hierarchical}: While using the same
	aggregation mechanism as \cite{hamiltoninductive},
	\cite{zhitao2018hierarchical} proposes a
	hierarchical approach to pool the node embeddings.
	Rather than a ``flat-pooling" step at the end, diff-pool repeatedly merges
	nodes into ``supernodes" until there is only one supernode whose embedding is
	treated as the graph embedding.
	\item \textsc{EigenGCN}~\cite{yao2019graph}: Attempting to implement hierarchical
	pooling like \cite{zhitao2018hierarchical}, \cite{yao2019graph}
	formulates a different way to combine nodes and their respective embeddings
	making use of the eigenvectors of the Laplacian matrix.
	\item \textsc{SAGPool}~\cite{lee2019self}: Hierarchical graph pooling employing
	self-attention mechanisms.
\end{itemize}

\section{Hyperparameter search}
\label{appendix:hyperparameter}
For each GNN, the hyperparameters regarding the network
architecture were tuned in the
same search space for the three settings: original, \textsc{2STG}\xspace,
\textsc{2STG+}\xspace. The search space for the dimensions of the input vector, hidden
vector and output vector for all GNNs was $\{16, 32, 64, 96, 128\}$. For
\textsc{DiffPool}, we used three
layers of graph convolution and one \textsc{DiffPool} layer as described in the original
paper. For \textsc{EigenGCN}, we used three pooling operators as it was shown to
achieve
the best performance in the original paper. For \textsc{SAGPool}, we used  three
pooling
layers as explained in the original paper. Other hyperparameters that are
exclusive to each GNN architecture were set to the
default values provided in each paper's original code of each architecture's
authors.

In all three settings, the architecture of the final classifier was also tuned.
The number of fully-connected layers was up to 3 while the search space for
the hidden dimension was $\{2^h | 1 \leq h \leq \log_2(d)\}$, where $d$ is the dimension of the output vector.

The two settings \textsc{2STG}\xspace and \textsc{2STG+}\xspace require an additional
hyperparameter $\alpha$ as the margin in the triplet loss.
While~\cite{schroff2015facenet} found $\alpha=0.2$ to
be effective, we empirically found that this value was too small to separate
instances of different classes. Instead, the search space for $\alpha$ we used
was $\{0.5, 1.0, 1.5, 2.0, 2.5\}$.

\section{Computing specifications \& dependencies}
All experiments were performed on Ubuntu 18.04 LTS running on a machine with four RTX2080Ti GPUs, each of which has $11$GB memory.

The following libraries were used to run the code:
\begin{itemize}
	\item torch 1.7.0,
	\item networkx 2.5,
	\item sklearn 0.23.2,
	\item numpy 1.16.4,
	\item torch-geometric 1.6.3.
\end{itemize}

\section{Introduction}
\label{sec:intro}
\input{010introduction}

\section{Related Works}
\label{sec:background}
\input{020background}

\section{Proposed Method}
\label{sec:method}
\input{030method}

\section{Experiments}
\label{sec:experiments}
\input{040experiments}

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

\bibliographystyle{named}