Let $L$ be the longest sequence that there exists $\epsilon'>\epsilon$ and $\mu_1\cdots \mu_L\in\cD$, for all $t\in[L]$, there exists $\phi\in\Phi$,  $\abs{\EE_{ \mu_t}[\phi(x)]}\geq\epsilon$ and $\sum_{i=1}^{t-1} \bracket{\EE_{ \mu_i}[\phi(x)]}^2\leq \epsilon^2$. We denote $L$ as the bellman eluder dimension $\operatorname{d_{DE}}(\Phi,\cD,\epsilon)$.

Given the function set have low eluder dimension, we have the standard elliptical potential lemma as below:

Given a function class $\phi\in\Phi$ in domain $\mathcal{X}$ with $\abs{\phi(x)}\leq1$. Let $\cD$ be a families of distributions on $\mathcal{X}$. Suppose $\sets{\phi_k}_{[K]}\subset \Phi$ and $\sets{\mu_k}_{[K]}\subset \cD$ be two sequences. If for any $k\in [K]$, $\sum_{i=1}^k\bracket{\EE_{\mu_i}[\phi_k]}^2\leq\beta $, then for any $k\in[K]$, $\sum_{i=1}^k \abs{\EE_{\mu_i}[\phi_i]}\leq\mathcal{\Tilde{O}}\bracket{\sqrt{\operatorname{d_{DE}}(\Phi,\cD,1/K)\beta K}}$.

In this section, we define our bellman eluder dimension as the distributional eluder dimension for the specific function class below:

Given a policy class $\Pi^\dag$, and a PDF function class $\cF$. For $h\in[H]$, we define the function class $\Phi_h$ as $\sets{\TV\bracket{f_h||\cT_{h,\mpi}f_{h+1}}}$, and the distribution family $\cD_h$ as $\mu_h^\mpi$ with domain $\gX_h$ as $s_h^\dag,a_h$. We define the bellman eluder dimension of our problem as:

\begin{align*}
\operatorname{d_{BE}}=\max_{h\in [H]}\operatorname{d_{DE}}(\Phi_h,\cD_h,1/K)\,.
\end{align*}

§.§ Setting

§.§.§ Notation

In this setting, we assume the density function of $Z_h^\mpi\in\mathcal{Z}_h$ belongs to a function class $\cF_h$. We denote the density function of $Z_h^\mpi(s^\dag,a)\in\mathcal{Z}_h$ at point $z\in\RR$ as $f_h^\mpi(z|s^\dag,a)\in\cF_h$.

Consider an upper and lower $\epsilon$-bracketing of $\cF$ under $\norm{\cdot}_1$, denoted as ${\cF}^{\uparrow}$ and ${\cF}^{\downarrow}$. We denote the corresponding lower bracket of $g$ as $g^{\downarrow}$, and the upper bracket of $f$ as $f^{\uparrow}$. Since $g^{\downarrow}$ may not be a valid distribution, we denote the normalized version as: $\underline{g}$, where $\underline{g}=g^{\downarrow}/\int_z g^{\downarrow}(z)$, and $1-\epsilon\leq\int_z g^{\downarrow}(z)\leq 1$. Thus, we have: $\underline{g}\leq(1+2\epsilon)g^{\downarrow}$.

§.§.§ Bellman Completeness

Define the augmented bellman operator:

\begin{equation*}
\mathcal{T}_{h,\mpi}^\dag f_{h+1}(z|s_h^\dag,a_h)=\int \TT({s_{h+1}^\dag|s_h^\dag,a_h})\pi_{h+1}(a_{h+1}|s_{h+1}^\dag) f_{h+1}(z-r_h|s_h^\dag,a_h)\,.
\end{equation*}

Here for completeness we  restate the distributional bellman completeness assumption, and give a corresponding example.
[augmented distributional bellman completeness]
For the density function class $\bm{\cF}=\cF_1\times\cdots\cF_H$ corresponding to the class of random variables $\bm{\gZ}=\mathcal{Z}_1\times\mathcal{Z}_2\cdots\mathcal{Z}_H$, we have for any $h\in[H]$, such that for any $f_{h+1}\in\cF_{h+1}$, we have for any $\pi\in\Pi^\dag$, $\mathcal{T}_{h,\mpi}^\dag f_{h+1}\in\cF_h$.

Instances: for linear MDP with $\PP_h(s_{h+1}|s_h,a_h)=\innerproduct{\phi(s_h,a_h)}{\mu_h(s_{h+1})}$ and suppose the reward is dicretized into a uniform grid of $M$ points $\sets{z_i}_{i=1}^M$. Then we can write the reward distribution as
$\RR(z_i|s_h,a_h)=\innerproduct{\mathbf{1_M}}{\theta(z_i)}$, where $\mathbf{1_M}$ is a $M$ dimensional vector with all the entries being $1$, and $\theta(z_i)$ is a $M$ dimensional vector with all $0$ but the $i$ th entry equal to $\RR(z_i|s_h,a_h)$. Then we have:

\begin{align*}
&\mathcal{T}_{h,\mpi}^\dag f_{h+1}(z|s_h^\dag,a_h)\\
=&\sum_{s_{h+1}}\PP(s_{h+1}|s_h,a_h)\sum_{z_i} \RR(z_i|s_h,a_h)\sum_{s}\pi_{h+1}\bracket{a_{h+1}|(s_{h+1},y_h+z_i)} f_{h+1}\bracket{z-z_i|(s_{h+1},y_h+z_i),a_{h+1}}\\
=& \sum_{s_{h+1}} \sum_{z_i}\phi(s_h,a_h)^\top \mu_h(s_{h+1})\theta_h(z_i)^\top \mathbf{1_M}\sum_{s}\pi_{h+1}\bracket{a_{h+1}|(s_{h+1},y_h+z_i)} f_{h+1}\bracket{z-z_i|(s_{h+1},y_h+z_i),a_{h+1}}\\
=&\phi(s_h,a_h)^\top W_h^\mpi(z,y_h)\mathbf{1_M}\,,
\end{align*}

\begin{align*}
W_h^\mpi(z,y_h)=\sum_{s_{h+1}} \sum_{z_i} \mu_h(s_{h+1})\theta_h(z_i)^\top \sum_{a}\pi_{h+1}\bracket{a_{h+1}|(s_{h+1},y_h+z_i)} f_{h+1}\bracket{z-z_i|(s_{h+1},y_h+z_i),a_{h+1}}
\end{align*}

depends only on $y_h$ and $z$.
We can also write the distribution function in a linear form, with $f_h^\mpi(z|s_h^\dag,a_h)=\Bar{\phi}(s_h,a_h)^\top w_h^\mpi(z,y_h)$, with $\Bar{\phi}(s_h,a_h)=\phi(s_h,a_h)\otimes\mathbf{1_M}\in\RR^{d\times M}$, and $w_h^\mpi(z,y_h)$ is the flattened version of $W_h^\mpi(z,y_h)$ with $w_h^\mpi(z,y_h)[d\times i+m]=W_h^\mpi(z,y_h)[i,m]$. Thus the function class has a linear structure similar to the case of risk neutral setting in linear MDPs [27].

§.§ Estimation of confidence set

Here, we estimate the confidence set $\mathcal{Z}_{k,\mpi}$ via MLE using the density functions $f$.
$(\gH_{k}, \gZ, \mpi,\gamma^{\MLE})$

Input: History information $\gH_k$, density function class $\bm\gF$ of random variable class $\bm\gZ$, and policy $\mpi$.
For all $i\in[k]$, and $(h,\bm{f},\mpi)\in[H]\times\cF\times\Pi^\dag$, sample
$z_{i,h+1}^{\underline{{f}},\underline{\mpi}}\sim \underline{f}_{h+1}\bracket{s_{i,h+1}^\dag,\underline{\pi}_{i,h+1}(s_{i,h+1}^\dag)}$, and let $z_{i,h}^{\underline{{f}},\underline{\mpi}}=z_{i,h+1}^{\underline{{f}},\underline{\mpi}}+r_{i,h}$. We define the version space as:

\begin{align*}
\widehat{\cF}_{k+1,{\mpi}}^{\MLE}=\sets{\bm{f}\in\cF:
\sum_{i=1}^k \log f_h(z_{i,h}^{\underline{{f}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})\geq \max_{f'\in\cF}\sum_{i=1}^k \log f'_h(z_{i,h}^{\underline{{f}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})-\gamma^{\MLE},\forall h\in[H]}\,.
\end{align*}

Return $\widehat{\gZ}_{k+1,\mpi}$ as the set of random variables corresponding to $\widehat{\cF}_{k+1,\mpi}^{\MLE}$.

We now describe the sampling procedure for our target function $g_{h+1}\in\cF$ and $\mpi\in\Pi$. Define $\underline{g}_{h+1}$ and $\underline{\mpi}$ as the normalized lower bracket in Section <ref>.
For $1\leq i\leq k$,
we sample $z_{i,h+1}^{\underline{{g}},\underline{\mpi}}\sim \underline{g}_{h+1}(s_{i,h+1}^\dag,\underline{\pi}_{h+1}(s_{i,h+1}^\dag))$, then we construct a one-step-back sample as $z_{i,h}^{\underline{{g}},\underline{\mpi}}=z_{i,h+1}^{\underline{{g}},\underline{\mpi}}+r_{i,h}$ where $r_{i,h}=y_{i,h+1}-y_{i,h}$. We estimate the likelihood of $f_h$ for $\mathcal{T}_{h,\mpi}^\dag g_{h+1}$ as $\log f_h(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h}) $. Then we can define the MLE confidence set as:

\begin{align*}
\cF_{h,k,g,{\mpi}}^{\MLE}=\sets{f_h\in\cF_h:
\sum_{i=1}^k \log f_h(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})\geq \max_{f'\in\cF}\sum_{i=1}^k \log f'_h(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})-\gamma^{\MLE}}\,,
\end{align*}
where $\gamma^{\MLE}=\log\bracket{\cN_{[\cdot]}\bracket{\cF,\epsilon,\norm{}_{\infty}}\cN_{[\cdot]}(\Pi^\dag,\epsilon,\norm{}_{1})/\delta}+K\epsilon$.
Then, we can show that w.h.p, we have

\begin{align*}
\mathcal{T}_{h,\mpi}^\dag g_{h+1}\in \cF_{h,k,g,\mpi}^{\MLE}
\end{align*}

We can define our version space as:

\begin{align*}
\widehat{\cF}^{\MLE}_{k,\mpi}=\sets{f\in\cF: f_h\in\cF_{h,k,f,\mpi}^{\MLE}, h\in[H]}\,.
\end{align*}

Thus, we can prove that $f^\mpi\in\widehat{\cF}^{\MLE}_{k,\mpi}$, and we have that $Z^\mpi\in\widehat{\gZ}_{k,\mpi}$

§.§.§ Proof of Condition <ref>

In this section we prove the concentration Condition <ref> (Condition <ref> in the Appendix)

Following standard MLE concentration analysis, we have:

For any $f_h\in\cF_h\,,~h\in[H]$, there exists a constant $c$ such that:

\begin{align*}
\sum_{i=1}^k \log f_h(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})\leq \sum_{i=1}^k \log (\cT_{h,{\mpi}}^{\dag}g_{h+1})(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})+c\gamma^{\MLE}\,.
\end{align*}

Consider an upper and lower $\epsilon$-bracketing of $\cF$ under $\norm{\cdot}_1$, denoted as ${\cF}^{\uparrow}$ and ${\cF}^{\downarrow}$. We denote the corresponding lower bracket of $g_{h+1}$ as $g^{\downarrow}$, and the upper bracket of $f$ as $f^{\uparrow}$. Since $g^{\downarrow}$ may not be a valid distribution, we denote the normalized version as: $\underline{g}$, where $\underline{g}=g^{\downarrow}/\int_z g^{\downarrow}(z)$, and $1-\epsilon\leq\int_z g^{\downarrow}(z)\leq 1$

