\subsection{Proof of Lemma \ref{lemma:recoveredsafeset}}\label{appendix:lemma_recoveredsafeset}
\begin{proof}
Assume for a contradiction that we compute $\Tilde{\mathcal{S}}$ as in Equation \eqref{eq:reachrecoveredsafeset} and that
\vspace{-0.25em}
\begin{equation}
\begin{aligned}
\small{
\Tilde{\mathcal{S}} \not\subseteq \text{BRT}^C}
\end{aligned}
\vspace{-0.25em}
\end{equation}
which implies
\vspace{-0.25em}
\begin{equation}
\begin{aligned}
\small{
\exists x^{'}: x^{'} \in \Tilde{\mathcal{S}}, x^{'} \not\in \text{BRT}^C}
\end{aligned}
\vspace{-0.25em}
\end{equation}
$x^{'} \in \Tilde{\mathcal{S}}$ implies that
\vspace{-0.25em}
\begin{equation}
\begin{aligned}
\small{
\Tilde{V}(x^{'}, 0) > \delta}
\end{aligned}
\vspace{-0.25em}
\end{equation}
by Equation \eqref{eq:reachrecoveredsafeset}.

On the other hand, $x^{'} \not\in \text{BRT}^C$ implies that $V(x^{'}, 0) \le 0$.
We know that $\forall (x,t), J_{\Tilde{\policy}}(x,t) \le V(x, t)$ from Equation \eqref{eq:valuefunc}.
Thus, $J_{\Tilde{\policy}}(x^{'},0) \le V(x^{'}, 0) \le 0$.
Then
\vspace{-0.25em}
\begin{equation}
\begin{aligned}
\small{
\Tilde{V}(x^{'}, 0) \le \delta}
\end{aligned}
\vspace{-0.25em}
\end{equation}
by Equation \eqref{eq:reachSafetyMetric}.

We have found the contradiction that $\Tilde{V}(x^{'}, 0) > \delta$ and $\Tilde{V}(x^{'}, 0) \le \delta$.
Thus, our assumption must be false, and it must be true that if we compute $\Tilde{\mathcal{S}}$ as in Equation \eqref{eq:reachrecoveredsafeset}, then $\Tilde{\mathcal{S}} \subseteq \text{BRT}^C$.

To prove the equivalent of Lemma \ref{lemma:recoveredsafeset} but in the case when $\mathcal{L}$ represents a set of desirable states, we use $\text{BRT}$ instead of $\text{BRT}^C$ and flip the value and cost inequalities.
\end{proof}

\subsection{Lemma \ref{lemma:scenarioOptimization}}\label{sec:lemma_scenarioOptimization}
\begin{lemma}
\label{lemma:scenarioOptimization}
Select a violation parameter $\epsilon \in (0, 1)$ and a confidence parameter $\beta \in (0, 1)$. Pick $N$ such that
\vspace{-0.25em}
\begin{equation}
\begin{aligned}
\small{
N \ge \frac{2}{\epsilon}\left(\ln{\frac{1}{\beta}}+1\right)}
\end{aligned}
\vspace{-0.25em}
\end{equation}
Then, with probability at least $1-\beta$, the solution $\hat{\delta}$ to Algorithm \ref{alg:scenarioOptimization} executed with $N$ satisfies the following condition.
\vspace{-0.25em}
\begin{equation} \label{eq:metricbound}
\begin{aligned}
\small{
\underset{\{x: x\in X, \Tilde{V}(x,0) > \hat{\delta}\}}{\mathbb{P}}\left(\left(\Tilde{V}(x,0): J_{\Tilde{\policy}}(x,0) \le 0\right) > \hat{\delta}\right) \le \epsilon}
\end{aligned}
\vspace{-0.25em}
\end{equation}
\end{lemma}
\begin{proof}
Lemma \ref{lemma:scenarioOptimization} follows from Theorem 1 in Scenario Optimization \cite{campi2009scenario}. To use the theorem, we need to prove the following conditions:
\begin{enumerate}
    \item computing $\delta$ can be converted into a standard Scenario Optimization problem
    \item Algorithm \ref{alg:scenarioOptimization} obtains $\hat{\delta}$ by sampling i.i.d from the same space that $\delta$ is optimized over, $X$, but with a probability distribution that is uniform over $\{x: x\in X, \Tilde{V}(x,0) > \hat{\delta}\}$
\end{enumerate}
$\delta$ can be formalized as the following optimization problem:
\begin{align*}
    &\min{g} \\
    &\text{s.t. } \forall x\in X, \left(\Tilde{V}(x,0): J_{\Tilde{\policy}}(x,0) \le 0\right) \le g
\end{align*}
This is a semi-infinite optimization problem where the constraints are linear, and thus convex, in the optimization variable $g$ for any given $x$. Thus, Lemma \ref{lemma:scenarioOptimization} follows from Theorem 1 in Scenario Optimization \cite{campi2009scenario} by replacing $c$ by $1$, $\gamma$ by $g$, $\Delta$ by the state space of interest $X$, and $f$ by $\left(\Tilde{V}(x,0): J_{\Tilde{\policy}}(x,0) \le 0\right) - g$. To use Theorem 1, we also require that i.i.d samples are chosen according to the uniform distribution over $\{x: x\in X, \Tilde{V}(x,t) > \hat{\delta}\}$. This can be proven by observing that, in Algorithm \ref{alg:scenarioOptimization}, the last iteration of the while loop samples $x$ randomly and independently from $\{x: x\in X, \Tilde{V}(x,t) > \hat{\delta}\}$ where $\hat{\delta}$ is exactly the final metric returned. Importantly, the algorithm does not update $\hat{\delta}$ in the final loop (since the loop subsequently breaks).

To prove the equivalent of Lemma \ref{lemma:scenarioOptimization} but in the case when $\mathcal{L}$ represents a set of desirable states, we flip the value and cost inequalities. We also take a maximum instead of a minimum and flip the inequality with $g$ when we formalize $\delta$ as an optimization problem.
\end{proof}

\subsection{Proof of Theorem \ref{theorem:scenarioOptimization}}\label{appendix:theorem_scenopt}
\begin{proof}
We rewrite the LHS of Inequality \eqref{eq:metricbound} of Lemma \ref{lemma:scenarioOptimization} in Appendix\ref{sec:lemma_scenarioOptimization} by substituting in $\hat{\mathcal{S}}$ by Equation \eqref{eq:reachapproximatelyrecoveredsafeset}:
\vspace{-0.25em}
\begin{align*}
\underset{\{x: x\in X, \Tilde{V}(x,0) > \hat{\delta}\}}{\mathbb{P}}&\left(\left(\Tilde{V}(x,0): J_{\Tilde{\policy}}(x,0) \le 0\right) > \hat{\delta}\right) \\
= \underset{\{x: x\in X, \Tilde{V}(x,0) > \hat{\delta}\}}{\mathbb{P}}&\left(\left(\Tilde{V}(x,0) > \hat{\delta}\right) \cap \left(J_{\Tilde{\policy}}(x,0) \le 0\right)\right) \\
= \underset{\{x: x\in X, \Tilde{V}(x,0) > \hat{\delta}\}}{\mathbb{P}}&\left(J_{\Tilde{\policy}}(x,0) \le 0\right) \\
= \underset{x \in \hat{\mathcal{S}}}{\mathbb{P}}&\left(J_{\Tilde{\policy}}(x,0) \le 0\right)
\end{align*}
\vspace{-0.25em}
where the second-to-last line follows because the probability is over the set of states that are already given to have a learned value greater than $\hat{\delta}$.

We know that $\forall (x,t), J_{\Tilde{\policy}}(x,t) \le V(x, t)$ from Equation \eqref{eq:valuefunc}. Thus, it follows that:
\vspace{-0.25em}
\begin{equation}
\begin{aligned}
\small{
\underset{x \in \hat{\mathcal{S}}}{\mathbb{P}}\left(V(x,0) \le 0\right) \le \underset{x \in \hat{\mathcal{S}}}{\mathbb{P}}\left(J_{\Tilde{\policy}}(x,0) \le 0\right) \le \epsilon}
\end{aligned}
\vspace{-0.25em}
\end{equation}

To prove the equivalent of Theorem \ref{theorem:scenarioOptimization} but in the case when $\mathcal{L}$ represents a set of desirable states, we flip the value and cost inequalities.
\end{proof}

\section{Background: Hamilton-Jacobi (HJ) Reachability and DeepReach} \label{sec:hj_reachability}
In this work, we build upon Hamilton-Jacobi reachability analysis to compute an approximation of the safe set. Here, we provide a quick overview of Hamilton-Jacobi reachability analysis and a specific toolbox to compute high-dimensional reachable sets, DeepReach.

\subsection{Hamilton-Jacobi (HJ) Reachability.}
In HJ reachabilty, the computation of $\text{BRT}$ is formulated as an optimal control problem. We will explain it in the context of $\mathcal{L}$ being a set of undesirable states. In the end, we will comment on the case where $\mathcal{L}$ represents a set of desirable states.

To construct this problem, we first define a target function $l(x)$ such that the sub-zero level of $l(x)$ yields $\mathcal{L}$:
\vspace{-0.25em}
\begin{equation}
\small{
    \mathcal{L} = \{x: l(x) \le 0\}}
    \vspace{-0.25em}
\end{equation}
$l(x)$ is commonly a signed distance function to $\mathcal{L}$. For example, we can choose $l(x) = \sqrt{p_x^2+p_y^2}-R$ for our Dubins3D Avoid running example.

Next, we define the cost function of a state corresponding to some policy $u(\cdot)$ to be the minimum of $l(x)$ over its trajectory:
\vspace{-0.75em}
\begin{equation}
    \label{eq:costfunctional}
    \small{
    J_{u(\cdot)}(x,t) = \min_{\tau \in [t, T]} l(\traj_{\state,\tvar}^{\ctrl}(\tau)).}
    \vspace{-0.25em}
\end{equation}
\normalsize
Since the system wants to avoid $\mathcal{L}$, our goal is to maximize $J_{u(\cdot)}(x,t)$.
Thus, the value function corresponding to this optimal control problem is:
\vspace{-0.25em}
\begin{equation}
\small{
    \label{eq:valuefunc}
    V(x,t) = \sup_{u(\cdot)} \{J_{u(\cdot)}(x,t)\}.}
    \vspace{-0.25em}
\end{equation}
\normalsize

By defining our optimal control problem in this way, we can recover $\text{BRT}$ using the value function. In particular, the value function being sub-zero implies that the target function is sub-zero somewhere along the optimal trajectory, or in other words, has reached $\mathcal{L}$. Thus, $\text{BRT}$ is given as the sub-zero level set of the value function at the initial time:

\small
\vspace{-0.8em}
\begin{equation}
    \label{eq:BRT_from_valfunc}
    \text{BRT} = \{x: x\in X, V(x,0) \le 0 \}
\end{equation}
\normalsize
The value function in Equation \eqref{eq:valuefunc} can be computed using dynamic programming, resulting in the following final value Hamilton-Jacobi-Bellman Variational Inequality (HJB-VI):

\small
\vspace{-0.5em}
\begin{equation}
\begin{aligned}
    \label{eq:HJBVI}
    \min\Big\{D_t V(x,t)+ H(x,t), l(x)-V(x,t)\Big\} = 0,
    \end{aligned}
    \vspace{-0.25em}
\end{equation}
\normalsize
with the terminal value function $V(x,T) = l(x)$.
$D_t$ and $\nabla$ represent the time and spatial gradients of the value function.
$H$ is the Hamiltonian that encodes the role of dynamics and the optimal control.

\small
\vspace{-0.5em}
\begin{equation}
    \label{eq:ham}
    \begin{aligned}
    H(x,t) = \max_u & \langle \nabla V(x,t), \dyn(x,u)\rangle.
        \end{aligned}
    \vspace{-0.25em}
\end{equation}
\normalsize

The value function in Equation \eqref{eq:valuefunc} induces the optimal safety controller:

\small
\vspace{-0.25em}
\begin{equation}
    \label{eq:opt_ctrl}
    \begin{aligned}
    u^*(x,t) = \underset{u}{\arg\max} & \langle \nabla V(x,t), \dyn(x,u)\rangle.
        \end{aligned}
    \vspace{-0.25em}
\end{equation}
\normalsize
Intuitively, the safety controller aligns the system dynamics in the direction of the value function's gradient, thus steering the system towards higher-value states.

We have just explained the case where $\mathcal{L}$ represents a set of undesirable states. When the system instead wants to reach $\mathcal{L}$, a infimum is used instead of supremum in Equation $\eqref{eq:valuefunc}$.
The control wants to reach $\mathcal{L}$, hence there is a minimum instead of a maximum in Equation \eqref{eq:ham} and Equation \eqref{eq:opt_ctrl}.

Traditionally, the value function is computed by solving the HJB-VI over a discretized grid in the state space. Unfortunately, doing so involves computation whose memory and time complexity scales exponentially with respect to the system dimensionality, making these methods practically intractable for high-dimensional systems, such as those beyond 5D. Fortunately, a recent deep learning approach, \text{DeepReach}{}, has been proposed to enable Hamilton-Jacobi reachability for high-dimensional systems.

\subsection{\text{DeepReach}{} Approximate Solutions.}
Instead of solving the HJB-VI over a grid, DeepReach represents the value function as a sinusoidal deep neural network (DNN) to learn a parameterized approximation of the value function \cite{bansal2021deepreach}. Thus, memory and complexity requirements for training scale with the value function complexity rather than grid resolution.
To train the DNN, DeepReach uses self-supervision on the HJB-VI itself.
Ultimately, it takes as input a state $x$ and time $t$, and it outputs a learned value $\Tilde{V}(x,t)$.
$\Tilde{V}(x,t)$ also induces a corresponding policy $\Tilde{\pi}(x,t)$.
We refer interested readers to \cite{bansal2021deepreach} for further details.

The learned $\Tilde{V}(x,t)$ from DeepReach can be used to obtain a BRT as in Equation \eqref{eq:BRT_from_valfunc}, but it will only be as accurate as $\Tilde{V}(x,t)$. Unfortunately, like any learning method, $\Tilde{V}(x,t)$ can deviate substantially from the true $V(x,t)$.

Our goal is to recover the biggest safe set $\hat{\mathcal{S}}$ using $\Tilde{V}(x,t)$ that we can (probabilistically) guarantee to be fully contained within the true safe set.
Although we work with \text{DeepReach}{} solutions in particular for this problem setup, our proposed approach in Section \ref{sec:approach} can verify any general $\Tilde{V}(x,t)$ and $\Tilde{\pi}(x,t)$, regardless of whether \text{DeepReach}{}, a level-set method, or some other tool is used to obtain them.

\section{Discussion and Future Work}
In this work, we present an approach to compute an error bound for DeepReach solutions to recover a provably safe approximation of the true reachable tube.
This allows us to utilize the power of learning-based reachability methods for providing safety assurances for high-dimensional dynamical systems.
We apply our method to obtain reachable tubes for high-dimensional rocket-landing and multi-vehicle collision-avoidance problems which traditional methods struggle with.

In the future, we will explore a more refined approach to error correction (as opposed to a uniform correction), as discussed briefly in the end of Section \ref{sec:approach}.
Other directions include additionally considering worst-case disturbances in the system dynamics and using the verification approaches proposed here for a targeted refinement of DeepReach solutions.

\section{Introduction}
Imagine the real-world motivating example of an aerial vehicle that wants to reach a destination while avoiding any collisions. There might even be adversarial disturbances, such as wind or another aircraft nearby, that have the potential to hinder the vehicle from achieving its objective. Since such autonomous systems are becoming increasingly common in our everyday lives, it is vital that we can compute controllers with provable safety guarantees. There already exist many approaches in the literature for computing these. For example, ... However, they are limited by ... The popular Hamilton-Jacobi reachability approach, on the other hand, addresses many of these limitations.

\textbf{Hamilton-Jacobi Reachability Methods.} Hamilton-Jacobi (HJ) reachability analysis is a popular formal verification method for autonomous systems that computes both a safe controller and its corresponding safe state configurations. Its advantages over other approaches includes its ability to handle general nonlinear system dynamics, bounded adversarial system disturbances, and state and input constraints.

In reachability analysis, we are interested in the \textit{Backward Reachable Tube (BRT)} of a dynamical system, which is defined in the following Section \ref{sec:problem_setup} in Equation \eqref{eq:reachBRT}, since it directly relates to system safety.
Hamilton-Jacobi reachability analysis computes a system's BRT by solving a Hamilton-Jacobi (HJ) partial differential equation (PDE) \eqref{eq:HJIVI}, as explained in Section \ref{sec:problem_setup}. Level set methods can compute BRTs accurately up to numerical precision on a grid representing a discretization of the state space \cite{mitchell2005time, lygeros2004reachability, Margellos11, fisac2015reach, bansal2017hamilton}.
Unfortunately, this process results in an exponential scaling of memory and computational complexity with respect to system dimensionality, making it challenging to scale beyond 5D systems \cite{bansal2017hamilton}.

Several methods have been proposed to overcome this challenge by trading off between the class of dynamics they can handle, the approximation quality of the BRT, and the required computation. For example, the methods presented in \cite{Greenstreet1998, Frehse2011, Kurzhanski00, Kurzhanski02, Maidens13, girard2005reachability, althoff2010computing} analyze high-dimensional affine systems using sets of pre-specified shapes, sum-of-squares optimization \cite{Majumdar13, majumdar2017funnel} and parallelotopes \cite{Dreossi16} rely on specific dynamics for computing BRTs for polynomial dynamics \cite{henrion2014convex}; leverage monotonicity of dynamics \cite{coogan2015efficient}; use decoupled or reduced order dynamics to decompose the BRT computation \cite{chen2018decomposition, holmes2020reachable, liu2020leveraging, Chen2016d, somil2020sequential}; leverage the Hopf-Lax formula to compute BRTs for convex dynamics \cite{chow2017algorithm}; or deal with affine dynamics \cite{bak2019numerical, Nilsson2016}.
Control Lyapunov (and Barrier) functions are also used to design safety controller \cite{ames2016control, ames2019control}; however, finding feasible Lyapunov functions can be challenging, especially in the presence of control or disturbance bounds.

Sampling and learning-based methods have also been explored for computing BRTs \cite{liebenwein2018sampling}, including scenario optimization \cite{sutter2017data, devonport2020estimating, ghosh2019new} and adversarial sampling \cite{lew2020sampling}, Support Vector Machines \cite{allen2014machine} and DNNs \cite{rubies2019classification, fisac2019bridging, darbon2020overcoming, niarchos2006neural, djeridane2006neural} leveraging data to compute approximate BRTs for nonlinear systems; however, they often struggle in the presence of adversarial disturbances and state constraints.
Other methods require explicit supervision for learning which is hard to generate for general high-dimensional systems.

Unfortunately, all these methods either need strong assumptions on dynamics, cannot handle presence of exogenous inputs or hard constraints, or lead to a very conservative approximation of the reachable set.
Additionally, very few of these approaches can handle adversarial inputs; yet, considering adversarial inputs (e.g., other agents or the model uncertainty) and state constraints is critical in analysis
and synthesis of safe controllers. Thus, in this work, we focus on a recently developed deep learning-based approach to approximately solve high-dimensional Hamilton-Jacobi reachability problems: \text{DeepReach}{}, which overcomes all of these limitations.

\textbf{DeepReach.} What sets \text{DeepReach}{} apart is its ability to compute BRTs as well as the \textit{corresponding safety controller} for \textit{general} nonlinear dynamical systems in the presence of \textit{disturbances} and \textit{state and input constraints}.
In contrast to level-set methods, DeepReach learns continuous, rather than discretized, solutions to high-dimensional systems which would otherwise be intractable.

Instead of solving this PDE over a grid, \text{DeepReach}{} draws inspiration from the recent progress in solving PDEs using deep learning, and represents the value function as a sinusoidal deep neural network (DNN) to learn a parameterized approximation of the value function via self-supervision on the PDE.
\text{DeepReach}{} can also be interpreted as an approximate dynamic programming method \cite{bertsekas1995dynamic} where a DNN is used to \textit{learn} the basis features for the value function.
Thus, the computation and memory requirements for obtaining the value function do not scale with the grid resolution, but rather the complexity of the value function.
This enables \text{DeepReach}{} to scale to general, high-dimensional dynamical systems.

Even if we can now obtain an approximate BRT for high-dimensional systems for the first time with DeepReach, they are limited in usefulness by their approximate nature. Recall that one of the primary use cases of Hamilton-Jacobi reachability is to provide formal, not approximate, safety guarantees. Thus, our work bridges this gap by proposing a verification method to transform a DeepReach solution into a formal guarantee.

\textbf{Existing Verification Methods.}

\textbf{\text{Scenario Optimization Verification}{} Method.}
\begin{enumerate}
    \item list specific contributions of this work
\end{enumerate}

\section{Approach}\label{sec:approach}
Our approach is to compute a minimal error correction to $\Tilde{V}(x,t)$ such that the corrected value function can be used to extract a (probabilistically) safe set.
We will explain our approach in the context of $\mathcal{L}$ being a set of undesirable states. In the end of each subsection, we will comment on the case where $\mathcal{L}$ represents a set of desirable states. In Section \ref{subsec:errormetric}, we propose an error metric for correcting the value function. In Section \ref{subsec:scenopt}, we propose a practical method to compute this error metric through scenario optimization.

\subsection{Error Metric}\label{subsec:errormetric}
We propose a new error metric $\delta_{\Tilde{\vfunc},\Tilde{\policy}}$ for the value function correction. It is defined as the maximum learned value of an empirically unsafe initial state under the induced policy $\Tilde{\pi}$:
\vspace{-0.25em}
\begin{equation}
\label{eq:reachSafetyMetric}
\begin{aligned}
\small{
\delta_{\Tilde{\vfunc},\Tilde{\policy}} \coloneqq \max_{x\in X}\{\Tilde{V}(x,0): J_{\Tilde{\policy}}(x,0) \le 0\}},
\end{aligned}
\end{equation}
where $J_{\Tilde{\policy}}(x,0)$ is the cost function associated with the trajectory obtained by using the policy $\Tilde{\pi}(x,t)$ from an initial state $x$ and initial time $t=0$ until $T$ (see Equation \eqref{eq:costfunctional}).
Here on, we use $\delta$ as a shorthand for $\delta_{\Tilde{\vfunc},\Tilde{\policy}}$ for brevity purposes.

Intuitively, $\delta$ finds the tightest level of the learned value function that separates the states that are safe under the induced policy $\Tilde{\pi}(x,t)$ from the ones that are not.
Thus, any initial state within the super-$\delta$ level set of $\Tilde{V}(x,0)$ is guaranteed to be safe under the (possibly sub-optimal) policy $\Tilde{\pi}(x,t)$.
Lemma \ref{lemma:recoveredsafeset} formalizes this claim.
\begin{lemma}
\label{lemma:recoveredsafeset}
If we compute $\Tilde{\mathcal{S}}$ as:
\vspace{-0.25em}
\begin{equation}
\label{eq:reachrecoveredsafeset}
\begin{aligned}
\small{
\Tilde{\mathcal{S}} = \{x\in X: \Tilde{V}(x,0) > \delta\}}
\end{aligned}
\vspace{-0.25em}
\end{equation}
then $\Tilde{\mathcal{S}} \subseteq \text{BRT}^C$.
\end{lemma}

The proof for Lemma \ref{lemma:recoveredsafeset} can be found in Appendix\ref{appendix:lemma_recoveredsafeset}.
Intuitively, Lemma \ref{lemma:recoveredsafeset} states that if we could exactly compute $\delta$, then we can recover a $\Tilde{\mathcal{S}}$ that is guaranteed to be a subset of the true safe set.
Furthermore, it is easy to see that, by definition, $\delta$ is the smallest (uniform) value adjustment required on $\Tilde{V}(x,t)$ to guarantee its safety.
Thus, $\Tilde{\mathcal{S}}$ is the largest safe set we can recover from the learned value function by such a uniform error correction procedure.

Unfortunately, computing $\delta$ is a challenging optimization problem, since both the value function and the cost function in the metric definition \eqref{eq:reachSafetyMetric} are typically non-convex functions of $x$.
In the next section, we propose a scenario optimization-based approach to compute a high-confidence bound on $\delta$ and thereby generate a high-confidence safe set.

\begin{remark}
We assume there exists some unsafe states in the system so that $\delta_{\Tilde{\vfunc},\Tilde{\policy}}$ is well-defined. $\delta_{\Tilde{\vfunc},\Tilde{\policy}}$ otherwise can trivially be defined as $-\infty$.
\end{remark}

So far, we have discussed an error correction metric for the case where $\mathcal{L}$ represents a set of undesirable states.
When $\mathcal{L}$ represents a set of desirable states, we use a minimum instead of a maximum and flip the cost inequality in the metric definition \eqref{eq:reachSafetyMetric}.
To recover the safe set, we extract the sub-$\delta$ level set of the value function in Equation \eqref{eq:reachrecoveredsafeset}.

\subsection{\text{Scenario Optimization Verification}{} Method}\label{subsec:scenopt}
We will compute a high-confidence bound on $\delta$ by utilizing a random sampling procedure referred to as \textit{scenario optimization} in the systems and control design literature \cite{campi2009scenario}. The \text{Scenario Optimization Verification}{} Method to compute an approximation $\hat{\delta}$ is summarized in Algorithm \ref{alg:scenarioOptimization}.

\setlength{\textfloatsep}{0pt}
\begin{algorithm}
\caption{\text{Scenario Optimization Verification}{}}
\label{alg:scenarioOptimization}
\begin{algorithmic}[1]
\REQUIRE $X$, $N$, $M$, $\Tilde{V}(x,0)$, $J_{\Tilde{\policy}}(x,0)$
\STATE $\delta_0 \gets -\infty$
\FOR{$i= 0, 1, \ldots, M-1$}
    \STATE $\mathcal{D}_i \gets$ Sample $N$ states IID from $\{x:x\in X, \Tilde{V}(x,0) > \delta_i\}$ \label{line:sample}
    \IF{$\exists x \in \mathcal{D}_i : J_{\Tilde{\policy}}(x,0) \le 0$}
        \STATE $\delta_{i} \gets \max_{x \in \mathcal{D}_i}\{\Tilde{V}(x,0): J_{\Tilde{\policy}}(x,0) \le 0\}$ \label{line:assign}
    \ELSE
        \STATE \textbf{break}
    \ENDIF
\ENDFOR
\RETURN $\hat{\delta} := \delta_i$
\end{algorithmic}
\end{algorithm}

At a high-level, Algorithm \ref{alg:scenarioOptimization} computes a converging sequence of $\delta_i$ that approximates $\delta$ via random sampling until we no longer find any safety violations or reach a maximum number of iterations, $M$.
Specifically, at each iteration $i$, we randomly sample $N$ initial states within the super-$\delta_i$ level set of $\Tilde{V}$ (Line \ref{line:sample}) and compute the costs of associated trajectories $J_{\Tilde{\policy}}(x,0)$ obtained by using $\Tilde{\pi}(x,t)$ as the controller (Equation \eqref{eq:costfunctional}). We next compute the maximum learned value among the states that violate safety and use that as the new estimate of $\delta$ (Line \ref{line:assign})), essentially discarding a level region of $\Tilde{V}(x, 0)$ that were empirically found to be unsafe under $\Tilde{\pi}$.
Thus, with each iteration, we obtain a tighter approximation of $\delta$.
We repeat these steps until no more safety violations are encountered or a maximum number of iterations are achieved.

While the method described in Algorithm \ref{alg:scenarioOptimization} is simple to understand and execute, scenario optimization provides a formal guarantee for bounding $\delta$ with $\hat{\delta}$.
Thus, we can use the approximated $\hat{\delta}$ instead of $\delta$ to compute an approximated safe set $\hat{\mathcal{S}}$ similar to the one defined in Lemma \ref{lemma:recoveredsafeset}:
\vspace{-0.25em}
\begin{equation}
\label{eq:reachapproximatelyrecoveredsafeset}
\begin{aligned}
\small{
\hat{\mathcal{S}} = \{x: x\in X, \Tilde{V}(x,0) > \hat{\delta}\}}
\end{aligned}
\vspace{-0.25em}
\end{equation}
Crucially, we can make a formal probabilistic safety guarantee for $\hat{\mathcal{S}}$. This is summarized in Theorem \ref{theorem:scenarioOptimization} and proven in Appendix\ref{appendix:theorem_scenopt}.
\begin{theorem}[Scenario Optimization Verification Theorem]
\label{theorem:scenarioOptimization}
Select a violation parameter $\epsilon \in (0, 1)$ and a confidence parameter $\beta \in (0, 1)$. Pick $N$ such that
\vspace{-0.25em}
\begin{equation} \label{eq:prescribedN}
\begin{aligned}
\small{
N \ge \frac{2}{\epsilon}\left(\ln{\frac{1}{\beta}}+1\right)}
\end{aligned}
\vspace{-0.25em}
\end{equation}
Then, with probability at least $1-\beta$, the recovered safe set $\hat{\mathcal{S}}$ from Equation \eqref{eq:reachapproximatelyrecoveredsafeset} satisfies the following condition:
\begin{equation}\label{eq:guarantee}
\begin{aligned}
\small{
\underset{x \in \hat{\mathcal{S}}}{\mathbb{P}}(V(x,0) \le 0) \le \epsilon}
\end{aligned}
\vspace{0.50em}
\end{equation}
\end{theorem}

Note that here we state our results for the case where Algorithm \ref{alg:scenarioOptimization} terminates prior to reaching the maximum number of iterations (i.e., the sequence $\delta_i$ converges).
Similar results can be proved for the other case as well (and are not here, for the lack of space).
In practice, we find that Algorithm \ref{alg:scenarioOptimization} converges within 3-4 iterations.

Disregarding the confidence parameter $\beta$ for a moment, the above theorem states that the volume of the unsafe states within $\hat{\mathcal{S}}$ is smaller than or equal to the prescribed $\epsilon$ value, as long as we sample enough states (to satisfy Inequality \eqref{eq:prescribedN}) during the computation of $\hat{\delta}$.
As $\epsilon$ approaches $0$, the number of safety violations in the recovered safe set also approaches $0$. In
turn, the simulation effort grows unbounded since $N$ is inversely
proportional to $\epsilon$.

To interpret the confidence parameter $\beta$, note that $\hat{\mathcal{S}}$ is a random variable that depends on a randomly sampled set of initial states. It may be the case that we just happen to draw a poorly representative sample, in which case the $\epsilon$ bound does not hold. $\beta$ controls the probability of this adverse event happening, which regards the correctness of the final guarantee in Equation \eqref{eq:guarantee}. Fortunately, $N$ only grows logarithmically with $\frac{1}{\beta}$, so $\beta$ can be chosen to be an extremely small value such as $10^{-12}$. $1-\beta$ should then be so close to $1$ that it does not have any practical importance.

When $\mathcal{L}$ represents a set of desirable states, we initialize $\delta_0$ to be infinity instead of negative infinity, flip the value and cost inequalities, and take a minimum instead of maximum in Algorithm \ref{alg:scenarioOptimization}. Also, we flip the value inequalities in Equation \eqref{eq:reachapproximatelyrecoveredsafeset} and Equation \eqref{eq:guarantee}.
\vspace{0.2em}

\noindent \textit{\textbf{Running example: Dubins3D Avoid.}} We now demonstrate the \text{Scenario Optimization Verification}{} method on the Dubins3D Avoid running example introduced in Section \ref{sec:problem_setup}. We choose system parameters $v=0.6 m/s, u_{\min}=-1.1 rad/s, u_{\max}=1.1 rad/s, R=0.25m$, and use DeepReach to learn the system's value function.
The used neural network parameters and architecture are the same as in the DeepReach paper \cite{bansal2021deepreach}, i.e., a DNN with 3 hidden layers and 512 neurons in each layer.

For this first example, we purposefully sabotage the training of DeepReach by permitting it to train for only 8,000 epochs, far less than $100$K epochs used in the DeepReach paper.
We do this to illustrate our method's utility even with low-quality learned solutions.
We apply our method to a well-trained DeepReach solution in Section \ref{sec:dubins3davoid}.

The overall training took less than 30 minutes on an NVIDIA GeForce RTX 3090Ti.
We execute our verification approach choosing $\beta = 10^{-16}, \epsilon = 10^{-3}$, resulting in an $N = 75683$ as per Theorem 1.
Thus, we will be $1-10^{-16}$ confident that at least $1-10^{-3}$ of our recovered set will be provably safe.
Algorithm \ref{alg:scenarioOptimization} converges in 4 iterations in this case, resulting in a $\hat{\delta} = 0.3728$.
Slices of the trained BRT (yellow), the recovered BRT (blue), and the ground truth BRT boundary (black) are shown in Figure \ref{fig:baddubinsavoid}.
The ground truth is computed by a state-of-the-art PDE solver, \textit{Level Set Toolbox (LST)} \cite{LSToolbox}, which computes the value function over a discrete grid of size $101\times 101\times 101$.

\begin{figure}[ht]
    \vspace{-0.65em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_1.png}
    \vspace{-1.5 em}
    \caption{Dubins3D Avoid: Slices of the trained, recovered, and ground truth border BRTs for three values of $\theta$ from a sabotaged DeepReach solution. The recovered BRT completely encompasses the ground truth BRT and provides a provably safe approximation of the safe set.}
    \vspace{-0.65em}
    \label{fig:baddubinsavoid}
\end{figure}

As evident from the figure, the recovered BRT slices completely encompass the true unsafe states (the blue region encompasses the states within the black boundary).
Hence, the complement of the recovered BRT is safe, as we expect based on Theorem \ref{theorem:scenarioOptimization}.
In the rightmost slice, observe that the trained BRT is smaller than the ground truth BRT and is thus unsafe as it is.
Our verification approach correctly expands the trained BRT past the ground truth boundary, ensuring safety.
Validation by sampling 1M states within the recovered safe set reveals a violation rate of $3\times 10^{-6} << \epsilon$.

Note that the trained solution's safety violations are concentrated near $\theta=\pi$.
It is likely these errors which are responsible for the large expansion of the recovered BRT even for slices that are almost accurate (left and middle slices in Figure \ref{fig:baddubinsavoid}).
This is because our method expands the BRT uniformly.
However, too much expansion is undesirable, since we want to recover as much of the safe set as possible.

One way to overcome this challenge is to apply refined error correction to $\Tilde{V}(x,t)$ by independently verifying separate regions of the states-space (i.e., separately consider each region as our $X$).
This allows us to recover sets with the same safety guarantee but expanded differently depending on their own region's error.
However, this method's efficacy depends heavily on the choice of sub-regions, which are hard to prespecify, especially for high-dimensional systems.
Instead, we can resort to a data-driven approach.

Specifically, we trained a simple multilayer perceptron (MLP) on 1M randomly sampled initial states and their empirical safety violation costs (the training took an hour on a standard GPU).
Even though the predictor may be inaccurate (after all, predicting the safe set is the original challenging problem), we hypothesize that it will have learned something useful enough about the general distribution of errors throughout the state space.
We divided the output of the MLP in 10 bins, each corresponding to a different range of empirical cost.
Thus, we expect the last bin to correspond to the states with maximal safety violations (in this case, that will correspond to states with $\theta=\pi$).
We next ran Algorithm \ref{alg:scenarioOptimization} independently for each bin to compute an error correction.
The corresponding results are shown in Figure \ref{fig:binneddubinsavoid}.
As evident, we get a much tighter approximation of the safe set near the slices that have small errors, allowing us to recover a much bigger safe set overall.
In this case, the recovered safe set volume increased by 35\%.

\begin{figure}[ht]
    \vspace{-0.65em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_3.png}
    \vspace{-1.5em}
    \caption{Dubins3D Avoid: Slices of the trained, recovered, and ground truth border BRTs for three values of $\theta$ from a sabotaged DeepReach solution when MLP binning is used to refine recovery. The recovered BRT completely encompasses the ground truth BRT and provides a provably safe approximation of the safe set that is larger than when binning is not used.}
    \vspace{-0.65em}
    \label{fig:binneddubinsavoid}
\end{figure}

As discussed earlier, the performance improvement of this approach depends heavily on the selection of subregions.
We defer this to future work, and for the rest of this paper, we only focus on uniform value function correction.
\section{Problem Setup} \label{sec:problem_setup}
Consider a dynamical system with state $x \in X \subseteq \mathbb{R}^n$, control $u \in \mathcal{U}$, and dynamics $\dot{x} = \dyn(x, u)$ governing how $x$ evolves over time until a final time $T$.
Let $\traj_{\state,\tvar}^{\ctrl}(\tau)$ denote the state achieved at time $\tau \in [t, T]$ by starting at initial state $x$ and time $t$ and applying control $u(\cdot)$ over $[t,\tau]$.
Let $\mathcal{L}$ represent a target set of states that the agent wants to either reach (e.g. goal states) or entirely avoid (e.g. obstacles).
\vspace{0.2em}

\noindent \textit{\textbf{Running example: Dubins3D Avoid.}} As a running example, consider a simple low-dimensional system in the literature known as Dubins3D. It involves a car with position $(p_x, p_y)$, heading $\theta$, velocity $v$, and steering control $u_1 \in [u_{\min}, u_{\max}]$. The car's state $x$ evolves according to
\begin{align*}
    \dot{p_x} = v\cos{\theta}, \quad \dot{p_y} = v\sin{\theta}, \quad
    \dot{\theta} = u_1
\end{align*}
It wants to avoid a circle with radius $R$ centered at the origin. Thus, we define the target set $\mathcal{L}$ of this system to be:
\begin{align*}
    \mathcal{L} = \{x: \sqrt{p_x^2 + p_y^2} \le R\}
\end{align*}
We use Dubins3D as a benchmark because it is tractable for traditional reachability methods to compute a solution for, and we can thus compare our results with a ground truth.
In Section \ref{sec:examples}, we present high-dimensional problems which traditional methods struggle with.

In this setting, we are interested in computing the system's initial-time Backward Reachable Tube within our state space of interest $X$, which we denote as $\text{BRT}$.
We define $\text{BRT}$ as the set of all initial states in $X$ from which the agent will eventually reach $\mathcal{L}$ within the time horizon $[0,T]$, despite best control efforts:
\vspace{-0.25em}
\begin{equation}
\label{eq:avoidBRT}
\begin{aligned}
\small{
\text{BRT} = \{x: x\in X, \forall u(\cdot), \exists \tau \in [0, T], \traj_{\state,0}^{\ctrl}(\tau) \in \mathcal{L}\}}
\end{aligned}
\vspace{-0.25em}
\end{equation}
When $\mathcal{L}$ represents unsafe states for the system, staying outside of $\text{BRT}$ is desirable.
When $\mathcal{L}$ instead represents the states that the agent wants to reach (e.g., a goal set), $\text{BRT}$ is defined as the set of all initial states in $X$ from which the agent, acting optimally, can eventually reach $\mathcal{L}$ within $[0,T]$. Thus, staying within $\text{BRT}$ is desirable.

Our goal in this work is to compute a provable approximation of the safe set. Specifically, we want to compute an approximation $\hat{\mathcal{S}}$ such that $\hat{\mathcal{S}}\subseteq \text{BRT}^C$ ($\hat{\mathcal{S}}\subseteq \text{BRT}$ when $\mathcal{L}$ represents goal states).
We are especially interested in settings where the system is high-dimensional, with which current state-of-the-art methods struggle.

\section{Case Studies} \label{sec:examples}
We will evaluate our approach on various reachability problems.
First, we show results for the avoid version and the reach version of the low-dimensional Dubins3D system for which we have a ground truth BRT available for comparison.
Next, we recover safe sets for a 9D multivehicle collision avoidance problem and for a 6D rocket landing problem, with which traditional methods struggle.
Like our running example, we use $\beta = 10^{-16}, \epsilon = 10^{-3}, N = 75683$.

\subsection{Dubins3D Avoid} \label{sec:dubins3davoid}
In Figure \ref{fig:dubinsavoid}, we show the results of applying our method to a well-trained DeepReach solution for the same Dubins3D Avoid system as in the running example from Section \ref{sec:approach}.
The solution is trained to consider the periodicity of $\theta$ and for 100K epochs (as opposed to 8K for the running example solution).
The proposed algorithm results in a very small error correction of -0.0016. Interestingly, this shrinks the trained BRT to recover a larger safe set than initially learned.
Validation by sampling 1M states within the recovered safe set reveals a violation rate of $3.2\times 10^{-5} << \epsilon$.

The fact that a much smaller error correction is found for the trained solution shown in Figure \ref{fig:dubinsavoid} than the one in Figure \ref{fig:baddubinsavoid} indicates it is of much higher quality.
This illustrates how the error correction metric $\delta$ can also be used to evaluate the quality of different approximate value function solutions.
For example, we can evaluate the relative performance of different DNN hyperparameters, which is especially helpful when the ground truth BRT is not available for comparison.

\begin{figure}[ht]
    \vspace{-0.65em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_2.png}
    \vspace{-1.5em}
    \caption{Dubins3D Avoid: Slices of the trained, recovered, and ground truth border BRTs for three values of $\theta$ from a well-trained DeepReach solution. The trained BRT completely encompasses the recovered BRT, so we plot the recovered BRT border on top for visualization purposes.}
    \vspace{-0.65em}
    \label{fig:dubinsavoid}
\end{figure}

\subsection{Dubins3D Reach} \label{sec:dubins3dreach}
In Figure \ref{fig:dubinsreach}, we show the results of applying our method to a DeepReach solution trained with the same scheme as for the Dubins3D Avoid solution in Section \ref{sec:dubins3davoid}, but for the reach version.
That is, $\mathcal{L}$ is now a set of desirable states, so a safe approximate BRT should be fully contained by the true BRT.
The trained BRT is just barely crossing the ground truth boundary, and the proposed method cuts it down to just behind the boundary, ensuring system safety.
\begin{figure}[ht]
    \vspace{-0.65em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_4.png}
    \vspace{-1.5em}
    \caption{Dubins3D Reach: Slices of the trained, recovered, and ground truth border BRTs for three values of $\theta$ from a well-trained DeepReach solution. The recovered BRT is correctly a subset of the ground truth BRT.}
    \vspace{-0.65em}
    \label{fig:dubinsreach}
\end{figure}

\subsection{Multivehicle Collision Avoidance} \label{sec:multicollision}
We now consider a 9D collision avoidance system involving 3 independent Dubins3D cars.
The $i$th car has position $(p_{xi}, p_{yi})$, heading $\theta_i$, velocity $v$, and steering control $u_i \in [u_{\min}, u_{\max}]$.
The dynamics of vehicle $i$ are given as:
\begin{align*}
    \dot{p}_{xi} = v\cos{\theta_i}, \quad \dot{p}_{yi} = v\sin{\theta_i}, \quad
    \dot{\theta_i} = u_i
\end{align*}
The (undesirable) target set is given by the states where any of the vehicle pairs is in collision:
\begin{equation*}
    \mathcal{L} = \{x: \min\{d(Q_1, Q_2), d(Q_1, Q_3), d(Q_2, Q_3)\} \le R\}
\end{equation*}
where $d(Q_i, Q_j)$ is the distance between cars $i, j$.
We choose $v=0.6, u_{\min}=-1.1, u_{\max}=1.1, R=0.25$ for our case study.
The results of applying our method to a DeepReach solution are shown in Figure \ref{fig:multi}.

\begin{figure}[ht]
    \vspace{-0.65em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_5.png}
    \vspace{-1.5em}
    \caption{Multivehicle Collision Avoidance: Slices of the trained and recovered BRTs for three values of $\theta_1$ from a DeepReach solution.}
    \vspace{-0.65em}
    \label{fig:multi}
\end{figure}

This high-dimensional example is typically difficult to compute with traditional methods, yet here we demonstrate successful recovery of a safe set with formal safety guarantees.
While the recovery preserves much of the safe set that is learned, a significant amount is pruned off.
The percent reduction in safe set size (calculated from 1M samples) is much larger in this example (26\%) than for the Dubins3D solutions in Section \ref{sec:dubins3davoid} (0\%) and Section \ref{sec:dubins3dreach} (10\%), suggesting that the learned value function is more inaccurate for higher dimensional systems.

To validate safety of the recovered BRT, in Figure \ref{fig:multihist}, we plot the minimum pairwise spare distance between vehicles along trajectories
spawning from states sampled within the trained and recovered safe sets.
Note that the trained safe set contains states which result in a collision, as shown by the mass to the left of the dotted line representing the collision distance. After verification, the recovered safe set shows no such mass.
Validation by sampling 1M samples in the recovered safe set reveals a violation rate of $5\times 10^{-6} << \epsilon$.

\begin{figure}[ht]
    \vspace{0.5em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_6.png}
    \vspace{-1.5em}
    \caption{Multivehicle Collision Avoidance: Histogram of the minimum pairwise spare distance between vehicles along trajectories which spawn from states sampled within the trained and recovered safe sets. The collision distance is indicated by the dotted line.}
    \vspace{0.25em}
    \label{fig:multihist}
\end{figure}

\subsection{Rocket Landing} \label{sec:rocketlanding}
Consider a 6D rocket landing system with position $(p_{x}, p_{y})$, heading $\theta$, velocity $(v_{x}, v_{y})$, angular velocity $\omega$, and torque controls $\tau_1, \tau_2\in [-250, 250]$.

\noindent The dynamics are:
\begin{gather*}
    \dot{p_x} = v_x,~\dot{p_y} = v_y, ~\dot{\theta} = \omega,~\dot{\omega} = 0.3\tau_1,\\ \dot{v_x} = \tau_1 \cos{\theta} - \tau_2 \sin{\theta},~
    \dot{v_y} = \tau_1 \sin{\theta} + \tau_2 \cos{\theta} - g,
\end{gather*}
where $g = 9.81$ is acceleration due to gravity.
The target set is the set of states where the rocket reaches a rectangular landing zone of side length 20 centered at the origin:
\begin{align*}
    \mathcal{L} = \{x: |p_x| < 20.0, p_y < 20.0\}
\end{align*}

In Figure \ref{fig:rocket}, we show the results of applying our method to a DeepReach solution.
We also plot the trajectories emanating from 20 randomly sampled initial states within the recovered safe set, demonstrating that they do indeed reach the target set (the green region) by following the induced policy. Validation by sampling 1M states within the recovered safe set reveals a violation rate of $5\times 10^{-6} << \epsilon$.

\begin{figure}[ht]
    \vspace{-0.65em}
    \centering
    \includegraphics[width=1 \columnwidth]{initial_sub/figures/Figure_7.png}
    \vspace{-1.5em}
    \caption{Rocket Landing: Slices of the trained and recovered BRTs for three values of $v_y$ from a DeepReach solution. The system trajectories safely reach the target set (green) from the initial states within the recovered BRT.}
    \vspace{-0.65em}
    \label{fig:rocket}
\end{figure}

\section{Introduction}
As autonomous systems become increasingly common in our everyday lives, it is vital that we can design provably safe controllers for these systems.
Hamilton-Jacobi (HJ) reachability analysis provides a powerful framework to design such controllers for general nonlinear dynamical systems \cite{lygeros2004reachability, mitchell2005time}.
In reachability analysis, the safe states for the system are characterized through the \textit{Backward Reachable Tube (BRT)} of the system.
This is the set of states that spawn trajectories that will eventually reach some given target set despite the best control effort.
Thus, if the target set represents the set of undesirable states, the BRT represents unsafe states for the system and should be avoided.
Along with the BRT, reachability analysis also provides a safety controller to keep the system outside the BRT.

Traditionally, the BRT computation in HJ reachability is formulated as an optimal control problem.
The BRT can then be obtained as a sub-zero level solution of the corresponding value function.
Obtaining the value function requires solving a partial differential equation (PDE) over a state-space grid, resulting in an exponentially scaling computation complexity with the number of states \cite{bansal2017hamilton}.
To overcome this challenge, a variety of solutions have been proposed that trade off between the class of dynamics
they can handle, the approximation quality of the BRT,
and the required computation.
These include specialized methods for linear and affine dynamics \cite{Greenstreet1998, Frehse2011, Kurzhanski00, Kurzhanski02, Maidens13, girard2005reachability, althoff2010computing, bak2019numerical, Nilsson2016}, polynominal dynamics \cite{Majumdar13, majumdar2017funnel, Dreossi16, henrion2014convex}, monotonic dynamics \cite{coogan2015efficient}, and convex dynamics \cite{chow2017algorithm}
(see \cite{bansal2017hamilton, bansal2021deepreach} for a survey).
Owing to the success of deep learning, there has also been a surge of interest in approximating high-dimensional BRTs \cite{rubies2019classification, fisac2019bridging, djeridane2006neural, niarchos2006neural, darbon2020overcoming} and optimal controllers \cite{onken2022neural} through deep neural networks (DNN).
Building upon this line of work, the authors in \cite{bansal2021deepreach} have proposed DeepReach -- a toolbox that leverages recent advances in neural implicit representations and neural PDE solvers to compute a value function and a safety controller for high-dimensional systems.
Compared to the aforementioned methods, DeepReach can handle general nonlinear dynamics, the presence of exogenous disturbances, as well as state and input constraints during the BRT computation.
However, even if we can now obtain an approximate BRT for high-dimensional systems, it is limited in usefulness by its approximate nature.
In particular, the learned BRT can be overly optimistic; thus, no formal safety guarantees can be provided on the obtained BRT and controller.

In this work, our goal is to compute a BRT and a safety controller for high-dimensional systems with provable safety guarantees.
Specifically, we build upon DeepReach and propose a verification method to provide safety assurances on the DeepReach solution.
The key insight of our work is to rely on the consistency between the learned value function and the implicit safety controller induced by this value function to compute a \textit{uniform correction bound} for the value function.
Essentially, if a particular state is outside the BRT (i.e., ``safe'' as per the learned value function), then starting from this state, the corresponding controller should keep the system trajectory outside the BRT at all times (i.e., the state should also be ``safe'' under the prescribed controller).
Otherwise, we cannot ensure safety of this state and the value function should be corrected such that this state is inside the BRT.
Once all such states have been added to the learned BRT, we have obtained a provably safe approximation of the true BRT.

However, in general, it is challenging to tractably compute the required value function correction across the entire state space, since the learned value function can be highly nonlinear and hard to optimize.
We propose a scenario optimization-based method for computing this correction.
Scenario optimization is a sampling-based method to solve semi-infinite optimization problems, and has been widely used for system and control design \cite{campi2011sampling, calafiore2006scenario, campi2009scenario}.
The proposed method has general applicability and is not restricted to a specific class of
system dynamics or value functions.
Given that the value function correction is obtained via randomization and, hence,
is a random quantity, we provide probabilistic guarantees on the safety of the recovered BRT.
However, this confidence is a design parameter and can be chosen as close to 1 as desired (within a simulation budget).

To summarize, the key contributions of this paper are:
\begin{itemize}
    \item an error correction mechanism for DeepReach that allows us to compute a provable BRT and safety controller for general dynamical systems;
    \item a practical method to compute this error correction that is not restricted to a specific class of systems; and
    \item a demonstration of the proposed approach for various dynamical systems, inspired by rocket landing and multi-vehicle collision avoidance problems.
\end{itemize}

\section*{Appendix}
\input{initial_sub/sections/appendix}

\end{document}


\section{Introduction}
Steganography and watermarking techniques aim to embed digital information into visual media in an invisible manner.
A watermarking scheme typically consists of two components---an \textit{encoder} that embeds a given digital message in the image pixels and a \textit{decoder} that extracts the message from a watermarked image.
Such technology has several applications, such as transmitting secret messages, copyright protection, establishing media ownership, and embedding invisible QR codes into images~\cite{tancik2020stegastamp}.
Another important use case of digital watermarking lies in media authentication or integrity verification by embedding semi-fragile watermarks into digital media at the source.
Semi-fragile watermarks require the property of being fragile to malicious manipulations or tampering while being robust to benign image-processing operations such as image compression, scaling, and color adjustments.

\begin{figure}[h]
\centering
\includegraphics[width=\linewidth]{figs/Final_Full_ICCAD_FastStamp.pdf}
\caption{Schematic diagram of FastStamp watermarking pipeline.
The pipeline is divided into two steps: watermark insertion using a DNN encoder on FPGA (top) and watermark extraction using a DNN decoder on a cloud server (bottom).}
\label{fig:overview}
\end{figure}
Conventional image watermarking systems use hand engineered pipelines~\cite{lin2000detection,ho2004semi,yang2009semi,li2015semi} to embed information in the spatial or frequency domain of an image.
However, the major limitations of traditional approaches lie in the higher visibility of the embedded watermarks, limited message capacity, and low robustness to digital compression techniques like JPEG transforms.
As a result, these traditional approaches are being quickly replaced with deep learning based techniques, which are achieving state-of-the-art results in image watermarking and steganography tasks.

Deep learning based watermarking techniques relies on encoder and decoder convolutional neural networks (CNNs).
These networks are trained end-to-end for the task of embedding and retrieving a given message in an image.
While deep learning techniques significantly outperform hand engineered watermarking pipelines, the improvement comes at the cost of increased computational overhead and memory requirement of these models.
The best performing neural image watermarking encoders are parameterized by around half a million floating point parameters, which makes it challenging to deploy such systems on resource constrained hardware such as FPGAs or handheld devices.
Such techniques have only been implemented at a software level which results in significant latency between image capture and transmission.
Our work aims to embed watermarks in images and videos in real-time as they are being captured.
Embedding the watermarks at the hardware level can not only reduce the latency of the watermarking process but also enable media authentication and provenance by leveraging unique hardware signatures from PUFs~\cite{gu2017improved} or secure enclaves as the watermarking data.

In this work, we propose FastStamp---a light-weight yet robust neural image watermarking framework to enable real-time watermarking on hardware platforms.
Keeping resource constraints in mind, we develop a parameter efficient CNN based watermarking model that can match and even outperform the success metrics of state-of-the-art neural image watermarking models.
Our watermarking model leverages efficient neural blocks such as depthwise separable convolutions, spatial upsampling operations, and linear layers to reduce the memory requirement and inference latency without compromising the watermark retrieval accuracy, message capacity, and imperceptibility.
We then design and verify an FPGA implementation of our watermarking encoder model to allow image watermarking directly on hardware.
Our most optimized design achieves real-time image watermarking and only requires $3$ milliseconds to watermark a given image while achieving the same results as the software implementations.
To the best of our knowledge, we propose the first framework to accelerate deep neural networks for steganography and watermarking of images using FPGAs.
FastStamp can be customized for semi-fragile and robust image watermarking.
While semi-fragile watermarking is useful for checking media integrity, robust image watermarking is useful for establishing media ownership and copyright protection.

\noindent \textbf{Summary of Contributions:}
\begin{itemize}
    \item We develop a parameter efficient and computationally inexpensive watermarking model that can embed recoverable watermarks at a significantly lower cost as compared to prior neural watermarking techniques.
    \item We develop the first FPGA accelerator platform for DNN based image watermarking and steganography.
    We design reusable and reconfigurable basic blocks pertinent to such models as separable convolutions, 2D upsampling, and skip connections.
    We deploy the FastStamp encoder model on Xilinx XCVU13P FPGA, which achieves $68 \times$ faster inference than prior neural watermarking models.

    \item Our framework is end-to-end and supports two types of watermarking schemes: \textit{robust} and \textit{semi-fragile}. FastStamp learns to be robust to a wide range of real-world digital image processing operations such as lighting, color adjustments, and compression techniques.
    Our semi-fragile watermarking scheme learns to be robust to above benign transformations while being fragile to media forgery and local tampering.
\end{itemize}

\vspace{-2mm}
\section{Background}
\subsection{Digital Watermarking}
Digital watermarking techniques broadly seek to generate three different types of watermarks: fragile~\cite{di2019fragile}, robust~,\cite{cox1997secure,shehab2018secure,zhu2018hidden} and semi-fragile~\cite{lin2000detection,yu2017review,ho2004semi}.
Fragile and semi-fragile watermarks are primarily used to certify the integrity and authenticity of image data.
Fragile watermarks are used to achieve accurate authentication of digital media, where even a one-bit change to an image will lead it to fail the certification system.
In contrast, \textit{robust} watermarks aim to be recoverable under several image manipulations to allow media producers to assert ownership over their content even if the video is redistributed and modified.
Semi-fragile watermarks combine the advantages of both robust and fragile watermarks and are mainly used for fuzzy authentication of digital images and identification of image tampering~\cite{yu2017review}.

Prior work has proposed hand-engineered pipelines to embed semi-fragile watermark information in the spatial\cite{xiao2008semi} and frequency domain of images and videos.
Particularly, in the frequency domain, the watermark can be embedded by modifying the coefficients produced with transformations such as the discrete cosine transform (DCT)~\cite{preda2015watermarking,ho2004semi} and discrete wavelet transform (DWT)~\cite{tay2002image,li2015semi,benrhouma2015tamper}.
For example, during 2D DWT an image is first decomposed into various frequency channels using a Haar filter.
A scaled image is used as the watermark and inserted into mid frequency wavelet channel.
Taking 2D inverse DWT of the altered wavelet decomposition produces the watermark embedded image.
The major drawbacks of traditional approaches lie in higher visibility of the embedded watermarks, increased distortions in generated images, and low robustness to compression techniques like JPEG transforms as discussed in Section~\ref{sec:evaluation_metrics}.
Specifically, since digital images shared over social media are highly compressed and undergo various lighting and color adjustments, the watermarks generated by DWT and DCT algorithms break under benign real-world transformations and compression.

More recently, CNNs have been used to provide an end-to-end solution to the watermarking problem. They replace hand-crafted hiding procedures with neural network encoding~\cite{baluja2017hiding,hayes2017generating,zhu2018hidden,zhang2019invisible,wang2021faketagger,tancik2020stegastamp,neekhara2022facesigns,luo2020distortion}.
These techniques train end-to-end encoder CNNs to embed and decode watermarks, which have resulted in lower imperceptibility and more robust recovery of the watermark data.
However, these models are memory intensive and much slower as compared to DCT based algorithms.
In our work, we address the performance limitations of neural watermarking systems and propose a light-weight framework for both robust and semi-fragile image watermarking suitable for hardware acceleration platforms.

\subsection{FPGA Accelerated Techniques}

There have been several efforts to accelerate neural networks on FPGAs~\cite{zhang2015optimizing,suda2016throughput,samragh2017customizing,fastwavepaper}. Equipped with the necessary hardware for basic DNN operations, FPGAs are able to achieve high parallelism and utilize the properties of neural network computation to remove unnecessary logic.
Some prior efforts have been made in FPGA acceleration of convolutional autoencoder architectures~\cite{liu2018optimizing,zhao2018fpga,govorkova2022autoencoders}.
While these works implement many sub-blocks used in neural network computation, we find that they cannot be directly used for an image watermarking framework because efficient implementations of sub-blocks like 2D upsampling and separable convolutions with skip-connections are absent. Moreover, existing neural architectures for image watermarking are not optimized for hardware-software co-design.

Prior work~\cite{kiran2013design,hazra2018fpga,hajjaji2019fpga,khoshki2014fpga} has made significant efforts in accelerating traditional image watermarking schemes that hide secret information in the frequency domain of images and rely on DCT and DWT based algorithms.
As discussed earlier, while such algorithms are vastly popular for hardware applications due to its simplicity and low computational overhead, the resulting watermarked image is often not robust to real-world image transformations. In our work, we develop the first FPGA accelerator platform for DNN based image watermarking and steganography that enables robustness to real-world digital image processing and compression while being selectively fragile to media tampering techniques.

\subsection{Countering Media Forgery}
With the widespread development of deep learning based image and video synthesis techniques~\cite{liu2019stgan,choi2020starganv2,dfsurvey,faceswap}, it has become increasingly easier and faster to generate high-quality convincing fake images/videos such as Deepfakes.
Such manipulated media can fuel misinformation, defame individuals and reduce trust in media.
While considerable research effort has been made in designing CNN based deepfake detectors~\cite{dolhansky2020deepfake,afchar2018mesonet}, these techniques have been shown to hold major security vulnerabilities and can be bypassed by attackers~\cite{advdeepfakes}.
To counter such threats, authors~\cite{9689555,neekhara2022facesigns,wang2021faketagger} have proposed semi-fragile watermarking as a solution to perform media authentication and distinguish deepfake media from real media by verifying a secret watermark embedded in the media.
For both fragile and semi-fragile techniques, a watermark must be inserted when the image is captured, which makes these techniques dependent on both algorithmic and hardware implementation.
If watermark information is embedded separately in images and videos after it is captured by a device, this method may fail in situations where tampering is carried out before inserting the signature or watermark.
While the proposed solution to media authentication is to add a verifiable digital signature or watermark to an image/video using neural networks, prior works~\cite{neekhara2022facesigns,9689555,wang2021faketagger} do not actively implement the technology in resource constrained settings or camera hardware.
Upon empirical study, we found that it is challenging to fit such off-the-shelf models~\cite{tancik2020stegastamp,zhu2018hidden,neekhara2022facesigns} on FPGAs since they were developed without any attention to hardware-software co-design practices and range from 500\,k to 2 million parameters for the encoder model.
To this end, we design our own DNN based watermarking system that utilizes depthwise separable convolutions and a parameter efficient message upsampler to reduce computational overhead while preserving required bit recovery accuracy, capacity, and imperceptibility.

\begin{figure*}
    \centering
    \includegraphics[width=1.0\linewidth]{figs/FastStamp_Training.pdf}
    \caption{FastStamp Encoder-Decoder Training framework.
    In both robust and semi-fragile schemes, the encoder model learns to embed a message in a given image by encouraging retrieval from the decoder model under benign transforms.
    For our semi-fragile design, we additionally encourage fragility to malicious transforms by maximizing message retrieval error on tampered watermark images.}
    \label{fig:model_diagram}
\end{figure*}

\section{Methodology}
\subsection{Training Framework}
\label{sec:training}
Our goal is to develop a learnable, end-to-end model for image steganography and watermarking such that the encoder model can embed a message as a visually invisible perturbation in the image, and the decoder network can extract the message from the watermarked image.
We develop two variants of our training framework to generate \textit{robust} and \textit{semi-fragile} watermarks.
A robust watermark is designed to be recoverable when real-world image transformations are applied.
A semi-fragile watermark is designed to be robust to benign image transformations such as compression, minor color, and contrast adjustments but it should be unrecoverable when malicious image transforms such as image tampering and face-swapping are applied.
A semi-fragile watermark is designed to be robust to benign image transformations such as compression and minor color, and contrast adjustments but it should be unrecoverable when malicious image transforms such as image tampering and face-swapping are applied.

The encoder network $E$ takes as input an image $x$ and a bit string $s \in \{0, 1\}^L$ of length $L$, and produces an encoded (watermarked) image $x_w$.
That is, $x_w=E(x, s)$.
Depending on our task of either \textit{semi-fragile} or \textit{robust} watermarking, the encoded image goes through the following operations:
\begin{enumerate}
\item \textit{\textbf{Robust watermarking:}}
In this setting, the image goes through a benign image transformation $g_b \sim G_b$ to produce $x_b=g_b(x_w)$.
The benign image is then fed to the decoder network, which predicts the message $s_b=D(x_b)$.
For optimizing secret retrieval during training, we use the $L_1$ distortion between the predicted and ground-truth bit strings.
The decoder is encouraged to be robust to benign transformations by minimizing the message distortion $L_1(s, s_b)$.
Therefore the secret retrieval error for an image $L_M(x)$ is obtained as follows:
    \begin{equation}
    L_M(x) = L_1(s, s_b)
    \end{equation}
\item \textit{\textbf{Semi-fragile watermarking:}} In this setting, the watermarked image goes through two image transformation functions---one sampled from a set of benign transformations ($g_b \sim G_b$) and the other sampled from a set of malicious transformations  ($g_m \sim G_m$) to produce a benign image $x_b=g_b(x_w)$ and a malicious image $x_m=g_m(x_w)$.
The benign and malicious watermarked images are then fed to the decoder network, which predicts the messages $s_b=D(x_b)$ and $s_m=D(x_m)$ respectively.
The decoder is encouraged to be robust to benign transformations by minimizing the message distortion $L_1(s, s_b)$; and fragile for malicious manipulations by maximizing the error $L_1(s, s_m)$.
Therefore the secret retrieval error for an image $L_M(x)$ is obtained as follows:
    \begin{equation}
    L_M(x) = L_1(s, s_b) - L_1(s, s_m)
    \end{equation}
\end{enumerate}

The watermarked image is encouraged to look visually similar to the original image by optimizing three image distortion metrics: $L_1$, $L_2$ and $L_\text{pips}$~\cite{zhang2018unreasonable} distortions.
\begin{equation}
\begin{split}
& L_\text{img}(x, x_w) = L_1(x, x_w) + L_2(x, x_w) + c_p L_\text{pips}(x, x_w)
\end{split}
\end{equation}

Therefore, the parameters $\alpha,\beta$ of the encoder and decoder network are trained using mini-batch gradient descent to optimize the following loss over a distribution of input messages and images:
\begin{equation}
\mathbb{E}_{x, s, g_b, g_m} [ L_\text{img}(x, x_w) + c_M L_M(x)  ]
\end{equation}

In the above equations, $c_p$, and $c_M$ are scalar coefficients for the respective loss terms.
We refer the readers to our supplementary material for the values we use for these coefficients and other implementation details.

\subsection{Message encoding} The input of our encoder network is a bit string $s$ of length $L$. This watermarking data includes a secret message or a hardware signature generated by trusted execution environments or PUFs.
To further ensure message secrecy, we can encrypt the message using a stream cipher with a secret key that is shared between the encoding and decoding devices.
In our work we embed messages of length 128 bits in an image of size $128{\times}128$, which allows embedding $2^{128}$ unique messages in each $128{\times}128$ image patch.

\subsection{Model Architecture and Optimization}
\label{sec:architecture}
The encoder model takes as input an image $x$ and a message bit-string $s$ to produce a watermarked image $x_w$.
The encoder model of a typical neural watermarking system follows a convolutional U-Net architecture comprising several downsampling and upsampling layers with skip-connections.
In prior work~\cite{tancik2020stegastamp,neekhara2022facesigns}, the secret message bit-string is first projected using a learnable linear layer reshaped as a matrix to have the same height and width as the input image; and then attached as the fourth channel of the input image.
The combined input and secret image then undergo the downsampling and upsampling operations of the U-Net to produce the watermarked image.

The above described encoder model architecture in prior work~\cite{tancik2020stegastamp,zhu2018hidden,neekhara2022facesigns} has a large memory footprint and is unsuitable for deployment in resource-constrained settings. To reduce the model size without compromising on the watermarking performance, we propose the following architectural optimizations:

\subsubsection{\textbf{Secret Message Upsampler}}
The secret message upsampler projects the message string $s$ to a matrix that gets attached as the fourth channel of the input image.
A na\"{i}ve implementation using a linear layer can result in a large memory footprint since the number of parameters is given by $hwL$, where $h$ and $w$ are the input image height and width, respectively, and $L$ is the secret message length.
An input image size of $128{\times}128$ and a message length of $128$, would result in more than two million parameters.
To optimize the number of parameters, we perform the secret upsampling operation as follows:
\begin{enumerate}
    \item Project the message $s$ to a vector $s_\text{proj}$ of size $h'w'$ using a linear layer.
    \item Reshape $s_\text{proj}$ to a matrix $s_\text{proj}{}_M$ of dimensions $(h', w')$
    \item Upsample $s_\text{proj}{}_M$ to a matrix of size $(h, w)$ using nearest-neighbour upsampling.
\end{enumerate}

\begin{figure}[htp]
\centering
\includegraphics[width=1.0\linewidth]{figs/Upsampling_Secret_Faststamp.pdf}
\caption{An example of optimized secret message upsampling using linear layer projection followed by nearest neighbor 2D upsampling.}
\label{fig:overview}
\end{figure}

The nearest-neighbour upsampling operation is parameter-free and computationally more efficient as compared to matrix-vector multiplication.
Through our experiments, we find that using an $(h', w')$ that are much smaller than $(h, w)$ can achieve the same watermarking while being significantly efficient in both time and memory.
For an image of size $(128, 128)$ and message length $L=128$, we use $h'=w'=16$ thus requiring only $32768$ parameters.
The upsampled secret gets attached as the fourth channel of the input image and undergoes the U-Net downsampling and upsampling operations described below.

\subsubsection{\textbf{U-Net Downsampling}}

The downsampling network in U-Net architecture typically comprises $5$ to $8$ convolutional blocks.
Each block contains a strided convolutional layer, a batch normalization layer, and a non-linear activation like ReLU or leaky ReLU.
The number of output channels of each convolutional layer increase with the depth of the network, doubling at each step until a maximum value is reached.
To optimize this architecture, we first replace the convolutional layers with depthwise and point-wise separable convolutional layers~\cite{depthwise}.
Not only does this optimization reduce the number of parameters but also reduces the number of floating point operations required in each layer computation.
Next, we optimize a number of output channels of each convolutional block.
In our experiments, we perform a design-space exploration to find that we can substantially reduce the number of output channels in each layer without compromising on secret retrieval accuracy and imperceptibility.
Our most optimized design uses $5$ downsampling layers with $64$ output channels in the final separable convolutional layer.
Figure~\ref{fig:unet} details the network architecture and output tensor sizes after each downsampling step.

\begin{figure}[htp]
\centering
\includegraphics[width=\linewidth]{figs/unet.pdf}
\caption{FastStamp Encoder Architecture.
Our encoder network takes as input an image $x$ and the output of the secret message upsampler $s_\text{proj}{}_M$ and generates the watermarked image}
\label{fig:unet}
\end{figure}

\subsubsection{\textbf{U-Net Upsampling}}
The upsampling network in U-Net architecture follows a mirror image of the downsampling network.
Instead of a regular convolutional layer, each upsampling block typically contains a transposed convolutional layer. However, transposed convolutional layers have been shown to introduce unwanted visual artifacts in the generated images~\cite{odena2016deconvolution} and we see the same effect empirically in our work.
To remove such artifacts, we replace the transposed convolution layer with a separable convolution layer followed up by nearest neighbor 2D upsampling following the recommendations given by past work~\cite{odena2016deconvolution}.
At each upsampling step, the output of the corresponding downsampling step is concatenated with the block input to provide skip-connections which are known to improve the performance of encoder-decoder models.
The output of the last upsampling layer undergoes a $\tanh$ activation function to normalize output values between $-1$ and $1$, which are then scaled between $0$ and $1$ to produce an RGB image.
Figure~\ref{fig:unet} details our encoder network architecture.

FastStamp's Decoder follows a similar architecture as the encoder but contains 8 downsampling and 8 upsampling layers. After the upsampling U-net, the decoder network follows the inverse architecture of the secret message upsampler to predict the 128-bit message.

\section{Accelerator Design}

\subsection{Design Overview}

Figure~\ref{fig:overview} gives a high-level overview of our FPGA accelerator design.
Each neural network layer is treated as a separate dataflow stage.
To be low latency and high throughput, all weights and biases are stored on-chip.
Complex activation functions are implemented via precomputed lookup tables.
The design uses task-level pipelining (i.e., HLS dataflow) for each layer and streams the data between each dataflow stage using first-in-first-out buffers (FIFOs).
As FIFOs can only be read once, to implement the skip connections, an additional dataflow stage is used to clone the skip connection data from its input FIFO into two other FIFOs so that it can be read twice for its two datapaths.

\begin{figure}[htp]
\centering
\includegraphics[width=\linewidth]{figs/Faststamp_accelerator_overview.pdf}
\vspace{-6mm}
\caption{Design overview of FastStamp Accelerator Platform.}
\label{fig:overview}
\end{figure}

\vspace{-4mm}
\subsection{Implementation Details}
\label{sec:implementation}
To implement FastStamp on FPGA, we began with using the standard \textit{hls4ml}~\cite{Duarte:2018ite,Aarrestad:2021zos,Fahim:2021cic} framework and added additional libraries to this to support the necessary functionalities such as depthwise separable convolutions, and nearest-neighbor 2D upsampling, concatenation layers (i.e., skip connections) which were previously not supported.
We create custom sub-blocks and modify existing implementations with optimizations that better support our model.
Our design is composed of the following $6$ main architectural sub-blocks:

\subsubsection{\textbf{Linear Layer}}
The linear layer in the \textit{Secret Message Upsampler} module is implemented using efficient matrix-vector multiplication based on product and sum tree.
In order to optimize the design and make efficient use of the DSP blocks, we use a parallelized approach to convert layer computations into multiple MAC operations.
Multiple rows of the weight matrix are processed simultaneously by dividing it into chunks.
In each round, a chunk of the weights matrix is copied to one of the weight buffers while the other weight buffer is fed into the \textit{dot product} modules together with a copy of the input vector.
The iterations end when all rows of the weight matrix have been processed.
Then each \textit{dot-product} function partitions its input vectors into chunks and concurrently executes MAC operations over the partitioned subsets.
The accumulated results of the subsets are then added together within the \textit{reduce\_sum} function to compute the final output.
The \textit{reduce\_sum} module performs a tree-based reduction algorithm which further reduces latency.

\subsubsection{\textbf{Nearest-neighbor 2D Upsampling}}
2D upsampling layers are used in the secret \textit{Secret Message Upsampler} and the \textit{U-Net Upsampling layers}.
To upsample a matrix, we iterate through every element in the matrix and start building the new resized matrix.
We use nearest neighbour as our interpolation method, as it is the least hardware intensive computation available for resizing purposes.
We do not utilize pipelining in this sub-block, as it would cause a drastic increase in resource utilization for this resizing task while not having a substantial effect on latency.
Due to these problems, pipelining makes this sub-block unscalable. Instead, loop unrolling is utilized to reduce loop iterations, giving us latency reductions with negligible impact on resource utilization.

\subsubsection{\textbf{Skip Connections}}
Concatenate layers are used in FastStamp to implement skip connections between upsampling and downsampling layers.
These layers contributed towards some of the highest resource utilization in our design.
To realize FastStamp on hardware, efficient support for concatenation of different sized inputs was implemented.
We also distribute resource utilization to BRAM within these layers due to a large allocation of LUTs.
Alongside this, we are able to minimize latency by pipelining the operations with complete unrolling in these layers.

\subsubsection{\textbf{Separable Convolution}}
Due to resource limitations, traditional convolution techniques in machine learning settings are often not feasible on hardware.
Recent works~\cite{bai2018sepconv, yoo2018cnn} have shown that utilizing separable convolutions significantly reduces the number of multiplications needed.
Instead of traditional convolution,  separable convolution layers are used in both U-Net downsampling and upsampling layers.
A separable convolution comprises of a depthwise separable convolution and point-wise convolution operation.
The depthwise separable convolution is performed independently across the input channels and results in the same number of output channels.
A point-wise convolution operation reduces to matrix-vector multiplication along the channel axis for each spatial cell of the input.
We also use a streaming implementation of depthwise convolution utilizing pipelining with complete loop unrolling, array partitioning, and loop flattening to achieve minimal latency.
We optimize standard point-wise convolution in a similar manner as our linear layer.

\subsubsection{\textbf{Batch-normalization}}
Batch-normalization~\cite{pmlr-v37-ioffe15} layers are used in U-Net downsampling and upsampling layers.
These layers are used after the separable convolution layers to stabilize the network.
Rather than using this optimization, pipelining, loop unrolling, and array partitioning are enforced to accelerate the normalization layer.
This allows for the batch-normalization output to be used as an input to skip connections, a feature that cannot be done when fusing layers.
Our optimizations achieve low latency and minimal resource utilization.

\subsubsection{\textbf{Non-linear activations}}
ReLU and $\tanh$ are the two non-linear activation functions are used in FastStamp.
The ReLU function, which is used in the U-net downsampling and upsampling layers computes the following function for each input $x$, $y = \max(0, x)$.
Due to the simplicity of the ReLU activation function, we use loop unrolling to optimize latency for this operation.
The $\tanh$ activation function computes the following function for each input $x$, $y=(e^x - e^{-x})/(e^x + e^{-x})$.
We use the default $\tanh$ implementation in \textit{hls4ml} which is performed using a pre-computed lookup table to reduce latency.

\section{Experiments and Results}

\setlength\tabcolsep{4pt}
\begin{table*}[htp]
\centering
\begin{tabular}{@{}l|r|rrr|cc|ccc|c@{}}
\multicolumn{1}{r}{} & \multicolumn{1}{|r|}{\emph{Model Size}} & \multicolumn{3}{c|}{\emph{Capacity}} & \multicolumn{2}{c|}{\emph{Imperceptibility}} &
\multicolumn{3}{c|}{\emph{BRA (\%) -- Benign}} &
\multicolumn{1}{c}{\emph{BRA (\%) -- Tampering}}
\\
\midrule
Method & \# Params & $H$,$W$ & $L$ & $\text{BPP}$ & PSNR & SSIM & None & JPG-75 & Filtering & FaceSwap \\
\midrule

DCT (Semi-Fragile)~\cite{ho2004semi} & --- & $128$ & $256$ & $5.2\times10^{-3}$ & $22.49$ & $0.871$ & $99.81$ & $56.65$ & $94.62$ & $85.51$ \\
HiDDeN (Robust)~\cite{zhu2018hidden} & 411\,k & $128$ & $30$ & $6.1\times10^{-4}$ & $27.57$ & $0.934$ & $97.06$ & $72.71$ & $94.52$ & --- \\
StegaStamp (Robust)~\cite{tancik2020stegastamp} & 528\,k & $400$ & $100$ & $2.0\times10^{-4}$  & $29.39$ & $0.925$ & $99.92$ & $99.91$ & $99.84$ & --- \\
\midrule
FastStamp (Robust) & 45\,k & $128$ & $128$ & $2.3\times10^{-3}$ & $30.65$ & $0.942$ & $100.00$ & $99.84$ & $99.78$ & --- \\
FastStamp (Semi-Fragile) & 45\,k & $128$ & $128$ & $2.3\times10^{-3}$ & $30.64$ & $0.940$ & $100.00$ & $99.74$ & $99.72$ & $51.11$ \\

\bottomrule
\multicolumn{1}{c}{}
\end{tabular}
\vspace{-2mm}
\caption{Capacity, imperceptibility, and BRA metrics of different watermarking systems. $H,W$ indicates the height and width of the input image. For both robust and semi-fragile watermarking systems a high BRA is desirable for benign transforms. For semi-fragile watermarking systems, a low BRA is desirable for tampering transforms.}

\label{tab:basicompare}
\end{table*}

\subsection{Dataset}
We conduct experiments on the CelebA dataset~\cite{liu2015faceattributes} which is a large database of over 200,000 face images of 10,000 unique celebrities.
We set aside 1000 images for testing the watermarking models and split the remaining data into 80\% training and 20\% validation.
All FastStamp models are trained using images of size $128{\times}128$, which are obtained after center-cropping and resizing the CelebA images.
We conduct experiments with message bit length $L=128$.

\subsection{Evaluation Metrics}
\label{sec:evaluation_metrics}
For the evaluation of our watermarking techniques, we investigate the following metrics based on prior works~\cite{hore2010image,hajjaji2019fpga}.

\noindent \textbf{1. Imperceptibility:} We compute \textbf{peak signal to noise ratio (PSNR)} and \textbf{structural similarity index (SSIM)} between the watermarked and original images.
A higher value of both these metrics indicates a more imperceptible watermark.

\noindent \textbf{2. Capacity:} Capacity measures the amount of information that can be embedded in the image.
We use \textbf{bits per pixel (BPP)} which is calculated as $L/(HWC)$ where $L$ is the message length and $H,W,C$ indicate the height, width, and channels of the image.
Higher BPP values indicate higher capacity.

\noindent \textbf{3. Bit-Recovery Accuracy
(BRA):} BRA calculates the recovery accuracy of the bit string $s$.
For robustness, we aim to have a high BRA when benign or intended transformations such as JPEG compression or color and contrast adjustments are applied.
For semi-fragile watermarking systems, the goal is to have a low BRA when image tampering operations such as local tampering or FaceSwap are applied while maintaining robustness against benign transformations.

We compare our watermarking framework against three prior works on image watermarking---a DCT based semi-fragile watermarking system~\cite{yang2009semi} and two robust neural image watermarking systems HiDDeN~\cite{zhu2018hidden} and StegaStamp~\cite{tancik2020stegastamp}.
To evaluate the robustness and fragility of our watermarking systems, we perform the following transformations that are unseen during training:
\begin{enumerate}
\item
\textbf{JPEG Compression: }
Digital images are usually stored in a lossy format such as JPEG.
We compress the watermarked images using JPEG-75 compression and measure the decoding BRA.
\item
\textbf{Filtering:}
We apply a set of real-world image filtering operations using the Pilgram library~\cite{pilgramfilter} that simulates photo editing filters that are common on social media. These include color, contrast, and lighting adjustments.
\item
\textbf{Face Swapping: }
For evaluating semi-fragile watermarking systems, we simulate image tampering by performing face swapping using the open source implementation of FaceSwap~\cite{faceswap}.
A low BRA is desirable against this transform to detect tampering.
\end{enumerate}

\vspace{-2mm}
\subsection{Training and Architecture Optimization}
Our encoder model follows the depthwise separable convolutional U-Net architecture as discussed in Section~\ref{sec:architecture}.
To find the optimum architecture, we create different-sized versions of the baseline U-Net architecture by reducing the number of channels in each layer by a factor of 2, 4 and 8.
We find that reducing the number of channels by a factor of 4 does not compromise model performance while being significantly lighter than the base U-Net model.

We train two variants of this optimized design in the robust and semi-fragile settings using the training technique described in Section~\ref{sec:training}.
In the robust setting, we simulate differentiable JPEG compression, Gaussian blur, color, and contrast adjustment as the benign transforms $g_b$ during training.
In the semi-fragile setting, in addition to the benign transforms, we simulate differentiable localized tampering as the malicious transform $g_m$ during training.
We train our models for 200\,k mini-batch iteration with a fixed learning rate of $1.5\times 10^{-4}$ using Adam optimizer.
Table~\ref{tab:basicompare} compares our optimized models \textit{ FastStamp (Robust) } and \textit{ FastStamp (Semi-Fragile) } against prior neural and DCT based image watermarking frameworks.
As compared to prior neural image watermarking and steganography models, FastStamp is significantly smaller and achieves a similar BRA with slightly improved imperceptibility as compared to StegaStamp~\cite{tancik2020stegastamp} and HiDDeN~\cite{zhu2018hidden}.
The improvement in imperceptible metrics is achieved by using nearest neighbor 2D upsampling in the U-Net architecture as opposed to transposed convolutions in prior work.

\setlength\tabcolsep{3pt}
\begin{table*}[htp]
\centering
\begin{tabular}{@{}l|rrrr|ccc|rrr@{}}
\multicolumn{1}{r}{} &
\multicolumn{4}{c|}{\emph{Resource Utilization (\%)}} &
\multicolumn{3}{c|}{\emph{Performance}} & \multicolumn{3}{c}{\emph{Correctness}}
\\
\midrule
  Design & BRAM & FF & LUT & DSP & Clock&  Latency & Throughput & BRA & PSNR & SSIM \\
   \textit{\small{Available Resources}} & 94 Mb & 3456K & 1728K & 12288 & Period (ns)&   (\# Cycles) & (Hz) &  &  &  \\
  \midrule
FixedPoint-32 & $>$100\% & 51\% & $>$100\% & $>$100\% & --- & --- & --- & 100.0 & 30.67 & 0.942\\
FixedPoint-16 & 89\%  & 18\% & $>$100\% & 54\% & --- & --- & ---  & 100.0 & 30.64 & 0.941\\
\textbf{FixedPoint-16-Optimized} & 59\% & 14\% & 72\% & 53\% & 5 & 596823 & 335 & 100.0 & 30.64 & 0.941\\
\bottomrule
\multicolumn{1}{c}{}
\end{tabular}
\caption{Design-space exploration for FPGA implementation of FastStamp on Xilinx XCVU13P FPGA board.
We report the resource utilization, timing performance, and implementation correctness for various designs.
Our optimized 16-bit fixed point implementations fit within the available resources while maintaining the same correctness metrics as the 32-bit implementation.}
\label{tab:fpga}
\end{table*}

\vspace{-2mm}
\subsection{Design Space Exploration}
We implement the individual submodules described in Section~\ref{sec:implementation} using Vivado HLS for the Xilinx XCVU13P FPGA board.
First, we perform a search over the bit-width of the fixed point representation of the network weights and intermediate outputs.
Our goal is to find the lowest bit-width that does not compromise on message recovery and imperceptibility.
Figure~\ref{fig:bitwidth} indicates the BRA and PSNR of different bit-width implementations of FastStamp.
Based on this analysis, we use a 16-bit representation with 6 bits for the integer and 10 bits for the decimal representation.
\begin{figure}[htp]
\centering
\includegraphics[width=\linewidth]{figs/bitwidthupdated.pdf}
\vspace{-6mm}
\caption{Watermarking success metrics for different fixed-point representations.
A high value for both BRA and PSNR is desirable for accurate message recovery and imperceptibility. }
\label{fig:bitwidth}
\end{figure}
Next, we explore pipelining and loop unrolling options in our Vivado HLS implementation of various submodules.
Complete loop unrolling in submodule implementations was an infeasible design choice for FastStamp since it exceeded the available resources on the device.
We found that partial loop unrolling with factors of 8 and pipelining loops that were completely unrolled
were the most effective optimizations for FastStamp.

While effective pipelining and loop unrolling resulted in a significant reduction in resource utilization, the LUT requirement of our design still exceeded the available resources on the device.
Reducing the loop unrolling factor was an effective strategy to reduce LUT utilization but resulted in significant increases in encoding latency.
To avoid latency increase, the next strategy we applied was distributing variables, such as model weights, that were initially all implemented as LUTs, to the BRAM on our board.
We also force operations, such as multiplication in the separable convolution into DSP blocks, which results in lower LUT utilization.
We utilize per layer reuse factor to tune the inference latency versus utilization of FPGA resources and enable parallelization.
This allows us to process multiple MAC operations at every unit of time.
Using all of these optimizations, we are able to store our model entirely in the on-chip memory of the FPGA and perform low latency encoding while avoiding communication with off-chip memory.
Table~\ref{tab:fpga} lists the resource utilization, performance, and correctness of some of the design choices that led to our most optimized design, \textit{FixedPoint-16-Optimized}.
Figure~\ref{fig:sampleimage} shows sample outputs of this optimized FPGA implementation and compares them with the GPU implementation in PyTorch.

\begin{figure}[htp]
\centering
\includegraphics[width=\linewidth]{figs/FastStamp_Sample_images.pdf}
\caption{Sample image outputs of FastStamp optimized design and PyTorch implementation with the original image}
\label{fig:sampleimage}
\end{figure}

\vspace{-3mm}
\subsection{\textbf{Performance and Power Analysis}}\label{sec:power}
Table~\ref{tab:cpugpu} compares the inference time and power requirement of our optimized FPGA implementation with the highly optimized CPU and GPU implementation of FastStamp and the open source implementations of prior neural watermarking systems.
We benchmark the optimized PyTorch implementation of FastStamp on the Nvidia Tesla V100 GPU.
The CPU implementation is a NumPy inference program written by us and optimized fully.
We measure the power consumption for the GPU benchmarks using the Nvidia power measurement tool (\textit{Nvidia-smi}) running on \textit{Linux} operating system, which is invoked during program execution.
For our FPGA implementations, we synthesize our designs using Xilinx Vivado 2020.1.
We then integrate the synthesized modules accompanied by the corresponding peripherals into a system-level schematic using the Vivado IP Integrator.
The frequency is set to 200\,MHz, and power consumption is estimated using the synthesis tool.
Our FPGA implementation achieves ${\sim}68{\times}$ faster speed against prior work's GPU implementation and ${\sim}10{\times}$ faster speed against FastStamp's GPU implementation at a $3{\times}$ lower power requirement.

\setlength\tabcolsep{7pt}
\begin{table}[htp]
\centering
\begin{tabular}{l|rr}
\multicolumn{1}{c|}{Implementation} &
\multicolumn{1}{|r}{Time (ms)} &
\multicolumn{1}{r}{Power (W)} \\
\midrule
StegaStamp GPU~\cite{tancik2020stegastamp} & 205 & 76 \\
HiDDeN GPU~\cite{zhu2018hidden} & 234 & 65 \\
\midrule
FastStamp CPU  & 326 & --- \\
FastStamp GPU  & 30 & 59 \\
FastStamp FPGA & 3  &  19\\
\bottomrule
\multicolumn{1}{c}{}
\end{tabular}
\caption{Power consumption and wall-clock time (in milliseconds) required to generate a single watermarked image per implementation.
}
\label{tab:cpugpu}
\end{table}

\vspace{-8mm}
\section{Conclusion}
We propose an efficient image watermarking model that matches or even outperforms prior neural image watermarking and steganography models while utilizing significantly fewer parameters.
By leveraging an efficient secret message upsampling module, depthwise separable convolutions, and 2-D upsampling, we were able to train a smaller model while preserving success metrics for watermarking tasks.
Finally, we implement our encoder model on an FPGA to achieve $68{\times}$ higher throughput as compared to prior GPU implementations.
Our implementation allows watermark embedding directly at the hardware source, which not only secures the image capture and transmission pipeline but also reduces latency in embedding the watermark.
In the process of this implementation, we develop reconfigurable sub-modules which can accelerate convolutional downsampling and upsampling networks on hardware.

\section{Acknowledgments}
This work was supported in part by National Science Foundation (NSF) Trust-Hub under award number CNS2016737, NSF TILOS under award number CCF-2112665, NSF Cooperative Agreement OAC-2117997 (\href{https://a3d3.ai}{a3d3.ai}), DoD UCR W911NF2020267 (MCA S-001364), and ARO MURI under award number W911NF-20-S-0009.
Thanks to the Fast ML collective (\href{https://fastmachinelearning.org}{fastmachinelearning.org}) for support.

\bibliographystyle{ACM-Reference-Format}