Then, we have:

\begin{align*}
=&\int_{s_h^\dag,a_h,s_{h+1}^\dag,z} \mu_{h}^{\mpi^i}(s_{h}^\dag,a_{h},s_{h+1}^\dag) \frac{\int_{\cA}\underline{\pi}_{h+1}(a_{h+1}|s_{h+1}^\dag) \underline{g}_{h+1}\bracket{z-r_h|s_{h+1}^\dag,a_{h+1}}f^{\uparrow}_h(z|s_h^\dag,a_h)}{\int \TT(s_{h+1}^\dag|s_h^\dag,a_h)\int_{\cA}\pi_{h+1}(a_{h+1}|s_{h+1}^\dag)g^{\downarrow}_{h+1}\bracket{z-r_h|s_{h+1}^\dag,a_{h+1}}}\\
=& \int_{s_h^\dag,a_h}\mu_{h}^{\mpi^i}(s_h^\dag,a_h)\int_z f^{\uparrow}_h(z|s_h^\dag,a_h)\int_{s_{h+1}^\dag}\frac{\TT(s_{h+1}^\dag|s_h^\dag,a_h)\int_{\cA}\underline{\pi}_{h+1}(a_{h+1}|s_{h+1}^\dag) \underline{g}_{h+1}\bracket{z-r_h|s_{h+1}^\dag,a_{h+1}}}{\int_{s_{h+1}^\dag}\TT(s_{h+1}^\dag|s_h^\dag,a_h)\int_{\cA}\pi_{h+1}(a_{h+1}|s_{h+1}^\dag)g^{\downarrow}_{h+1}\bracket{z-r_h|s_{h+1}^\dag,a_{h+1}}}\\
\leq& \int_{s_h^\dag,a_h}\mu_{h}^{\mpi^i}(s_h^\dag,a_h)\int_z f^{\uparrow}_h(z|s_h^\dag,a_h)\int_{s_{h+1}^\dag}\frac{\TT(s_{h+1}^\dag|s_h^\dag,a_h)\int_{\cA}(1+2\epsilon){\pi}^{\downarrow}_{h+1}(a_{h+1}|s_{h+1}^\dag) \underline{g}_{h+1}\bracket{z-r_h|s_{h+1}^\dag,a_{h+1}}}{\int_{s_{h+1}^\dag}\TT(s_{h+1}^\dag|s_h^\dag,a_h)\int_{\cA}\pi_{h+1}(a_{h+1}|s_{h+1}^\dag)g^{\downarrow}_{h+1}\bracket{z-r_h|s_{h+1}^\dag,a_{h+1}}}\\
\leq&\int_{s_h^\dag,a_h}\mu_{h}^{\mpi^i}(s_h^\dag,a_h)\int_z f^{\uparrow}_h(z|s_h^\dag,a_h) (1+2\epsilon)^2\\
\leq&1+\frac{6}{K}\,.
\end{align*}

Thus we obtain the result via Markov inequality:

\begin{align*}
& \PP\bracket{\sum_{i=1}^k f^{\uparrow}_h(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})-\log (\cT_{h,\mpi}^{\dag}g^{\downarrow}_{h+1})(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})\geq \log\bracket{1/\delta}}\\
\leq & \EE\mbracket{\exp\bracket{\sum_{i=1}^k \log\frac{f^{\uparrow}(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})}{\cT_{h,\mpi}^{\dag}g^{\downarrow}(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})}} }\exp(-\log (1/\delta))\\
\leq& e^6\delta\,.
\end{align*}

Applying a union bound, for all $f^{\uparrow}$ and $g^{\downarrow}$, we have w.p. $1-\delta$, there exists a constant $c$,

\begin{align*}
\sum_{i=1}^k \log f^{\uparrow}_h(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})-\log (\cT_{h,{\mpi}}^{\dag}g^{\downarrow}_{h+1})(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})\leq c\beta\,.
\end{align*}

We conclude our result by the definition of upper and lower brackets:

\begin{align*}
&\sum_{i=1}^k \log f_h(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})-\log (\cT_{h,\mpi}^{\dag}g_{h+1})(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})\\
\leq&  \sum_{i=1}^k \log f^{\uparrow}_h(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})-\log (\cT_{h,\mpi}^{\dag}g^{\downarrow}_{h+1})(z_{i,h}^{\underline{{g}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})\\
\leq& c\beta\,.
\end{align*}

As a result, we have:

\begin{align*}
\cT_{h,\mpi}^\dag f_{h+1}\in\cF_{h,k,f,\mpi}^{\MLE}
\end{align*}

by the definition of $\cF_{h,k,f,\mpi}^{\MLE}$

We prove Condition <ref> via induction. If $f_{h+1}^\mpi\in \cF_{h+1,k,f^\mpi,\mpi}^{\MLE}$, From Lemma <ref>, $f_h^\mpi=\cT_{h,\mpi}^\dag f_{h+1}^\mpi\in \cF_{h,k,f^\mpi,\mpi}^{\MLE}$. Since $f_{H}^\mpi\in \cF_{H,k,f^\mpi,\mpi}^{\MLE}$ we have for every $h$, $f_h^\mpi\in\cF_{h,k,f^\mpi,\mpi} $. Thus $\bm{f}^\mpi\in\widehat{\cF}^{\MLE}_{k,\mpi}$. As a result, $Z^\mpi\in\widehat{\gZ}_{k,\mpi}$ for any $k\in[K]$ and $\mpi\in\Pi^\dag$.

§.§.§ Proof of Condition <ref>

In this section we prove the elliptical potential Condition <ref> (Condition <ref> in the Appendix).

The following lemma is the standard result for MLE generalization bound.

We can bound the square TV distance of the bellman error for any $\bm{f}\in\cF_{k,\mpi}$

\begin{align*}
\sum_{i=1}^k\EE_{s_h^\dag,a_h\sim\mu^{\mpi^i}_h}\mbracket{\TV^2\bracket{f_h(s_h^\dag,a_h)||\mathcal{T}_{h,\mpi}^\dag f_{h+1}(s_h^\dag,a_h)}}\leq \mathcal{O}(\gamma^{\MLE})\,.
\end{align*}

Since we have for any $\bm{f}\in\cF_{k,\mpi}$

\begin{align*}
&  \sum_{i=1}^k \log \bracket{\cT_{h,\underline{\mpi}}^\dag \underline{f}_{h+1}}(z_{i,h+1}^{\underline{{f}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})-\log f_h(z_{i,h+1}^{\underline{{f}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})\\
\leq&\sum_{i=1}^k \log \bracket{\cT_{h,{\mpi}^{\downarrow}}^\dag {f}^\downarrow_{h+1}}(z_{i,h+1}^{\underline{{f}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})-2\log(1-\epsilon)-\log f_h(z_{i,h+1}^{\underline{{f}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})
\end{align*}

which holds by the normalization constant of $\underline{f}$ and $\underline{\pi}$. By the definition of the lower bracket function $f^\downarrow_{h+1}$, we have

\begin{align*}
&\sum_{i=1}^k \log \bracket{\cT_{h,{\mpi}^{\downarrow}}^\dag {f}^\downarrow_{h+1}}(z_{i,h+1}^{\underline{{f}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})-2\log(1-\epsilon)-\log f_h(z_{i,h+1}^{\underline{{f}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h}) \\
\leq& \sum_{i=1}^k \log \bracket{\cT_{h,{\mpi}}^\dag {f}_{h+1}}(z_{i,h+1}^{\underline{{f}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})-2\log(1-\epsilon)-\log f_h(z_{i,h+1}^{\underline{{f}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})\\
\leq& \max_{f'\in\cF}\sum_{i=1}^k \log f'_h(z_{i,h}^{\underline{{f}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})-\log f_h(z_{i,h+1}^{\underline{{f}},\underline{\mpi}}|s_{i,h}^\dag,a_{i,h})+\mathcal{O}(k\epsilon)\\
\leq& \gamma^{\MLE}+\mathcal{O}(k\epsilon)\\
\end{align*}

where the second inequality is due to the distributional bellman completeness, which ensures that $\cT_{h,{\mpi}}^\dag {f}_{h+1}\in\cF_h$. The last two inequalities are due to the construction of the confidence set and the choice of $\epsilon=1/K$.

Since the conditional distribution of $z_{i,h+1}^{\underline{{f}},\underline{\mpi}}$ given $s_{i,h}^\dag,a_{i,h}$ is the same as $\mathcal{T}_{h,\underline{\mpi}}^\dag \underline{f}_{h+1}(s_h^\dag,a_h)$, from Lemma <ref>, we have:

\begin{align*}
\sum_{i=1}^k\EE_{s_h^\dag,a_h\sim\mu^{\mpi^i}_h}\mbracket{\TV^2\bracket{f_h(s_h^\dag,a_h)||\mathcal{T}_{h,\underline{\mpi}}^\dag \underline{f}_{h+1}(s_h^\dag,a_h)}}\leq\mathcal{O}(\gamma^{\MLE})\,.
\end{align*}

We also have:

\begin{align*}
&\TV\bracket{\mathcal{T}_{h,{\mpi}}^\dag f_{h+1}(s_h^\dag,a_h)||\mathcal{T}_{h,\underline{\mpi}}^\dag \underline{f}_{h+1}(s_h^\dag,a_h)}\\
\leq &\TV\bracket{\mathcal{T}_{h,{\mpi}}^\dag f_{h+1}(s_h^\dag,a_h)||\cT_{h,\underline{\mpi}}^\dag f_{h+1}(s_h^\dag,a_h)}+\TV\bracket{\mathcal{T}_{h,\underline{\mpi}}^\dag \underline{f}_{h+1}(s_h^\dag,a_h)||\cT_{h,\underline{\mpi}}^\dag f_{h+1}(s_h^\dag,a_h)}\,.
\end{align*}

We can bound the first term as:

\begin{align*}
&\TV\bracket{\mathcal{T}_{h,{\mpi}}^\dag f_{h+1}(s_h^\dag,a_h)||\cT_{h,\underline{\mpi}}^\dag f_{h+1}(s_h^\dag,a_h)}\\
\leq&\int_{s_{h+1}^\dag}\TT(s_{h+1}^\dag|s_h^\dag,a_h)\int_{a_{h+1}}\abs{\underline{\pi}_{h+1}(a_{h+1}|s_{h+1}^\dag)-\pi_{h+1}(a_{h+1}|s_{h+1}^\dag)}\int_z(f_{h+1}(z-r_h|s_{h+1}^\dag,a_{h+1})\\
\leq&\max_{s_{h+1}^\dag}\norm{\underline{\pi}_{h+1}(\cdot|s_{h+1}^\dag)-\pi_{h+1}(\cdot|s_{h+1}^\dag)}_1\\
\leq&\epsilon\,.
\end{align*}

We can also bound the second term as:

\begin{align*}
&\TV\bracket{\mathcal{T}_{h,\underline{\mpi}}^\dag \underline{f}_{h+1}(s_h^\dag,a_h)||\cT_{h,\underline{\mpi}^\dag} f_{h+1}(s_h^\dag,a_h)}\\
\leq&\int_{s_{h+1}^\dag}\TT(s_{h+1}^\dag|s_h^\dag,a_h)\int_{a_{h+1}}\underline{\pi}_{h+1}(a_{h+1}|s_{h+1}^\dag)\int_z \abs{\bracket{f_{h+1}-\underline{f}_{h+1}}(z-r_h|s_{h+1}^\dag,a_{h+1})}\\
\leq & \max_{s_{h+1}^\dag,a_{h+1}}\norm{f_{h+1}(s_{h+1}^\dag,a_{h+1})-\underline{f}_{h+1}(s_{h+1}^\dag,a_{h+1})}_1\\
\leq& \mathcal{O}(\epsilon)\,,
\end{align*}

since $\underline{f}\leq (1+2\epsilon)f^\downarrow\leq(1+2\epsilon)f$ and $\underline{f}\geq f^\downarrow$ point wise. Then we can have $\norm{\underline{f}-f}_1=\int_z \abs{f(z)-\underline{f}(z)}\leq \int_z \max\sets{2\epsilon f, f-f^\downarrow}\leq \int_z 2\epsilon f+\int_z(f-f^\downarrow)\leq 3\epsilon$.

We can conclude that: $\sum_{i=1}^k\EE_{s_h^\dag,a_h\sim\mu^{\mpi^i}_h}\mbracket{\TV^2\bracket{\mathcal{T}_{h,{\mpi}}^\dag f_{h+1}(s_h^\dag,a_h)||\mathcal{T}_{h,\underline{\mpi}}^\dag \underline{f}_{h+1}(s_h^\dag,a_h)}}\leq\mathcal{O}(k\epsilon)$. Thus,

\begin{align*}
\sum_{i=1}^k\EE_{s_h^\dag,a_h\sim\mu^{\mpi^i}_h}\mbracket{\TV^2\bracket{f_h(s_h^\dag,a_h)||\mathcal{T}_{h,\mpi}^\dag f_{h+1}(s_h^\dag,a_h)}}\leq \mathcal{O}(\gamma^{\MLE})\,.
\end{align*}

Next, we present the distribution difference lemma for our model free analysis.

We can bound the distance between the CDFs of the return by the bellman error of each step as:

\begin{align*}
\norm{F_{Z}-F_{Z^{\mpi}}}_{\infty}\leq \sum_{h=1}^H \EE_{(s_h^\dag,a_h)\sim\mu^{\mpi}_h}\mbracket{\norm{f_h(s_h^\dag,a_h)-\mathcal{T}_{h,\mpi}^\dag f_{h+1}(s_h^\dag,a_h)}_1}\,.
\end{align*}

We begin by induction. By triangle inequality:

\begin{align*}
\norm{f_h(s_h^\dag,a_h)-f_h^{\mpi}(s_h^\dag,a_h)}_1\leq \norm{f_h(s_h^\dag,a_h)-\mathcal{T}_{h,\mpi}^\dag f_{h+1}(s_h^\dag,a_h)}_1+\norm{\mathcal{T}_{h,\mpi}^\dag f_{h+1}(s_h^\dag,a_h)-\mathcal{T}_{h,\mpi}^\dag f_{h+1}^{\mpi}(s_h^\dag,a_h)}_1\,.
\end{align*}

We derive a recursion for the second term.

\begin{align*}
&\norm{\mathcal{T}_{h,\mpi}^\dag f_{h+1}(s_h^\dag,a_h)-\mathcal{T}_{h,\mpi}^\dag f_{h+1}^{\mpi}(s_h^\dag,a_h)}_1\\
=&\int_z\abs{\int_{s_{h+1}^\dag,a_h} \TT(s_{h+1}^\dag|s_h^\dag,a_h)\pi_{h+1}(a_{h+1}|s_{h+1}^\dag)\bracket{f_{h+1}(z-r_h|s_{h+1}^\dag,a_{h+1})-f_{h+1}^{\mpi}(z-r_h|s_{h+1}^\dag,a_{h+1})} }\\
\leq& \int_{s_{h+1}^\dag,a_h} \TT(s_{h+1}^\dag|s_h^\dag,a_h)\pi_{h+1}(a_{h+1}|s_{h+1}^\dag) \int_z \abs{f_{h+1}(z-r_h|s_{h+1}^\dag,a_{h+1})-f_{h+1}^{\mpi}(z-r_h|s_{h+1}^\dag,a_{h+1})} \\
=& \EE_{(s_{h+1}^\dag,a_{h+1})\sim(s_h^\dag,a_h)}\norm{f_{h+1}(s_{h+1}^\dag,a_{h+1})-f_{h+1}^\mpi(s_{h+1}^\dag,a_{h+1})}_1\,.
\end{align*}

Then we have:

\begin{align*}
& \EE_{(s_h^\dag,a_h)\sim \mu^{\mpi}_h}\mbracket{\norm{f_h(s_h^\dag,a_h)-f_h^{\mpi}(s_h^\dag,a_h)}_1}\\
\leq& \EE_{(s_h^\dag,a_h)\sim\mu^{\mpi}_h}\mbracket{\norm{f_h(s_h^\dag,a_h)-\mathcal{T}_{h,\mpi}^\dag f_{h+1}(s_h^\dag,a_h)}_1}\\
&+\int_{s_h^\dag,a_h} \mu^{\mpi}_h(s_h^\dag,a_h)\int_{s_{h+1}^\dag,a_{h+1}} \mu^{\mpi}(s_{h+1}^\dag,a_{h+1}|s_h^\dag,a_h)\norm{f_{h+1}(s_{h+1}^\dag,a_{h+1})-f_{h+1}^\mpi(s_{h+1}^\dag,a_{h+1})}_1\\
\leq &\EE_{(s_h^\dag,a_h)\sim\mu^{\mpi}_h}\mbracket{\norm{f_h(s_h^\dag,a_h)-\mathcal{T}_{h,\mpi}^\dag f_{h+1}(s_h^\dag,a_h)}_1}+ \EE_{(s_{h+1}^\dag,a_{h+1})\sim\mu^{\mpi}_{h+1}}\mbracket{\norm{f_{h+1}(s_{h+1}^\dag,a_{h+1})-f_{h+1}^\mpi(s_{h+1}^\dag,a_{h+1})}_1}\,.
\end{align*}

Using the definition of $F_Z(x)=\int_{s_1^\dag,a_1}\mu_1^\mpi(s_1^\dag,a_1) \int_{z\leq x}f_1(z|s_1^\dag,a_1)$, we have:

\begin{align*}
\leq&\EE_{s_1^\dag,a_1\sim\mu^{\mpi}_1}\mbracket{\norm{f_1\bracket{s_1^\dag,a_1}-f_1^\mpi(s_1^\dag,a_1)}_1}\\
\leq& \sum_{h=1}^H \EE_{(s_h^\dag,a_h)\sim\mu^{\mpi}_h}\mbracket{\norm{f_h(s_h^\dag,a_h)-\mathcal{T}_{h,\mpi}^\dag f_{h+1}(s_h^\dag,a_h)}_1}\,.
\end{align*}

Combining the elliptical potential condition for low bellman eluder dimension (Lemma <ref>) and the concentration result Lemma <ref>, we have for all $h\in[H]$ and any $\hat{\bm{f}}^k\in\widehat{\cF}_{k,\mpi^k}^{\MLE}$:

\begin{align*}
\sum_{k=1}^K \EE_{s_h^\dag,a_h\sim\mu^{\mpi^k}_h}\mbracket{\TV\bracket{\hat{f}^k_h(s_h^\dag,a_h)||\mathcal{T}_{h,\mpi^k}^\dag \hat{f}_{h+1}^k(s_h^\dag,a_h)}}\leq \mathcal{O}(\sqrt{\operatorname{d_{BE}}\gamma^{\MLE} K})\,.
\end{align*}

Here we invoke Lemma <ref> by setting $\gX:\cS^\dag\times\cA$, $\Phi: \TV\bracket{f_h(s^\dag,a)||\cT_{h,\mpi}^\dag f_{h+1}(s^\dag,a)}$  for all $(f,\mpi)\in\cF\times\Pi^\dag$. $\cD:\mu^{\mpi},\,\mpi\in\Pi^\dag$ is the family of all the visitation measures defined on $(s^\dag,a)$.

Thus, using Lemma <ref>, we have:

\begin{align*}
\sum_{k=1}^K \norm{F_{\widehat{Z}^k}-F_{Z^{\mpi^k}}}_\infty\leq \mathcal{O}(H\sqrt{\operatorname{d_{BE}}\gamma^{\MLE} K})\,.
\end{align*}

§ AUXILIARY LEMMAS

Let $\gX$ be a feature space and $\gY$ be the output space. Given a dataset $D=\sets{(x_i,y_i)}_{i=1}^n$ which is collected from a martingale process: $x_i\sim\gD_i(x_{1:i-1},y_{1:i-1})$, and $y_i\sim p(\cdot|x_i)$. Given a function set $\gF:\gX\times\gY\rightarrow\RR$, we have the real conditional distribution $f^*(x,y)=p(y|x)\in\gF$. Then, there exists a constant $c$, for any $\delta>0$, with probability at least $1-\delta$, we have:

\begin{align*}
\sum_{i=1}^n \EE_{x\sim\cD_i}\mbracket{\TV\bracket{f(x,\cdot)||f^*(x,\cdot)}^2}\leq c\bracket{\sum_{i=1}^n \log\bracket{\frac{f^*(x_i,y_i)}{f(x_i,y_i)}}+\log \bracket{\cN_{[\cdot]}(\epsilon,\cF,\norm{\cdot}_1/)\delta}+n\epsilon}
\end{align*}

let $(X_p,Y_p)_{p=1,2\cdots}$ be a set of random variables, $X_p\in\gX$ for some measurable set $\gX$ and $Y_p\in\RR$. Let $\cF$ be a set of real valued measurable function with domain $\mathcal{X}$. Let $\mathbb{F}=\bracket{\mathbb{F}_p}_{p=0,1,2\cdots}$ be a filtration such that for all $p\geq1$, we have $\bracket{X_1,Y_1,\cdots,X_{p-1},Y_{p-1},X_p}$ is $\mathbb{F}_{p-1}$ measurable, and such that there exists some function $f_*\in\cF$ such that $\EE[Y_p|\mathbb{F}_{p-1}]=f_*(X_p)$ for all $p\geq1$. Let $\hat{f_t}=\argmin_{f\in\cF}\sum_{p=1}^{t}\bracket{f(X_p)-Y_p}^2$. Let $\cN(\cF,\alpha)$ be the $\alpha$-covering number of set $\cF$ under $\norm{\cdot}_{\infty}$ metric at scale $\alpha$. Define $\dist_t(f||f_t)=\sum_{p=1}^t\bracket{f(X_p)-f_t(X_p)}^2$.

If the functions in $\cF$ are bounded by some constant $C>0$. Assume that for each $p\geq1$, $\bracket{Y_p-f_*(X_p)}$ is conditionally $\sigma$-sub-gaussian given $\mathbb{F}_{p-1}$. Then, for any $\alpha>0$, with probability $1-\delta$ for all $t\geq1$, we have:

\begin{align*}
\dist\bracket{f_*||f_t}\leq 8\sigma^2\log\bracket{\cN(\cF,\alpha)/\delta}+4t\alpha\bracket{C+\sqrt{\sigma^2 \log\bracket{t(t+1)/\delta}}}
\end{align*}

§ LINEAR CVAR

Similarly to [36, 26], the general algorithms provided for general version space are information theoretic, which means that they cannot be implemented efficiently in general. This is because we consider the general risk measure LRM and the general function approximation settings. However, when specified to the CVaR risk measure under the discretized linear MDP, a distributional extension of natural linear MDP [27], we can design and implement an efficient model-free algorithm that achieves sub-linear regret.

An augmented MDP $\cM^\dag$ is a discretized linear MDP with feature map $\phi:(s,a)\rightarrow \RR^d$ and an uniform grid of $M$ points $\sets{z_i}_{i=1}^M$, if for any $h\in[H]$, $(r_h,y_h)\in\sets{z_i}_{i=1}^M$, and there exists unknown measures $\mu_h:\cS\rightarrow\RR^d$ and $\theta_h:\sets{z_i}_{i=1}^M\rightarrow\RR^d$, such that:

\begin{align*}&\PP_h(s_{h+1}|s_h,a_h)=\phi(s_h,a_h)^\top\mu_h(s_{h+1})\\&\RR_h(z_i|s_h,a_h)=\phi(s_h,a_h)^\top\theta_h(z_i)~,
\end{align*}

for all $s_h,a_h,s_{h+1},z_i$.

This discretized linear MDP is the natural extension of the linear MDP assumption in [27], where we consider the discretized distributional reward instead of determined reward. and we generalize the linear expected reward to its distributional counterpart. Another important ingredient in our definition is the discretized reward space, which is commonly used in practice. C51 and Rainbow [8, 22] both set $M=51$ and achieved empirical success in Atari games. We need the discretized reward space mainly to bound the covering number of the value distribution, similar to [47].

In a discretized linear MDP, we have that the distribution function have a quadratic structure: $f_h^\pi(z_i|(s_h,y_h),a_h)=\phi(s_h,a_h)W_h^\pi(z_i,y_h)\phi(s_h,a_h)$. Thus, we can use linear regression in estimating statistical functionals of $Z_h^\mpi(s_h^\dag,a_h)$. We present our regret bound as follows:

If MDP $\cM^\dag$ is a discretized linear MDP satisfying Definition <ref>, and the risk measure $\rho(Z^\mpi)=CVaR_{\tau}(Z^\pi)=\max_b\sets{b-\tau^{-1}\EE[(b-Z^\pi)^+]}$, we can bound the regret as:

\begin{align*}
\operatorname{Regret}(K)\leq \mathcal{\Tilde{O}}(\tau^{-1}d^3H^2\sqrt{MK})
\end{align*}

We highlight that this is the first efficient model free algorithm for linear MDPs, and the $\sqrt{M}$ dependency is due to the covering number of the value distribution class which also appears in other model free distributional RL algorithms, such as [47].

Define $\operatorname{CVaR}_\tau(Z^\pi)=\argmax_{b}\sets{b-\tau^{-1}\EE[(b-Z^\pi)^{+}]}$. We can define the statistical functionals $Q$ and $V$ as: $Q^\pi_h(s_h,y_h,a_h)=\sum_{z_i}f_h^\pi(z_i|s_h,y_h,a_h)[(-z_i-y_h)^{+}]$, and $V^\pi_h(s_h,y_h)=\sum_{z_i}f_h^\pi(z_i|s_h,y_h,\pi(s_h,y_h))[(-z_i-y_h)^{+}]$. Then we can write the CVaR objective as $\operatorname{CVaR}_\tau(Z^\pi)=\argmax_{b\in[0,H]}\sets{b-\tau^{-1}V_1^{\pi}(s_1,-b)}$.

Define recursively $\pi^*_h(s_h,y_h)=\argmin_a Q_h^{\pi^*}(s_h,y_h,a_h)$, then by [46] we have for any $b\in[0,H]$,

\begin{align*}
\end{align*}

We denote $\operatorname{CVaR}_{\tau}^*=\max_{\pi\in\Pi^\dag}\operatorname{CVaR}_{\tau}^\pi=b^*-\tau^{-1}V_1^*(s_1,-b^*)$

§.§ Linear Augmented MDPs

When consider the linear function approximation (Definition <ref>), we can also linearize the augmented MDPs. We ahve

where $\psi(s_h,a_h)$ and $\chi(r_h,s_{h+1})$ are the flattened versions of $\phi(s_h,a_h)\phi(s_h,a_h)^\top$ and $\theta(r_h)\mu\bracket{s_{h+1}}^\top$. Also, we assume the reward space is discretized into $M$ points $z_1\cdots z_M$ such that for all $h\in[H]$, $y_h\in \sets{z_i}_{i=1}^M$. We have that $-H\leq z_1\leq z_2\cdots z_M\leq H$.
We highlight that this discretization is standard in practice as in C51 [8] and Rainbow [22]. Here we need this assumption to bound the complexity of the function class.

§.§ Linear Completeness

Because the density function of the rewards satisfy the distributional bellman equation:

\begin{align*}
\end{align*}

We have the statistical functionals satisfy the augmented bellman equation:

\begin{align*}
=&\bracket{\TT_h V_{h+1}}(s_h,y_h,a_h)
\end{align*}

Where we denote the augmented transition operating on a function $V$ as:

\begin{align*}
\TT_hV_{h+1}(s_h,y_h,a_h)=\sum_{s_{h+1},r_h}\TT_h(s_{h+1},r_h|s_h,a_h)V_{h+1}(s_{h+1},y_h+r_h)
\end{align*}

Since $f_h^\pi(z_i|s_h,y_h,a_h)=\phi(s_h,a_h)^\top W_h^\pi(y_h,z_i)\phi(s_h,a_h)$, we have $Q_h^{\pi}(s_h,y_h,a_h)=\psi(s_h,a_h)^\top w_h^\pi(y_h)$. Where $w_h^\pi[i\cdot d+j]=\sum_{z_i}W_h^\pi(y_h,z_i)[i,j][(-z_i-y_h)^{+}]$

§.§ Algorithm

In this section, we present our computationally efficient algorithm . Notice that we present Upper Confidence Bound Value Iteration (UCV-VI) in this algorithms instead of the general optimistic planning used in previous frameworks, due to its computation-tractable property.

Input: Features $\Psi:\psi(s,a)$ Bonus $b_h^k(s,a)=\beta \norm{\psi(s,a)}_{\Lambda_{k,h}^{-1}}$ with $\beta=c_\beta d^2H\sqrt{M\log(d^2UMK/\delta)}$
$k=1$ to $K$
Set $V_{H+1}(s_{H+1},y)=[(-y)^+]$ for all $y\in\sets{z_i}_{i=1}^M$
$h=H$ to $1$

\begin{align*}
&\Lambda_{k,h}=\sum_{i=1}^{k-1}\psi(s_h^i,a_h^k)\psi(s_h^i,a_h^i)^\top+\lambda \mathbf{I}\\
&\pi^k(s_h,y_h)=\argmin_a Q_h^k(s_h,y_h,a_h)\\
\end{align*}

\begin{equation*}
\end{equation*}

Start at state $s_1^\dag=(s_1,-\lambda^k)$ and execute policy $\pi^k$ in the augmented MDP $\cM^\dag$, collect information $\sets{(s_h^k,a_h^k,r_h^k)}_{h\in[H]}$

§.§ Concentration and Covering

In this section we provide the concentration and covering arguments needed for our linear analysis.

Let $\sets{x_\tau}_{\tau=1}^{\infty}$ be a stochastic process on domain $\gX$ with corresponding filtration $\sets{\cF_\tau}$. Let $\sets{\psi_\tau}$ be an $\RR^{d^2}$ valued stochastic process  stochastic process such that $\psi_\tau\in\cF_{\tau-1}$ and $\norm{\psi_\tau}\leq 1$. Let $\Lambda_k=\sum_{\tau=1}^{k-1}\psi_\tau\psi_\tau^\top$. Then, for any $\delta>0$, with probability at least $1-\delta$, for all $k>0$, and any $V\in\cV$ such that $\abs{\sup_x V(x)-\inf_x V(x)}\leq H$, we have:

\begin{align*}
\norm{\sum_{\tau=1}^{k-1}\psi_\tau\sets{V(x_\tau)-\EE[V(x_\tau|\cF_{\tau-1}]}}_{\Lambda_k^{-1}}^2\leq H^2 \mbracket{\frac{d}{2}\log\bracket{\frac{k+\lambda}{\lambda}}+\log\bracket{\cN_C(\epsilon,\cV,\norm{\cdot}_{\infty}}/\delta}+\frac{8k^2\epsilon^2}{\lambda}~.
\end{align*}

Let $\gV$ denote a class of functions mapping from domain $\gX$ to $\RR$ with the following parametric form:

\begin{align*}
V(\cdot)=w^\top\phi(\cdot)+\beta \norm{\phi(\cdot)}_{\Lambda^{-1}}
\end{align*}

where $\phi\in\RR^d$ are features on domain $\gX$.
The parameters satisfy $\norm{w}\leq L$, $ \lambda_{\min}(\Lambda)\geq \lambda$, $\beta\in[0,B]$, and $\norm{\phi(\cdot)}\leq1$. Then, the log covering number can be bounded as:

\begin{align*}
\log\bracket{(\cN_C(\epsilon,\gV,\norm{\cdot}_\infty)}\leq d \log (1+4L/\epsilon)+d^2\log(1+8d^{1/2}B^2/\lambda\epsilon^2)
\end{align*}

For $h\in[H]$ and $k\in[K]$, if $V_{h+1}^k(s_{h+1},y_{h+1})\leq H$, then $\norm{w_h^k}_2\leq Hd\sqrt{k/\lambda}$

For any $h\in[H]$, $\norm{w_h^\pi}\leq 2Hd$

If $V_{h+1}^k\leq H$ for any $k\in[K]$, then there exists a constant $C$, for any $y_h\in\sets{z_i}_{i=1}^M$, $\delta>0$, with probability $1-\delta$, we can bound the self normalized martingale process as:

\begin{equation*}
\norm{\sum_{i=1}^{k-1}\psi_{i,h}\bracket{V_{h}^k(s_{h}^i,y_{h-1}+r_{h-1}^i)-\TT_{h-1} V_{h}^k(s_{h-1}^i,y_{h-1},a_{h-1}^i)}}_{\Lambda_{k,h}^{-1}}^2\leq C d^4 H^2M{\log(c_\beta d^2 KM/\delta)}
\end{equation*}

According to Lemma <ref>, we have that $\norm{w_h^k}_2\leq Hd\sqrt{k/\lambda}$. According to Lemma <ref>, we have for a fixed $z_i$ the covering number of the function class $V_h^k(s_h,z_i)$ can be bounded as (notice that $\psi\in\RR^{d^2}$):

\begin{align*}
\log\bracket{(\cN_C(\epsilon,\gV(\cdot,z_i),\norm{\cdot}_\infty)}\leq d^2 \log (1+4Hd\lambda/k\epsilon)+d^4\log(1+8d\beta^2/\lambda\epsilon^2)
\end{align*}

Then we can bound the entire function class by regarding each $V(\cdot,z_i)$ an individual function. Thus, the total covering number can be bounded as:

\begin{align*}
\log\bracket{(\cN_C(\epsilon,\gV,\norm{\cdot}_\infty)}\leq  Md^2 \log (1+4Hd\lambda/k\epsilon)+d^4\log(1+8d\beta^2/\lambda\epsilon^2)
\end{align*}

Thus, we can apply Lemma <ref>, and we have the result by choosing $\epsilon=1/k$, $\lambda=1$ and $\beta=c_\beta d^4HM\sqrt{\log(d^2UMK/\delta)}$.

§.§ Proof of Theorem <ref>

Define the event
<ref> at $h\in[H]$ and $k\in[K]$ as:

\begin{equation}\label{eq:linear optimism at h}\tag{Optimism}
Q_{h}^{k}(s_{h},y_{h},a_{h})\leq Q_{h}^{*}(s_{h},y_{h},a_{h})
\end{equation}

Also, define

\begin{align*}
\psi(s,a)w_h^k(y)-Q_h^\pi(s,y,a)-\TT_h (V_{h+1}^k-V_{h+1}^\pi)(s,a,y)=\Delta_h^k(s,a)
\end{align*}

Then we have the following concentration result:

If <ref> holds at $h\in[H]$ and $k\in[K]$. We have that $\Delta_{h}^k(s,a)\leq\beta \norm{\psi(s,a)}_{\Lambda_{k,h}^{-1}}$.

Since we have:

\begin{align*}
\psi(s_h^k,a_h^k)\TT_h V_{h+1}^\pi(s_h^k,a_h^k,y_h)=& \sum_{s_{h+1},r_h}\psi(s_h^k,a_h^k)\psi(s_h^k,a_h^k)^\top\chi(s_{h+1},r_h)V_{h+1}^\pi(s_{h+1},y_h+r_h)=\psi(s_h^k,a_h^k)\psi(s_h^k,a_h^k)^\top w_h^\pi(y_h)
\end{align*}

We bound the bias as:

\begin{align*}
w_h^k(y_h)-w_h^\pi(y_h)=&\Lambda_{k,h}^{-1}\sets{-\lambda w_h^\pi+\mbracket{\sum_{i=1}^{k-1}\psi_{i,h}\bracket{V_{h+1}^k(s_{h+1}^i,y_h+r_h^i)-\TT_h V_{h+1}^\pi(s_h^i,a_h^i,y_h)}}}\\
=&\underbrace{-\lambda \Lambda_{k,h}^{-1}w_h^{\pi}}_{q_1}+\underbrace{\Lambda_{k,h}^{-1}\sum_{i=1}^k\psi_{i,h}\bracket{V_{h+1}^k(s_{h+1}^i,y_h+r_h^i)-\TT_h V_{h+1}^k(s_h^i,a_h^i,y_h)}}_{q_2}\\
&+\underbrace{\Lambda_{k,h}^{-1}\sum_{i=1}^k \psi_{i,h}\TT_h(V_{h+1}^k-V_{h+1}^\pi)(s_h^i,a_h^i,y_h)}_{q_3}
\end{align*}

We have according to Lemma <ref>:

\begin{align*}
\abs{\psi(s_h,a_h)^\top q_1}\leq \sqrt{\lambda}\norm{w_h^\pi}\norm{\psi(s_h,a_h)}_{\Lambda_{k,h}^{-1}}\leq 2\sqrt{\lambda}Hd\norm{\psi(s_h,a_h)}_{\Lambda_{k,h}^{-1}}
\end{align*}

For the second term, on event <ref>, we have that $V_{h+1}^k\leq V_{h+1}^*\leq H$. Apply Lemma <ref> we have:

\begin{align*}
\abs{\psi(s_h,a_h)^{\top}q_2}\leq  C d^2 H\sqrt{M\log(c_\beta d^2 KM/\delta)}\norm{\psi(s_h,a_h)}_{\Lambda_{k,h}^{-1}}
\end{align*}

For the third term, we have:

\begin{align*}
{\psi(s_h,a_h)^\top q_3}=&{\psi(s_h,a_h)\Lambda_{k,h}^{-1}\sum_{i=1}^k \psi_{i,h}\psi_{i,h}^\top \sum_{s_{h+1},r_h}\chi(s_{h+1},r_h)(V_{h+1}^k-V_{h+1}^\pi)(s_{h+1},y_h+r_h)}\\
=&\underbrace{\psi(s_h,a_h)^\top\sum_{s_{h+1},r_h}\chi(s_{h+1},r_h)(V_{h+1}^k-V_{h+1}^\pi)(s_{h+1},y_h+r_h) }_{p_1}\\
\end{align*}

\begin{align*}
p_1= \TT_h(V_{h+1}^k-V_{h+1}^\pi)(s_h,a_h)
\end{align*}

\begin{align*}
\abs{p_2}\leq \sqrt{\lambda}d^2U(H-h+1))\norm{\psi(s_h,a_h)}_{\Lambda_{k,h}^{-1}}
\end{align*}

Thus, by choosing the appropriate constants, we have:

\begin{align*}
\abs{\psi(s_h,a_h)^\top w_h^k(y_h)-Q_h^\pi(s_h,y_h,a_h)-\TT_h(V_{h+1}^k-V_{h+1}^\pi)(s_h,a_h,y_h)}\leq \beta \norm{\psi(s_h,a_h)}_{\Lambda_{k,h}^{-1}}
\end{align*}

For any $h\in[H]$ and $k\in[K]$, event <ref> holds

We prove the Lemma via induction. For any $k\in[K]$, <ref> holds at step $H$. Suppose that <ref> holds at step $h+1$, then we have from Lemma <ref>,

\begin{align*}
=& Q_h^\pi(s_h,y_h,a_h)+\TT_h(V_{h+1}^k-V_{h+1}^\pi)(s_h,a_h,y_h)+\Delta_h^k(s_h,a_h)\\
\leq Q_h^\pi(s_h,y_h,a_h)
\end{align*}

So <ref> holds at step $h$.

On event <ref> for all $h\in[H]$ and $k\in[K]$, we have that:

\begin{align*}
b_k-\tau^{-1}V_1^k(s_1,-b_k)\geq b^*-\tau^{-1}V_1^k(s_1,-b^*)\geq b^*-\tau^{-1}V_1^*(s_1,-b^*)=\operatorname{CVaR}_\tau^*
\end{align*}

Also, we have that $\operatorname{CVaR}_\tau^{\pi^k}=\argmax_b\sets{b-\tau^{-1}V_1^{\pi^k}(s_1,-b)}\geq b_k-\tau^{-1}V_1^{\pi^k}(s_1,-b_k)$
Thus, the total regret can be bounded as:

\begin{align*}
\sum_{k=1}^K \operatorname{CVaR}_{\tau}^*-\operatorname{CVaR}_{\tau}^{\pi^k}\leq \sum_{k=1}^K
\tau^{-1}\bracket{V_1^{\pi^k}(s_1,-b_k)-V_1^k(s_1,-b_k)}    \end{align*}

Define $\delta_h^k=V_h^{\pi^k}(s_h^k,y_h^k)-V_h^k(s_h^k,y_h^k)$, and $\zeta_{h+1}^k=\EE[\delta_{h+1}^k|s_h^k,y_h^k,a_h^k]-\delta_{h+1}^k$.
According to Lemma <ref>, we have that:

\begin{align*}
\delta_h^k\leq\delta_{h+1}^k+\zeta_{h+1}^k+2\beta\norm{\psi_{k,h}}_{\Lambda_{k,h}^{-1}}
\end{align*}

Then we have:

\begin{align*}
\sum_{k=1}^K V_1^{\pi^k}(s_1,-b_k)-V_1^k(s_1,-b_k)\leq\sum_{h=1}^H \sum_{k=1}^K\zeta_h^k+2\sum_{h=1}^H\sum_{k=1}^K\beta  \norm{\psi_{k,h}}_{\Lambda_{k,h}^{-1}}
\end{align*}

Since we have event <ref>, we have $\abs{\delta_h^k}\leq H$. Thus, using the Hoeffding inequality, we have:

\begin{align*}
\sum_{h=1}^H\sum_{k=1}^K \zeta_{h}^k \leq \mathcal{O}\bracket{H\sqrt{K\log(KH/\delta)}}
\end{align*}

Using the standard linear elliptical potential lemma, we obtain:

\begin{align*}
\sum_{h=1}^H\sum_{k=1}^K \norm{\psi_{k,h}}_{\Lambda_{k,h}^{-1}}\leq \mathcal{O}(Hd\sqrt{K\log(KH/\delta)})
\end{align*}

Then we have the regret bounded as:

\begin{align*}
\sum_{k=1}^K \operatorname{CVaR}_{\tau}^*-\operatorname{CVaR}_{\tau}^{\pi^k}\leq{\mathcal{\Tilde{O}}}\bracket{\tau^{-1}H^2d^3\sqrt{MK}}
\end{align*}

§.§ Numerical Experiment Results

In this section we provide the details of our numerical experiments. We construct a zero-mean MDP where the expected return for all the state-action pairs are $0$, thus risk-neutral algorithms such as LSVI-UCB of [27] will learn nothing. We also compare our results with the optimistic MDP algorithm of [6]. For simplicity we constructed a toy MDP with $S=3$, $A=2$, $d=2$, $H=6$, $M=3$. The results are in Figure <ref>. From the figures we can see that the risk-neutral algorithm LSVI-UCB fails to learn anything, while the model-based algorithm of [6] converges much slower than ours. Also, with smaller $\tau$ we have larger regret, which is consistent with previous analysis.

Comparison for different algorithms for the CVaR objective $\operatorname{CVaR}_{\tau}$ under different risk parameter $\tau$.


# Xel-FPGAs: An End-to-End Automated Exploration Framework for Approximate
Accelerators in FPGA-Based Systems

Bharath Srinivas Prabakaran1,4, Vojtech Mrazek2,4, Zdenek Vasicek 2, Lukas
Sekanina2, Muhammad Shafique3 4These two authors have contributed to this work
equally. 1Institute of Computer Engineering, Technische Universität Wien (TU
Wien), Austria
<EMAIL_ADDRESS>2Faculty of Information Technology,
IT4Innovations Centre of Excellence, Brno University of Technology, Czech
Republic
{mrazek, vasicek<EMAIL_ADDRESS>3eBrain Lab, Division of
Engineering, New York University Abu Dhabi (NYUAD), United Arab Emirates (UAE)
<EMAIL_ADDRESS>
###### Abstract

Generation and exploration of approximate circuits and accelerators has been a
prominent research domain to achieve energy-efficiency and/or performance
improvements. This research has predominantly focused on ASICs, while not
achieving similar gains when deployed for FPGA-based accelerator systems, due
to the inherent architectural differences between the two. In this work, we
propose a novel framework, Xel-FPGAs, which leverages statistical or machine
learning models to effectively explore the architecture-space of state-of-the-
art ASIC-based approximate circuits to cater them for FPGA-based systems given
a simple RTL description of the target application. We have also evaluated the
scalability of our framework on a multi-stage application using a hierarchical
search strategy. The Xel-FPGAs framework is capable of reducing the
exploration time by up to $95\%$, when compared to the default synthesis,
place, and route approaches, while identifying an improved set of Pareto-
optimal designs for a given application, when compared to the state-of-the-
art. The complete framework is open-source and available online at
https://github.com/ehw-fit/xel-fpgas.

###### Index Terms:

Approximate Computing, Accelerator, FPGA, ASIC, Arithmetic Units, Regression,
Models, Synthesis.

##  I Introduction

Field Programmable Gate Arrays (FPGAs) act as a lucrative computing platform
in various applications, especially ones that need to be updated and/or
modified on-the-go, due to their ability of partial run-time reconfiguration,
i.e., field-programmability. Since their introduction in $1984$ [1], FPGAs
have gained a lot of popularity due to their decreased time-to-market and
lower prototype costs, when compared to Application Specific Integrated
Circuits (ASICs). Besides their deployment in embedded use-cases and cyber-
physical systems, FPGAs are also deployed in servers and high-performance
computing clusters to offer on-demand performance acceleration to complex
compute-intensive algorithms [2]. These reconfigurable platforms include hard
IPs (IC realization) of low-power ARM processor cores and commonly used
hardware accelerators like video codecs, besides the soft IP cores, which are
used to realize the custom accelerators [3]. Due to these wide-ranged
capabilities, FPGAs can act as a Programmable-System-on-chip (PSoC) for
embedded use-cases or as a service for on-demand compute acceleration.
However, FPGAs tend to be more power hungry, with decreased throughput, when
compared to ASICs.

Besides the use of conventional power-reduction approaches, like Dynamic
Voltage and Frequency Scaling (DVFS) [4] and clock-gating [5], Approximate
Computing appears to be quite suitable for increasing the energy-efficiency of
a system. The primary notion behind approximating the computational units of a
system, stems from the fact that these systems and their algorithms are
inherently error-resilient, i.e., they do not lead to “significant”
degradation in output quality [6]. This inherent error-resilience property is
exhibited by applications encompassing several domains like recognition,
mining, and synthesis, to include data analysis [7], audio, image, or video
processing [8], speech and image recognition [9], etc. The resilience
exhibited by these applications can be attributed to four key factors, namely,
(1) redundant input data,  (2) implicit error attenuation by the deployed
algorithms,  (3) non-distinguishable, user-level, differences in the output
quality, and  (4) lack of a single golden output.  These factors have been
exploited widely in various research works, spanning academia and industry,
focusing on both hardware and software to obtain either power, latency, or
energy benefits [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22].

The state-of-the-art tends to primarily focus on generating approximate
circuits (AC) and/or accelerators to obtain power, latency, and/or energy
benefits for ASIC-based systems. Several past works have illustrated that the
approximate computing principles and techniques developed for ASIC-based
systems tend to offer dissimilar benefits when realized for FPGA-based systems
[12, 13, 14, 15, 16, 17, 18, 19, 20]. For instance, approximate arithmetic
units developed for ASICs can achieve up to $70\%$ energy reduction when
synthesized for ASIC platforms; the same designs offer minimal savings, or at
times negative savings, when synthesized for FPGA platforms [12]. Ullah et al.
[14] proposed a methodology for systematically generating a wide-range of
approximate multiplier architectures that offer similar benefits when deployed
on FPGA platforms. The methodology also includes a small-scale design space
exploration stage to identify approximate accelerator designs suitable for the
application. On the other hand, Prabakaran et al. [15] illustrate the benefits
of exploring the vast architecture-space of approximate ASIC-based arithmetic
circuits to identify a set of Pareto-optimal111By Pareto-optimal, we mean the
best solution obtained by a given method. designs suitable for the target FPGA
platform. This is subsequently leveraged by a large-scale automated
approximate accelerator exploration strategy to identify near-optimal
approximate accelerators for an application. However, these works do not
holistically explore the large architecture-space of approximate accelerators
that can be constructed using several arithmetic ACs, which might cause the
exploration to overlook certain trade-offs that are clearly situated on the
Pareto front for the application.

###  I-A Motivational Analysis

To further illustrate the differences between ASICs and FPGAs, we explore
approximate accelerator variants of a Gaussian filter composed of nine $8$-bit
multipliers and eight $16$-bit adders. The approximate accelerator variants
are built using the Pareto-optimal ACs present in the evolutionary approximate
arithmetic circuit library [23] and the ACs identified as Pareto-optimal for
FPGAs [15]. The range of approximate accelerators obtained by identifying
exhaustive combinations of these designs are synthesized using the Vivado
$2017.2$ tool-chain for the Xilinx xc7vx485tffg1157-1 FPGA. We also disable
the use of DSP blocks on the FPGA to ensure that the accelerator is mapped
onto the reconfigurable logic. The quality of these variants, in terms of
average Peak Signal to Noise Ratio (PSNR), are determined using their software
models for a set of input signal samples. These experimental results are
illustrated in Fig. 1, from which we make the following key observations:

  1. (1)

Approximate accelerators that are Pareto-optimal for ASIC-based platforms are
not necessarily Pareto-optimal when synthesized for an FPGA, especially when
considering the different hardware parameters, such as power and latency.

  2. (2)

The time required for synthesizing a subset of approximate accelerators,
composed of $1000$ random designs, is $2.8$ days. This overhead can be
attributed primarily to the high-level mapping algorithms of the Vivado tool-
chain.

  3. (3)

Although quite a few research works focus on designing custom-made approximate
circuits and accelerators for FPGAs [13, 18], they are not as beneficial in
achieving similar improvements in hardware costs, when compared to automated
generation and exploration of approximate circuits, as reinforced by [15, 23].

Figure 1: Analysis of the Pareto-optimal Approximate Accelerators. Figure 2:
An Overview of the Xel-FPGAs Framework.

From these observations, we have identified the following research challenges:

  * •

Given the amount of time required to explore the large number of approximate
accelerator variants for a simple Gaussian filter, the time required for
exploring larger applications would be in the order of months or years.

    * –

How can we leverage statistical or machine learning models to efficiently
reduce the exploration time for FPGA-based approximate accelerators?

    * –

Can non-conventional features and tools be used, besides the traditional
hardware requirements obtained from an FPGA synthesis tool, to reduce the
exploration time?

  * •

There is no single tool that can identify approximate accelerator variants for
FPGAs given an application.

    * –

Can we build an end-to-end framework that can generate and explore approximate
accelerator variants for FPGAs given the application’s hardware and software
models?

To address these research challenges, we propose the following novel
contributions:

  * •

An end-to-end automated framework, Xel-FPGAs, for systematically exploring the
architecture-space of FPGA-based approximate accelerator variants given a
target application.

  * •

Our framework employs statistical learning models to traverse the
architecture-space of approximate accelerators for FPGAs without synthesizing
each individual design.

  * •

To further reduce the exploration time of the approximate accelerator
architecture-space, we use the ABC tool [24] for estimating the number of
lookup tables, power, and latency of the accelerator on our target FPGA.

  * •

The complete Xel-FPGAs framework, including the approximate accelerators
discussed in this work, are open-source and available online at
https://github.com/ehw-fit/xel-fpgas, to enable wide-spread adoption and
research in this area.

##  II The Automated Framework

Fig. 2 presents an overview of the proposed Xel-FPGAs framework. The framework
is made of three key stages, namely, Model Training, Architecture Exploration,
and Final Evaluation. In the first stage, we train and test various
statistical learning models that can be used to estimate the parameters of a
given hardware accelerator module. Next, the models are used to estimate the
potential output quality and hardware requirements of all approximate
accelerators in the architecture-space. Finally, the selected accelerator
designs are re-evaluated by simulating and synthesizing their software and
hardware models, respectively.

###  II-A Model Training

We start with a set of ACs that can replace the arithmetic units in the target
algorithms. We consider the wide-ranged evolutionary approximate component
library [23], which is composed of both adders and multipliers, to evaluate
the efficacy of our framework. The use of other adder and multiplier libraries
[14] are orthogonal to the use of this library and can be easily incorporated
into Xel-FPGAs.

The current architecture-space exploration pipeline requires the system
designers to generate all possible approximate accelerator variants and
synthesize them using an FPGA synthesis tool in order to determine their
hardware requirements. Similar exploration strategies need to be carried out
using software models for determining the quality-of-result (QoR) achievable
by these variants. As emphasized earlier, this can lead to an exhaustive
architecture-exploration stage that can potentially last for months or even
years. To circumvent this bottleneck, we propose to train and deploy
statistical learning models that can be used to estimate the hardware
requirements and QoR capabilities of an accelerator variant, thereby
drastically reducing the time needed for synthesizing and evaluating each
individual approximate accelerator. Previous works [15] have evaluated similar
models for approximate circuits and determined that statistical regression
models achieve maximum correlation, in comparison to other approaches like
Stochastic Gradient Descent, Adaptive Boosting, or even Multi-Layer Perceptron
(MLP). We focus primarily on statistical regression models, which are
lightweight, instead of larger machine learning models, to reduce the overall
time required for training the models and exploring the architecture-space.

Figure 3: Overview of the Possible Model Training Pipelines.

Fig. 3 illustrates the model training pipelines explored in our work to
determine the optimal setup that can quickly and efficiently explore the
architecture-space of approximate accelerators. (A) denotes the default
exhaustive exploration approach discussed earlier, which requires each
approximate accelerator to be synthesized by our FPGA tool-chain. Since this
is not feasible, we propose to generate a feature-set (composed of power,
lookup tables, and delay) by synthesizing the arithmetic ACs present in our
library. This feature-set is used to estimate the hardware parameters for a
set of $1000$ random approximate accelerators. These random approximate
accelerators need to be subsequently synthesized using the FPGA tool-chain for
the model to learn the relation between the features and its desired value.
This pipeline is depicted in (B). However, using the Vivado tool-chain for
synthesizing the complete library of ACs is a highly time-consuming process
(see Section IV). Therefore, we propose to replace the Vivado tool-chain with
the ABC tool [24] to estimate the features necessary for the model to learn.
This pipeline is presented in (C). Another variation of this version would be
to include more features that could potentially improve the accuracy of our
statistical model. This variant, illustrated in (D), includes an additional
step wherein the approximate accelerator variants’ features are extracted by
ABC and fed into the model training stage, with the vivado synthesis values of
the approximate accelerators. And finally, to evaluate all possible
combinations, we tested out two additional scenarios; one (i) wherein the
original feature extraction pipeline, with the Vivado tool-chain, stays in
place for the arithmetic ACs, while the ABC tool is used to extract the
approximate accelerators’ features (presented in (E)), and second (ii) where
the arithmetic ACs’ features are not estimated and only the approximate
accelerators’ features, extracted by the ABC tool, are used to train the
statistical model (illustrated in (F)). We have compared the experimental
results of these pipelines and analysed them to determine that configuration
(D) is most suitable for our framework (see Section IV).

###  II-B Architecture Exploration

Now that we have the statistical models trained, for estimating the
application-level QoR and its hardware requirements, we can efficiently
explore the architecture-space of our approximate accelerator variants.
Despite deploying statistical models to reduce the overhead time, we
investigated the use of several state-of-the-art evolutionary algorithms to
further decrease the time needed for exploring the architecture-space by
selecting specific models that need to be evaluated. For example, we can start
with a set of $1000$ random approximate accelerator variants, which are fed to
the ABC tool to determine their features. With this information, we can
estimate the approximate accelerators’ hardware requirements and its QoR. The
$200$ best, Pareto- and near-optimal, designs are subsequently fed back to the
evolutionary algorithms to act as the parent set. Using the parent set, we
obtain the next generation of $1000$ approximate accelerators that are to be
evaluated, from which we select the $200$ best designs that subsequently act
as the parent set for the next generation. This process is repeated until we
have $1000$ approximate accelerator variant generations. We consider the $200$
best approximate accelerator variants obtained from the last generation to be
the final set of Pareto-optimal designs, which we have observed to achieve
optimal results. Further generation and exploration of approximate accelerator
variants achieves sub-optimal results or minimalistic improvements when
compared to the time needed for evaluating these designs. Therefore, we
terminate the exploration stage and move to the final stage of the framework,
where we evaluate its efficacy.

Evolutionary Algorithms: Besides investigating the commonly used NSGA-II [25]
genetic algorithm for exploring our design space, we decided to evaluate two
other approaches based on local search to investigate their applicability on
our design space. Motivated by the results illustrated in [22], we chose to
implement and evaluate the Hill climber algorithm [26], which ultimately
turned out to be highly inefficient for exploring our design space. The
algorithm was able to identify only $1$ candidate solution in each iteration
stage, which is then fed into our statistical models. However, these models
are built in such a way that they make heavy use of vectorization and
parallelism, and the evaluation on a feature-by-feature basis is slow, which
is against the ultimate goal of our framework, i.e., to reduce exploration
time. Hence, we chose to explore the Evolution Strategy (ES) $\mu+\lambda$,
where $\lambda$ new candidate designs are generated from a set of best
individuals ($\mu$). In contrast to NSGA-II, which depends on “genetics” and
complex “crossover” mechanisms to identify the next generation of designs, ES
is only dependent on small random changes in the best individuals. We
iteratively run through each design in a given Pareto-front –under a time
constraint– to obtain new designs, which can be used to extract a new Pareto-
frontier. The designs from this new frontier act as the base for the next
generation of designs.

Our multi-objective evolutionary algorithm is described in Algorithm 1. The
algorithm starts with a set of random parents ($\mu$) using which we generate
$\lambda$ new candidates with random mutations on each parent. The new
candidates’ power and QoR are parallelly estimated by our statistical models
to determine if the design updates the Pareto-frontier $P$. If a new design is
added to $P$, the design becomes a new parent, which can be subsequently used
to generate other new candidates. Otherwise, we iteratively select a new
candidate from the existing Pareto-frontier $P$ in order to prevent getting
stuck.

Algorithm 1 Our Evolutionary Multi-objective Optimization

$RL$ – set of libraries, $RL=\\{RL_{1},RL_{2},\cdots,RL_{n}\\}$, $M_{HW}$ – HW
costs model, $M_{QoR}$ – quality model

Pareto set $P\subseteq RL_{1}\times RL_{2}\times\cdots\times RL_{n}$

function EvolutionaryStrategyMOO($RL$, $M_{QoR}$, $M_{C}$)

$ParentSet\leftarrow$

$\textsc{PickRandom}(RL_{1}\times RL_{2}\times\cdots\times RL_{n})$
$\triangleright$ $\mu$ parents

$P\leftarrow\emptyset$

while $\neg TerminationCondition$ do

for all $Parent\in ParentSet$ do

$C\leftarrow\\{\textsc{Mutate}(Parent),\cdots\\}$ $\triangleright$ $\lambda$
new candidates

$e_{QoR}\leftarrow M_{QoR}(C)$ $\triangleright$ Estimate the quality of
$C_{i}$

$e_{HW}\leftarrow M_{HW}(C)$ $\triangleright$ Estimate the HW costs of $C_{i}$

for all $c_{i}\in C$ do

if ParetoInsert($P,(e_{QoR},e_{HW}),C$) then

$Parent\leftarrow C$ $\triangleright$ Replace current parent

end if

end for

if $Parent\notin C$ then

$Parent\leftarrow\textsc{PickRandomlyFrom}(P)$ $\triangleright$ Prevent
stagnation

end if

end for

end while

return P

end function

###  II-C Final Evaluation

Once the final set of approximate accelerator variants are obtained, we
synthesize and evaluate them to determine their accurate hardware requirements
and QoR values, instead of the estimates obtained from the models. These
values are subsequently used to construct a Pareto-frontier, which acts as the
final set of Pareto-optimal variants for the target application. The system
designer can then analyze the quality-hardware trade-offs of these designs to
deploy the relevant approximate accelerator variant for the application based
on system specifications.

##  III Experimental Setup

The approximate arithmetic circuit library used in this work is open-source
and immediately accessible222https://github.com/ehw-
fit/evoapproxlib,333https://github.com/ehw-fit/approx-fpgas; they include both
the hardware (Verilog) and software models (C++) of the ACs. These can be
combinatorially deployed in the target application to generate different
possible approximate accelerator variants. Without loss of generality, in this
work we consider a sample set of $1000$ random approximate accelerator
variants to train our statistical regression model. Since our framework
follows pipeline (D) from Fig. 3, we perform a component- and accelerator-
level feature-extraction, from the RTL models, while also synthesizing and
mapping the designs, using the Vivado Design Suite 2017.2 for the
xc7vx485tffg1157-1 Xilinx FPGA, to obtain the relevant hardware parameters,
such as power, number of lookup tables, and latency of the approximate
accelerator. We also simulate the application with behavioral models of the
accelerator in C++, using the benchmark sequences discussed in [27], to
determine the application-level output quality. Using these values we train
two different regression models, one for estimating the hardware parameters
and the other for estimating the QoR of the approximate accelerator. The
hardware estimator requires the features of the ACs and the approximate
accelerators as input to predict the area, power, or latency of the variants,
while the software estimator requires the mean and average error of the
approximate circuits deployed in the application to estimate its QoR. Once the
models are trained, they are used to selectively and iteratively explore the
architecture-space of the accelerator using the NSGA-II and ES evolutionary
algorithms, with $\mu$ set to $1$ for the latter. The scripting for generating
the set of approximate accelerator variants, executing the genetic algorithm,
and training the statistical learning models are completed in Python. The
final generation of the accelerators obtained this way are synthesized and
simulated to determine their accurate hardware requirements and QoR, which are
then used to construct the set of Pareto-optimal hardware accelerators for the
target application. An overview of this experimental setup is presented in
Fig. 4.

Figure 4: An Illustration of our Experimental Setup and Tool-flow.

##  IV Results and Discussion

Figure 5: Pearson Correlation Coefficient (PCC) and Exploration Time Analysis
of the Pipelines Discussed in Fig. 3. Figure 6: Analyzing the PCC of Xel-
FPGAs’ $3$ Fastest Pipelines (C, D, and F) - Estimating the Power of the HEVC
Accelerator for Various MCM Modules Deployed in the DCT Block.

We illustrate the efficacy of the Xel-FPGAs framework using the High
Efficiency Video Coding application. HEVC has been shown to be resilient to
errors due to its Discrete Cosine Transform (DCT) function. The DCT block is
composed of $8$ to $14$ adders and multipliers, which vary based on the type
of Multiplierless Multiple Constant Multiplier (MCM) architecture deployed.
These arithmetic units can be approximated to increase the system’s energy-
efficiency and/or performance. More information regarding the HEVC accelerator
block can be found in [27].

###  IV-A Pipeline Analysis

First, we evaluate the effectiveness and time requirements of the pipelines
illustrated in Fig. 3. We do this by estimating the correlation between the
feature-set established by the pipeline with the desired hardware parameter
value (area or power or latency), which is obtained via Vivado Synthesis. We
also evaluate the time required for each pipeline to explore our architecture-
space to obtain a set of Pareto-optimal approximate accelerator variants for
the target application. While there is no dependency on the search algorithm
for pipelines (A), (B), and (C), the evolutionary algorithm associated
selection induces a non-zero evaluation overhead of roughly $30$ms for each
approximate accelerator explored using the ABC tool. Fig. 5 illustrates these
experimental results.

As anticipated, configuration (B) is successful in reducing the exploration
time from thousands of hours (in (A)) to just $41$ hours, which tends to
reduce the Pearson correlation coefficient (PCC) to $0.9$. We can further
reduce this to just $2$ minutes by replacing Vivado with ABC in (B), to obtain
(C), while further dropping the PCC by $2\%$. Although this may not seem like
much, when considering a substantially large architecture-space, such as the
HEVC application, the differences in the model estimates can be quite large,
leading the algorithm to explore non-optimal designs. Therefore, an approach
that retains the PCC closer to $0.9$ while further reducing the exploration
time would be preferable. Pipeline (D), which incorporates features from the
approximate accelerator variants, achieves a PCC more than $0.9$, while
reducing the overhead by nearly $95\%$ for exploring $1,000,000$ approximate
accelerator variants. While configuration (E) achieves a higher PCC, it also
increases the exploration time back to $41$ hours, due to the involvement of
the vivado tool-chain, whereas with pipeline (F), the PCC decreases while
reducing the exploration time. Hence, pipeline (D) serves as the ideal
configuration for Xel-FPGAs.

Besides, the reduction in exploration time, we have also performed experiments
to analyze the effectiveness of the three fastest pipelines in estimating the
power of the HEVC accelerator block when four different types of MCM
components are deployed. The results of these experiments are presented in
Fig. 6. Although pipeline (F) appears to achieve better solutions, on average,
we are interested in the pipeline configuration that obtains the best possible
PCC estimates. For the HEVC use-case, pipeline (D) achieves a trade-off better
than the other pipelines.

###  IV-B Ablation studies

Next, we analyze the capability of three different statistical regression
models, namely,  (1) Random Forest,  (2) Bayesian Regression, and  (3) Support
Vector Machines,  by deploying them in our HEVC use-case. We explore these
models’ effectiveness in accurately estimating the QoR and power of the
approximate accelerator variants using the PCC metric and pipeline (D), as
established earlier. Fig. 7 illustrates the results of these models for the
four different MCM accelerator components deployed in the application. As can
be observed, the different models of regression achieve varying levels of
correlation based on the parameter being estimated. For instance, Random
Forest achieves the best results when estimating the QoR of all the
accelerators in question, i.e., MCM$1$ through MCM$4$. Whereas Bayesian
Regression achieves the best power estimates for all the accelerators.
Therefore, we use these two models, of the nearly $20$ models pre-available in
the Xel-FPGAs framework, to estimate the relevant parameters when exploring
the architecture-space.

Figure 7: Pearson Correlation Coefficient (PCC) Analysis for Three Different
Statistical Regression Models that are Used to Estimate HEVC Accelerators’
Hardware and Quality Parameters.

###  IV-C Architecture Exploration

Figure 8: Analysing the Convergence of ES and NSGA-II. Figure 10: Analysis of
the Pareto-optimal HEVC Approximate Accelerators from Each Generation of ES
(top) and NSGA-II (bottom). Figure 9: Analyzing the Percentage of Total
Pareto-optimal Designs (Density) Generated by ES and NSGA-II in Each
Generation.

Since we have determined the optimal regression models and their corresponding
configurations for the use-case with Xel-FPGAs, we now move on to analyzing
the framework’s model generation and exploration capabilities. Our assumption
is that each new generation of the models explored tends to improve the trade-
offs, with respect to the previous generation, thereby leading to a set of
Pareto-optimal hardware accelerators by the end of our search. We do this by
first analyzing the capabilities of the evolutionary algorithms used in our
framework. Besides the Hill climber algorithm, which was not very beneficial
for our use-case, we have evaluated the NSGA-II and the Evolution Strategy
(ES) algorithms on the HEVC application, as stated in Section II-B. These
algorithms were executed with similar population settings and ran with the
same design time overhead. Fig. 8 illustrates the number of Pareto-optimal
points obtained by the ES and NSGA-II algorithms in each generation of their
execution. Fig. 9, on the other hand, illustrates the percentage of total
number of Pareto-optimal designs obtained, termed density, in each iteration
of algorithm execution. From the former, we can observe that the simpler
MCM$3$ and MCM$4$ modules converge faster, i.e., reach a maximum earlier in
the number of Pareto-optimal designs obtained, as opposed to the more complex
MCM$2$ and MCM$1$ modules. The simpler modules have hardly achieved better
solutions since generation $300$, which is more clearly visible in ES. On the
other hand, NSGA-II requires a few more generations of search before achieving
a similar number of Pareto-optimal designs when compared to ES. The reason
behind this becomes clearer when we analyze Fig. 9. A large portion of the
Pareto-optimal designs obtained by the NSGA-II algorithm are obtained in the
first $300$ iterations, especially for the more complex modules. Hence, a
second peak is observed for the MCM$1$ and MCM$2$ modules where-in another
large percentage of the Pareto-optimal designs are identified. Since the NSGA-
II algorithm starts with a large population set it is able to identify a
significant percentage of the Pareto-optimal designs in the earlier generation
of search. The same is not true for the ES algorithm, which has a more uniform
distribution in identifying the Pareto-optimal designs through the search
iterations. Hence, ES might be more suitable for identifying a large number of
design trade-offs for the target application, whereas NSGA-II might be better
suited to achieve good quality solutions quickly, which contradicts our
initial assumption.

Next, we move on to the main goals of our framework. Fig. 10 illustrates the
result of generating hardware accelerators using the ES and NSGA-II
algorithms, wherein we depict the estimated PSNR and power values of the
approximate accelerator variants searched by each generation. Although the
points identified by Generation $0$ illustrate a random distribution, as
expected, only Generation $1000$’s outcomes are clearly visible while the
others outcomes appear to be missing. This is because the outcomes of
Generation $1000$ are masking the outcomes of the other generations, i.e., the
newer generation outcomes mask the earlier generation outcomes, because their
outcomes are similar to the past generations. This offers the opportunity of
further reducing the exploration time by a factor of $10$, to reduce the
overall exploration time by $99\%$. Moreover, as anticipated, there are a lot
more Pareto-optimal points that are obtained by the ES algorithm when compared
to NSGA-II search.

###  IV-D Final Evaluation

Figure 11: Comparison of the Approximate Accelerators Obtained by our Xel-
FPGAs Framework and the State-of-the-Art [15]. Figure 12: Final Comparison
Between Xel-FPGAs and the SoA [15].

With all the intermediate evaluations complete, we move on to to evaluating
the overall effectiveness of the framework and comparing it with the current
state-of-the-art in approximate accelerator exploration. We have re-
implemented the technique discussed in ApproxFPGAs [15], to generate
approximate accelerator variants with a set of Pareto-optimal FPGA-based ACs
for our HEVC use-case. Fig. 11 illustrates the accelerator-level comparison
between the set of Pareto-optimal approximate accelerator variants obtained by
our Xel-FPGAs framework and the state-of-the-art ApproxFPGAs approach. As can
be observed, for each of the four accelerator types (MCM$1$ through MCM$4$),
our framework outperforms the state-of-the-art and random search to obtain a
better set of Pareto-optimal approximate accelerator variants for the HEVC
application. The use of ES increases the number of Pareto-optimal designs
obtained, i.e., it further pushes the Pareto-frontier by obtaining candidate
designs that lie on the Pareto-frontier. Similarly, the combination of
evolutionary algorithms in the Xel-FPGAs framework clearly identifies better
design trade-offs in comparison to ApproxFPGAs, when we perform an
application-level analysis of the Pareto-optimal approximate accelerators
obtained, as illustrated by the Pareto-fronts in Fig. 12. Moreover, we have
also managed to identify design points that reduce the area requirements of
the applications, in terms of the number of lookup tables required, which is
another benefit of using our Xel-FPGAs framework.

###  IV-E Scalability

Finally, in order to evaluate the ability of our framework to scale for a
multi-stage application, we consider an approximate version of the bio-signal
processing Pan-Tompkins algorithm discussed in [28]. The algorithm is made up
of five application stages each of which is a filter composed using a varying
combination of adders and multipliers, accumulating to $73$ arithmetic units.
Hence, the primary goal of this experiment is to evaluate whether the
framework is better suited for decomposing the exploration into multiple-
stages (n-stage), i.e., a search for each application stage followed by a
combined search with Pareto-optimal designs from each stage, or to design
everything at once ($1$-stage). As emphasized in the article [28], we use the
application-level quality metric –number of heartbeats– to evaluate the
quality of an application stage by embedding the behavioral model of an
approximate variant into the application. Fig. 13 presents an overview of the
results obtained when using the $1$-stage and $2$-stage (in this case,
$n$$=$$2$) exploration strategy. We can observe that the $2$-stage strategy
designs, obtained from both ES and NSGA-II, significantly out-perform the
designs obtained from the $1$-stage. Hence, we can conclude that the
hierarchical and successive approximation strategy ($n$-stage) is better at
identifying an improved set of Pareto-optimal designs, which are more suitable
for the application.

Figure 13: Analyzing the Scalability of Xel-FPGAs.

##  V Conclusion

We presented Xel-FPGAs, which is an end-to-end automated framework for
architecture-space exploration of approximate accelerators given a target
application that is to be deployed in FPGAs. We train and use statistical
regression models to estimate the hardware requirements and output quality of
an application, thereby effectively reducing the time needed for exploring the
architecture-space. We also propose to use the ABC tool as a feature extractor
for the ACs and our approximate accelerators, to further reduce the
exploration time of our framework while not lowering the quality of our
estimations. Using these features we efficiently explore the architecture-
space of approximate accelerators, iteratively using the NSGA-II algorithm, to
identify a set of $200$ Pareto-optimal approximate accelerators at the end of
$1000$ generations. These designs are re-synthesized and simulated to compute
their accurate hardware requirements and quality to determine the final set of
Pareto-optimal approximate accelerator variants for the target application.
The trade-offs between these approximate accelerator variants can be analyzed
by the system designer to select one that meets the system specification and
quality requirements of the application. The Xel-FPGAs framework outperforms
the state-of-the-art in identifying a better set of Pareto-optimal approximate
accelerator variants, while also reducing the exploration time by nearly
$95\%$. We have also illustrated the scalability of our framework using a
multi-stage application to discuss the benefits of using a hierarchical search
strategy. Our Xel-FPGAs framework is open-source and accessible online at
https://github.com/ehw-fit/xel-fpgas.

## Acknowledgment

This work was partially supported by the Doctoral College Resilient Embedded
Systems, which is run jointly by the TU Wien’s Faculty of Informatics and the
UAS Technikum Wien, and partially by the Czech Science Foundation project
21-13001S. This research is also supported by ASPIRE, the technology program
management pillar of Abu Dhabi’s Advanced Technology Research Council (ATRC),
via the ASPIRE Awards for Research Excellence.

## References

  * [1] S. M. S. Trimberger, “Three ages of fpgas: a retrospective on the first thirty years of fpga technology: this paper reflects on how moore’s law has driven the design of fpgas through three epochs: the age of invention, the age of expansion, and the age of accumulation,” _IEEE Solid-State Circuits Magazine_ , vol. 10, no. 2, pp. 16–29, 2018.
  * [2] R. Watanabe, S. Ura, Q. Zhao, and T. Yoshida, “Implementation of fpga building platform as a cloud service,” in _Proceedings of the 10th International Symposium on Highly-Efficient Accelerators and Reconfigurable Technologies_ , 2019, pp. 1–6.
  * [3] L. H. Crockett, R. A. Elliot, M. A. Enderwitz, and R. W. Stewart, _The Zynq book: embedded processing with the ARM Cortex-A9 on the Xilinx Zynq-7000 all programmable SoC_. Strathclyde Academic Media, 2014.
  * [4] P. Mantovani, E. G. Cota, K. Tien, C. Pilato, G. Di Guglielmo, K. Shepard, and L. P. Carloni, “An fpga-based infrastructure for fine-grained dvfs analysis in high-performance embedded systems,” in _Proceedings of the 53rd Annual Design Automation Conference_ , 2016, pp. 1–6.
  * [5] S. Huda, M. Mallick, and J. H. Anderson, “Clock gating architectures for fpga power reduction,” in _2009 International Conference on Field Programmable Logic and Applications_. IEEE, 2009, pp. 112–118.
  * [6] V. K. Chippa, S. T. Chakradhar, K. Roy, and A. Raghunathan, “Analysis and characterization of inherent application resilience for approximate computing,” in _Proceedings of the 50th Annual Design Automation Conference_ , 2013, pp. 1–9.
  * [7] M. Samadi, D. A. Jamshidi, J. Lee, and S. Mahlke, “Paraprox: Pattern-based approximation for data parallel applications,” in _Proceedings of the 19th international conference on Architectural support for programming languages and operating systems_ , 2014, pp. 35–50.
  * [8] W. El-Harouni, S. Rehman, B. S. Prabakaran, A. Kumar, R. Hafiz, and M. Shafique, “Embracing approximate computing for energy-efficient motion estimation in high efficiency video coding,” in _Design, Automation & Test in Europe Conference & Exhibition (DATE), 2017_. IEEE, 2017, pp. 1384–1389.
  * [9] V. K. Chippa, D. Mohapatra, K. Roy, S. T. Chakradhar, and A. Raghunathan, “Scalable effort hardware design,” _IEEE Transactions on Very Large Scale Integration (VLSI) Systems_ , vol. 22, no. 9, pp. 2004–2016, 2014.
  * [10] H. Saadat, H. Javaid, and S. Parameswaran, “Approximate integer and floating-point dividers with near-zero error bias,” in _Proceedings of the 56th Annual Design Automation Conference 2019_ , 2019, pp. 1–6.
  * [11] S. Hashemi, R. I. Bahar, and S. Reda, “Drum: A dynamic range unbiased multiplier for approximate applications,” in _2015 IEEE/ACM International Conference on Computer-Aided Design (ICCAD)_. IEEE, 2015, pp. 418–425.
  * [12] S. Ullah, S. Rehman, B. S. Prabakaran, F. Kriebel, M. A. Hanif, M. Shafique, and A. Kumar, “Area-optimized low-latency approximate multipliers for fpga-based hardware accelerators,” in _Proceedings of the 55th annual design automation conference_ , 2018, pp. 1–6.
  * [13] B. S. Prabakaran, S. Rehman, M. A. Hanif, S. Ullah, G. Mazaheri, A. Kumar, and M. Shafique, “Demas: An efficient design methodology for building approximate adders for fpga-based systems,” in _2018 Design, Automation & Test in Europe Conference & Exhibition (DATE)_. IEEE, 2018, pp. 917–920.
  * [14] S. Ullah, S. S. Sahoo, N. Ahmed, D. Chaudhury, and A. Kumar, “Appaxo: Designing app lication-specific a ppro x imate o perators for fpga-based embedded systems,” _ACM Transactions on Embedded Computing Systems (TECS)_ , vol. 21, no. 3, pp. 1–31, 2022.
  * [15] B. S. Prabakaran, V. Mrazek, Z. Vasicek, L. Sekanina, and M. Shafique, “Approxfpgas: Embracing asic-based approximate arithmetic components for fpga-based systems,” in _2020 57th ACM/IEEE Design Automation Conference (DAC)_. IEEE, 2020, pp. 1–6.
  * [16] J. Echavarria, S. Wildermann, A. Becher, J. Teich, and D. Ziener, “Fau: Fast and error-optimized approximate adder units on lut-based fpgas,” in _2016 International Conference on Field-Programmable Technology (FPT)_. IEEE, 2016, pp. 213–216.
  * [17] S. Ullah, S. S. Sahoo, and A. Kumar, “Clapped: A design framework for implementing cross-layer approximation in fpga-based embedded systems,” in _2021 58th ACM/IEEE Design Automation Conference (DAC)_. IEEE, 2021, pp. 475–480.
  * [18] S. Ullah, S. Rehman, M. Shafique, and A. Kumar, “High-performance accurate and approximate multipliers for fpga-based hardware accelerators,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , vol. 41, no. 2, pp. 211–224, 2021.
  * [19] S. Ullah, S. S. Murthy, and A. Kumar, “Smapproxlib: Library of fpga-based approximate multipliers,” in _Proceedings of the 55th Annual Design Automation Conference_ , 2018, pp. 1–6.
  * [20] Z. Vasicek and L. Sekanina, “Search-based synthesis of approximate circuits implemented into fpgas,” in _2016 26th International Conference on Field Programmable Logic and Applications (FPL)_. IEEE, 2016, pp. 1–4.
  * [21] S. Venkataramani, V. K. Chippa, S. T. Chakradhar, K. Roy, and A. Raghunathan, “Quality programmable vector processors for approximate computing,” in _Proceedings of the 46th Annual IEEE/ACM International Symposium on Microarchitecture_ , 2013, pp. 1–12.
  * [22] V. Mrazek, M. A. Hanif, Z. Vasicek, L. Sekanina, and M. Shafique, “autoax: An automatic design space exploration and circuit building methodology utilizing libraries of approximate components,” in _Proceedings of the 56th Annual Design Automation Conference 2019_ , 2019, pp. 1–6.
  * [23] V. Mrazek, R. Hrbacek, Z. Vasicek, and L. Sekanina, “Evoapprox8b: Library of approximate adders and multipliers for circuit design and benchmarking of approximation methods,” in _Design, Automation & Test in Europe Conference & Exhibition (DATE), 2017_. IEEE, 2017, pp. 258–261.
  * [24] A. Mishchenko _et al._ , “Abc: A system for sequential synthesis and verification,” _URL http://www. eecs. berkeley. edu/alanmi/abc_ , vol. 17, 2007.
  * [25] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: Nsga-ii,” _IEEE transactions on evolutionary computation_ , vol. 6, no. 2, pp. 182–197, 2002.
  * [26] S. H. Jacobson and E. Yücesan, “Analyzing the performance of generalized hill climbing algorithms,” _Journal of Heuristics_ , vol. 10, pp. 387–405, 2004.
  * [27] Z. Vasicek, V. Mrazek, and L. Sekanina, “Towards low power approximate dct architecture for hevc standard,” in _Design, Automation & Test in Europe Conference & Exhibition (DATE), 2017_. IEEE, 2017, pp. 1576–1581.
  * [28] B. S. Prabakaran, S. Rehman, and M. Shafique, “Xbiosip: A methodology for approximate bio-signal processing at the edge,” in _Proceedings of the 56th Annual Design Automation Conference 2019_ , 2019, pp. 1–6.