\subsection{Contributions}
Under the described setup, we present the following contributions.
\begin{enumerate}
\item Using $OOG$ as impact metric, and the maximum expected loss as a risk metric, we formulate the risk assessment problem. We observe that the risk assessment problem corresponds to an untractable infinite non-convex robust optimization problem which is NP-hard in general.
\item We propose a convex SDP which solves the risk assessment problem approximately by sampling the uncertainty set. Additionally, we provide  probabilistic guarantees on the original robust optimization problem.
\item We derive the necessary and sufficient conditions for the risk metric to be bounded.
\end{enumerate}

The remainder of the article is organized as follows. The uncertain system and the adversary is described in Section \ref{description}. The problem is formulated in Section \ref{PF}, to which an approximate solution is presented in Section \ref{rational}. The efficacy of the proposed optimization problem is illustrated through a  numerical example in Section \ref{Example}. We conclude the article in Section \ref{Conclusion}.

\textit{Notation}: Throughout this article, $\mathbb{R}, \mathbb{C}, \mathbb{Z}$ and $\mathbb{Z}^{+}$ represent the set of real numbers, complex numbers, integers and non-negative integers respectively. A positive (semi-)definite matrix $A$ is denoted by $A \succ 0 (A \succeq 0)$. Let $x: \mathbb{Z} \to \mathbb{R}^n$ be a discrete-time signal with $x[k]$ as the value of the signal $x$ at the time step $k$. Let the time horizon be $[0,N]=\{ k \in \mathbb{Z}^+|\; 0 \leq k\leq N \}$. The $\ell_2$-norm of $x$ over the horizon $[0,N]$ is represented as $|| x ||_{\ell_2, [0,N]}^2 \triangleq \sum_{k=0}^{N} x[k]^Tx[k]$. Let the space of square summable signals be defined as $\ell_2 \triangleq \{ x: \mathbb{Z}^+ \to \mathbb{R}^n |\; ||x||^2_{\ell_2, [0,\infty]} < \infty\}$ and the extended signal space be defined as $\ell_{2e} \triangleq \{ x: \mathbb{Z}^+ \to \mathbb{R}^n | \;||x||^2_{[0,N]} < \infty, \forall N \in \mathbb{Z}^+ \}$. For the sake of simplicity, we represent $||x||^2_{\ell_2,[0,\infty]}$ as $||x||^2_{\ell_2}$. For $x \in \mathbb{R}, \left \lceil {x} \right \rceil$ represents $x$ rounded to the nearest integer greater than or equal to $x$. Let $(\Omega, \mathcal{D}_a, \mathbf{P})$ represent a probability space with sample space $\Omega \subset \mathbb{R}^v$, $\sigma-$algebra $\mathcal{D}_a$, and probability measure $\mathbf{P}$. Let $(\Omega^w, \mathcal{D}_a^w, \mathbf{P}^w)$ represent the $w$-times Cartesian product of $\Omega$ with the $\sigma-$algebra $\mathcal{D}_a^w$ and the probability measure $\mathbf{P}^w = \mathbf{P} \times \dots \times \mathbf{P}$. A point drawn from $(\Omega^w, \mathcal{D}_a^w, \mathbf{P}^w)$ is thus $(\delta_1,\delta_2,\dots\delta_w)$, i.e., $w$ independent elements in $\mathbb{R}^v$ drawn independently from $\Omega$ according to the same probability $\mathbf{P}$.\footnote{In this article, the Cartesian product is considered over the same probability space ($\Delta$). But this can be generalized to arbitrary probability spaces.}. The relative importance of various uncertainty outcomes of an arbitrary function $f(\delta),\; \delta \triangleq \begin{bmatrix} \delta_1 & \dots & \delta_m \end{bmatrix}^T \in \Omega,$  is represented by $\mathbb{P}_{\Omega}(f(\delta))$. An empty set is denoted by $\emptyset$. For any $a, b \in \mathbb{R}$ and $a \geq b$, ${a \choose b} = \frac{a!}{b!(a-b)!}$. The non-zero elements of a vector $s \in \mathbb{R}^{n}$ is denoted by $\text{supp}(s)$. The cardinality of the non-zero elements in a vector $s \in \mathbb{R}^{n}$ is denoted by  $||s||_0 \triangleq |\text{supp}(s)|$, where $|\text{supp}(s)|$ represents the cardinality of the set $\text{supp}(s)$.
\subsection{Data injection attack scenario}\label{Attack scenario:sec}
In the system described in \eqref{P}-\eqref{D}, we consider an adversary injecting false data into the sensors or actuators of the plant but not both. Next, we discuss the resources the adversary has access to.
\subsubsection{Disruption and disclosure resources}\label{disclosure:sec}
The adversary can access (eavesdrop) the control or sensor channels and can inject data. This is represented by
\[ \begin{bmatrix}
\tilde{u}[k]\\
\tilde{y}[k]
\end{bmatrix}
=
\begin{bmatrix}
{u}[k]\\
{y}[k]
\end{bmatrix}
+
\begin{bmatrix}
B_a\\
D_a
\end{bmatrix}
a[k], \left[\begin{array}{c|c}
B_a^T & D_a^T
\end{array}\right] \triangleq \left[\begin{array}{c|c}
E_a^T & 0\\
0^T & F_a^T
\end{array}\right]\]
where $a[k] \in \mathbb{R}^{n_a}$ is the attack signal injected by the adversary. The matrix $E_a (F_a)$ is a diagonal matrix with $E_a(i,i)=1 \;(F_a(i,i)=1)$, if the actuator (sensor) channel $i$ is under attack and $0$ otherwise.

\subsubsection{System knowledge}
In general, the system operator knows about the parameters of the controller and detector as s/he is the one who design it. However, it is not known to the adversary. To this end, we consider a realistic adversary whose system knowledge is bounded. We defined this adversary as bounded-rational. Henceforth, we refer to the bounded-rational adversary as a rational adversary.
\begin{definition}[Rational adversary]\label{RA}
An adversary is defined to be rational if it knows all the matrices of \eqref{P}-\eqref{D} but with bounded uncertainty. $\hfill\triangleleft$
\end{definition}

That is, in contrary to the operator, a rational adversary has uncertainty also in the system matrices of the detector and the controller ($A,B,C,C_J,D_J,A_c,B_c,C_c,D_c,A_e,B_e,K_e,C_e,D_e,$ and $E_e$). Defining ${x}[k] \triangleq [ x_p[k]^T \; z[k]^T \; s[k]^T]^T$, the closed-loop system under attack with the performance output and detection output as system outputs becomes
\begin{equation}\label{system:CL:uncertain}

                \begin{array}{ll}
                            {x}[k+1] &= {A}_{cl}^{\Delta}{x}[k] + {B}_{cl}^{\Delta}a[k]\\
                            y_p[k] &= {C}_p^{\Delta}{x}[k] + {D}_p^{\Delta} a[k]\\
                            y_r[k] &= {C}_r^{\Delta}{x}[k] + {D}_r^{\Delta} a[k],\\
                \end{array}

\end{equation}
where the nominal matrices are given by $\left[ \begin{array}{c|c}
A_{cl} & B_{cl}
\end{array}\right] \triangleq $
\[\left[
\begin{array}{c|c}
\begin{matrix}
    A+BD_cC & BC_c & 0\\
    B_cC & A_c & 0\\
    (B_eD_c +K_e)C & B_eC_c & A_e
    \end{matrix} & \begin{matrix}
    BB_a + BD_cD_a \\ B_cD_a \\ (B_eD_c+K_e)D_a
    \end{matrix}
\end{array}\right]
\]
\begin{align*}
    {C}_p &\triangleq \begin{bmatrix}
    C_J+D_JD_cC & D_JC_c & 0
    \end{bmatrix},\quad {D}_p \triangleq D_J(D_cD_a+B_a),\\
    {C}_r &\triangleq \begin{bmatrix}
    (D_eD_c + E_e)C & D_eC_c & C_e
    \end{bmatrix}, {D}_r \triangleq (D_eD_c +E_e)D_a.
\end{align*}

Although the uncertainty is present in the adversarial system knowledge, the bounds of the uncertainty is chosen by the defender based on the his/her belief of the adversaries resources.  Next, we assume the following for clarity.
\begin{assumption}\label{assume_stable}
The closed-loop system \eqref{system:CL:uncertain} is stable $ \forall \;\delta \in \Omega$.$\hfill\triangleleft$
\end{assumption}
\begin{assumption}\label{Ba}
The input matrix has full column rank i.e., $ \nexists \;s \in \mathbb{R}^{n_a} \neq 0$ such that $B_{cl}^\Delta s =0$.$\hfill\triangleleft$
\end{assumption}
\subsubsection{Attack goals and constraints.}
Given the resources the adversary has access to, the adversary aims at disrupting the system's behavior while staying stealthy. The system disruption is evaluated by the increase in energy of the performance output whereas, the adversary is stealthy if the energy of the detection output is below a predefined threshold $\epsilon_r=1$.
\begin{definition}[Stealthy attack]
An attack vector $a$ is said to be stealthy if it satisfies $||y_r||_{\ell_2}^2 \leq 1$.$\hfill\triangleleft$
\end{definition}
We discuss the attack policy for a deterministic system next.

\subsection{Optimal attack policy for the nominal system}

From the previous discussions, it can be understood that the goal of the adversary is to maximize the performance cost while staying undetected. When the system \eqref{system:CL:uncertain} is deterministic, i.e., $\Omega=\emptyset$, \cite{teixeira2015strategic} formulates the attack policy of the adversary as the following non-convex optimization problem
\begin{equation}\label{opti:o2o:uncertain}
\begin{aligned}
||\Sigma||_{\ell_{2e},y_p \leftarrow y_r}^2 \triangleq \sup_{a \in \ell_{2e}} \; & ||y_p||^2_{\ell_{2}} \\
\textrm{s.t.} \quad & ||y_r||_{\ell_2}^2 \leq 1,\; x[0]=0,\;x[\infty]=0,
\end{aligned}
\end{equation}
where $||\Sigma||_{\ell_{2e},y_p \leftarrow y_r}^2$ represents the $OOG$ that characterizes the disruption caused by the attack signal $a$. In \eqref{opti:o2o:uncertain}, the constraint $x[0]=0$ is introduced because the system is assumed to be at equilibrium before the attack.

\begin{assumption}\label{equil_begin}
The closed-loop system \eqref{system:CL:uncertain} is at equilibrium $x[0]=0$ before the attack commences.$\hfill\triangleleft$
\end{assumption}

We assume that the adversary has finite energy (similar to $H_{\infty}/H_{\_}$ approaches) and thus attacks the system for a sufficiently long but finite amount of time, say $T$. Although this time $T$ is unknown \textit{a priori}, it can be ensured that $ a[\infty] \triangleq \lim_{k\to \infty} a[k] = 0$ by setting $x[\infty] \triangleq \lim_{k\to \infty} x[k] = 0$.
\begin{lemma}\label{vanish_lemma}
If $x[\infty]=0$, then it holds that $a[\infty]=0$.
\end{lemma}
\begin{proof}
If $x[\infty]=0$, it follows that
\begin{equation}
    \lim_{k\to \infty} \Vert x[k+1] - A_{cl}x[k]\Vert = 0 =     \lim_{k\to \infty} \Vert B_{cl}a[k] \Vert = \Vert B_{cl}a[\infty]\Vert.
\end{equation}
From \textit{Assumption \ref{Ba}}, we known that $B_{cl}a[\infty] = 0 $ holds only when $a[\infty]=0$. This concludes the proof.
\end{proof}

For the above discussed reason, the constraint $x[\infty]=0$ is introduced in the optimization problem \eqref{opti:o2o:uncertain}. In this article, we focus on attacks which satisfy the constraint $x[\infty]=0$, and define them as state-vanishing attacks.
\begin{definition}[State-vanishing attack]\label{vanish_attack_def}
An attack vector $a$ is defined to be state-vanishing if applying $a$ to \eqref{system:CL:uncertain} generates the state vector $x$ which satisfies $\lim_{k \to \infty}x[k] = 0$.$\hfill\triangleleft$
\end{definition}

In the literature, such characterization of impact of stealthy attacks \eqref{opti:o2o:uncertain} has only been studied for deterministic systems. The remainder of this article is focused on discussing methods to quantify the impact in terms of risk on the uncertain system \eqref{system:CL:uncertain}.

The concept of risk is conventionally adopted to address decision-making in an uncertain environment \cite[Chapter 14]{ferrari2021safety}. Since we also focus on decision-making (impact assessment) under uncertainty, it is useful to adopt tools from the risk management framework. Thus, before introducing the problem formulation, a brief introduction to risk management and risk metrics is provided as it helps in the problem formulation.
\subsection{Omniscient adversary}\label{problem1}

With the random variable defined in \textit{Definition \ref{RV}}, we next formulate the risk assessment problem in for a rational adversary.  Consider the data injection attack scenario where only the bounds of the parametric uncertainty set $\Omega$ is known to the adversary. Although, the adversary can determine the attack vector which maximizes the expected loss over the entire uncertainty set $\Omega$ i.e., the adversary maximizes the function
\begin{equation}
\begin{aligned}
 ||\bar{\Sigma}(\delta)||_{\ell_{2e},y_p \leftarrow y_r}^2 \triangleq &\; \mathbb{E}_{\Omega}\left\{X^A(a,\delta)\right\}.
\end{aligned}
\end{equation}

This setup is common in game-theoretic approaches \cite{lye2005game} where the players do not know the strategy of the other players and thus play by maximizing its expected return over all the strategies of the other players. Similarly, since the adversary has limited system knowledge, s/he chooses an attack policy which maximizes the expected loss of the system operator over the set $\Omega$ whilst remaining stealthy. This strategy of maximum expected loss can be defined as follows.
\begin{definition}[Maximum Expected Loss]\label{expect_loss_def}
The maximum expected loss ($\mathbb{EL}$) associated with the impact-random variable $X^A(\cdot)$, defined in \eqref{newimpact}, is defined as
\begin{align}
\mathbb{EL}[X^A] \triangleq \sup_{a \in \ell_{2e}} \; \mathbb{E}_{\Omega} \left\{X^A(\delta,a) \right\},
\end{align}
where $X^A(\delta,a)$ is the loss on scenario $\delta$ and $\mathbb{E}_{\Omega}$ represents the expectation operator over the set $\Omega$.$\hfill\triangleleft$
\end{definition}

Thus by determining the attack vector that solves for maximal expected loss, one can ensure that the attack vector is stealthy with respect to all uncertainties whilst maximizing the performance loss. Using \textit{Definition \ref{expect_loss_def}}, the risk associated with the impact caused by a bounded-\textbf{R}ational \textbf{A}dversary can be characterized as

\begin{equation}\label{w1}
\begin{aligned}
\gamma_{RA}  &\triangleq \sup_{a \in \ell_{2e}} \mathbb{E}_{\Omega} (X^A(a,\delta)).
\end{aligned}
\end{equation}

Since the operator $\mathbb{E}$ in \eqref{w1} operate over the continuous space $\Omega$, \eqref{w1} is computationally intensive or in general NP-hard \cite[Section 3]{ben1998robust}. Besides, the problem is also non-convex. In the remainder of this article, we discuss methods to solve the optimization problem approximately and efficiently.
\subsection{Discrete uncertainty set}

To begin with, we consider an uncertainty set $\Omega$ which is discrete. The results of this section is the basis for addressing a continuous uncertainty set. For brevity, given a sampled uncertainty $\delta_i \in \Omega$, we define $\tilde{\Sigma}_{p,i} \triangleq ({A}_{cl,i}, {B}_{cl,i}, {C}_{p,i}, {D}_{p,i})$ and  $\tilde{\Sigma}_{r,i} \triangleq ({A}_{cl,i}, {B}_{cl,i}, {C}_{r,i}, {D}_{r,i})$ with $y_p(\delta_i)=y_{pi},y_r(\delta_i)=y_{ri}$ and $x(\delta_i)=x_{i}$ as the outputs and states of $\tilde{\Sigma}_{p,i}$ and  $\tilde{\Sigma}_{r,i}$ correspondingly. For such an isolated uncertainty, the optimization problem \eqref{eqq4} can be rewritten as

\begin{equation}\label{pp1}
\begin{aligned}
 ||\tilde{\Sigma}(\delta_i)||_{\ell_{2e},y_p \leftarrow y_r}^2 &\triangleq  \sup_{a_i \in \ell_{2e}}\; X^A(a_i,\delta_i)\\
 & \begin{aligned}
  =\sup_{a_i \in \ell_{2e}} & \Vert y_{p,i}\Vert_{\ell_2}^2\\
\textrm{s.t.} \; & \Vert y_{r,i}\Vert_{\ell_2}^2 \leq 1,x_i[\infty]=0,
\end{aligned}
\end{aligned}
\end{equation}
where $a_i$ is the attack vector corresponding to the uncertainty $\delta_i$. Using the Lagrange dual function, the non-convex optimization problem \eqref{pp1} can be rewritten into a convex dual problem. This equivalency is depicted in the following lemma.

\begin{lemma}\label{dual_OA}
The optimization problems \eqref{pp1} and \eqref{t5} are equivalent.
\begin{equation}\label{t5}
\begin{aligned}
||\tilde{\Sigma}(\delta_i)||_{\ell_{2e},y_p \leftarrow y_r}^2 = \inf_{\gamma_i \geq 0 } \; & \gamma_i \\
\textrm{s.t.} \; & ||y_{p,i}||^2_{\ell_{2}} - \gamma_i||y_{r,i}||_{\ell_2}^2 \leq 0, \forall a_i \in \ell_{2e}\\
& x_i[0]=0,x_i[\infty]=0,
\end{aligned}
\end{equation}
\end{lemma}
\begin{proof}
The proof is postponed to Appendix \ref{proof_remark}.
\end{proof}

The constraint of the optimization problem \eqref{t5} lies in the signal space which is infinite-dimensional. It would be convenient if the inequality is in the form of a finite-dimensional matrix inequality. To this end, following \cite{Truncated} (also see Appendix \ref{diss_app} \textit{Proposition  \ref{prop1}, $3)$}), the constraint of \eqref{t5} can be characterized as a matrix inequality such that the following conditions are equivalent:
\begin{enumerate}
    \item Given that $x_i[0]=x_i[\infty]=0$, it holds $\forall a_i \in \ell_{2e}$ that \[||y_{p,i}||_{\ell_2}^2 - \gamma_i ||y_{r,i}||_{\ell_2}^2 \leq 0. \]
    \item $\exists P_i=P_i^T \; \text{s.t.}\; M(\gamma_i,\delta^i,P_i) \preceq 0,$ where
\end{enumerate}
\begin{multline*}
M (\gamma_i,\delta^i,P_i) \triangleq
\begin{bmatrix} A_{cl,i}^TP_iA_{cl,i}-P_i & A_{cl,i}^TP_iB_{cl,i} \\ B_{cl,i}^TP_iA_{cl,i} & B_{cl,i}^TP_iB_{cl,i} \end{bmatrix}
+\\
\begin{bmatrix}
C_{p,i}^T\\
D_{p,i}^T
\end{bmatrix}
\begin{bmatrix}
C_{p,i} & D_{p,i}
\end{bmatrix} - \gamma_i
\begin{bmatrix}
C_{r,i}^T\\
D_{r,i}^T
\end{bmatrix}
\begin{bmatrix}
C_{r,i} & D_{r,i}
\end{bmatrix}.
\end{multline*}
Using this equivalency, \eqref{t5} is reformulated as a convex SDP in the following lemma.
\begin{lemma}\label{lemm1}
The optimization problem \eqref{t5} and the following optimization problem are equivalent.
\begin{equation}\label{s2}
\begin{aligned}
||\tilde{\Sigma}(\delta_i)||_{\ell_{2e},y_p \leftarrow y_r}^2 = \min_{\gamma_i \geq 0 ,P_i = P_i^T} \; & \gamma_i \\
\textrm{s.t.} \; & M (\gamma_i,\delta^i,P_i) \preceq 0 \end{aligned}
\end{equation}
$\hfill\square$
\end{lemma}

Next, we discuss the the conditions for boundedness of $||\tilde{\Sigma}(\delta_i)||_{\ell_{2e},y_p \leftarrow y_r}^2$ in the following lemma.

\begin{lemma}[Boundedness]\label{bound_OA_single}
Consider the closed-loop system \eqref{system:CL:uncertain} with an uncertainty $\delta_i$ which is known to the adversary. Then, the optimal value of \eqref{s2} is bounded if and only if one of the following condition hold:
\begin{enumerate}
    \item The system $\tilde{\Sigma}_{r,i}$ has no zeros on the unit circle.
    \item The zeros on the unit circle of the system $\tilde{\Sigma}_{r,i}$ (including multiplicity and input direction) are also zeros of $\tilde{\Sigma}_{p,i}$.
\end{enumerate}
\end{lemma}
\begin{proof}
See Appendix \ref{bound_1}.
\end{proof}

In words, the above lemma states that the attack impact caused by an omniscient adversary on the closed-loop \eqref{system:CL:uncertain} with an uncertainty $\delta_i$ is bounded if either, there does not exist an attack vector which makes the output $y_r$ identically zero, or for all attack vectors which yields $y_r$ identically $0$, it also yields $y_p$ identically zero. To this end, we have formulated the results on characterizing the attack impact by a convex SDP and the condition for its boundedness when the uncertainty set is discrete. These results are the basis for addressing a continuous uncertainty set. Next, the approach discussed in this section is extended to the risk assessment framework when considering a continuum of uncertainties.

\subsection{Continuous uncertainty set}

The optimization problem \eqref{s2} can be directly extended to solve \eqref{uncertain:system} only when we consider a discrete set $\Omega$. But the optimization problem \eqref{uncertain:system} that we intend to solve has a set $\Omega$ which is a continuum. Using the framework of scenario-based reliability estimation \cite{calafiore2007probabilistic}, $\Omega$ can be approximated with a finite set. With this framework of scenario-based reliability estimation, we revisit \eqref{uncertain:system} in the theorem below.

\begin{theorem}\label{Thm1}
Let $\epsilon_1 \in (0,1)$ represent the accuracy with which the probability operator $\mathbb{P}_{\Omega}$ in \eqref{uncertain:system} is to be approximated. Let $\beta_1 \in (0,1)$ represent the confidence with which the accuracy $\epsilon_1$ is guaranteed, i.e.,
\[ \mathbb{P}\{ \vert \mathbb{P}_{\Omega}(||\tilde{\Sigma}(\delta)||_{\ell_{2e},y_p \leftarrow y_r}^2 \leq \gamma ) - \hat{\mathbb{P}}_{N_1} \vert \geq\epsilon_1\} \leq \beta_1.\]
Here $ \hat{\mathbb{P}}_{N_1}$ represents the approximation of the probability operator $\mathbb{P}_{\Omega}$ in \eqref{uncertain:system} defined as
\[
\hat{\mathbb{P}}_{N_1} \triangleq \frac{1}{N_1}\sum_{i=1}^{N_1}\mathbb{I}\left( ||\tilde{\Sigma}(\delta_i)||_{\ell_{2e},y_p \leftarrow y_r}^2 \leq \gamma \right),
\]
where
\begin{equation}\label{sam1}
N_1 \geq \frac{1}{2\epsilon_1^2}\text{log}\frac{2}{\beta_1}.
\end{equation}
Then, the $\text{VaR}_{\beta}$ defined in \eqref{uncertain:system} can be obtained with an accuracy $\epsilon_1$ and confidence $\beta_1$ by solving
\begin{equation}\label{problem:decoupled}
\begin{aligned}
\min & \;\;\gamma\\
\textrm{s.t.} \quad & \frac{1}{N_1} \sum_{i=1}^{N_1} \mathbb{I}\left( \gamma_i \leq \gamma \right) \geq 1-\beta,
\end{aligned}
\end{equation}
where $\mathbb{I}$ is the indicator function. The value of  $\gamma_i\;\forall i \in \{1,\dots, N_1\}$ is obtained by solving the convex optimization problem \eqref{s2}.
\end{theorem}
\begin{proof}
See Appendix \ref{App1}
\end{proof}

\textit{Theorem \ref{Thm1}} states that, to solve \eqref{problem:decoupled}, one could solve the optimization problem \eqref{s2} or equivalently \eqref{t5} for $N_1$ unique realizations of the uncertainty. Since strong duality holds between \eqref{t5} and \eqref{pp1}, the optimal value of the dual optimization problem \eqref{problem:decoupled} indeed provides the VaR$_{\beta}$ with an accuracy $\epsilon_1$. Next, by extending \textit{Lemma \ref{bound_OA_single}}, the condition for boundedness of \eqref{problem:decoupled} is stated in the following lemma.

\begin{lemma}[Boundedness]\label{bound_OA_multi}
Consider $N_1$ independent and identically distributed (i.i.d) realizations of the closed-loop system \eqref{system:CL:uncertain}, each with an uncertainty $\delta_i$. The optimal value of \eqref{problem:decoupled} with these $N_1$ system realizations is bounded if and only if the optimal value of \eqref{s2} is bounded for at least $\left \lceil{N_1(1-\beta)}\right \rceil $ system realizations.
\end{lemma}
\begin{proof}
Necessity: Let us multiply both sides of the constraint of \eqref{problem:decoupled} by $N_1$. Then, for a bounded solution $\gamma_{OA}$ to exist, we require that
\begin{equation}
\sum_{i=1}^{N_1} \mathbb{I}\left( \gamma_i \leq \gamma_{OA} , \; \forall a_i \in \ell_{2e}\right) \geq N_1-N_1\beta,
\end{equation}
should hold with a bounded $\gamma_{OA}$. This is satisfied only if $
\lceil N_1(1-\beta)\rceil$ values of the set $\{\gamma_i\}_{i=1,\dots,N_1}$ are bounded.\\
Sufficiency: Assume that the values of the set $\{\gamma_i\}_{i=1,\dots,N_1}$, for the first $\left \lceil N_1(1-\beta)\right \rceil$ realizations is bounded and is unbounded for the other realizations. Let us choose $\gamma_{OA} = \max_{i=1,\dots,\lceil N_1(1-\beta) \rceil} \gamma_i$. Substituting the definition of $\gamma_i, \gamma_{OA}$ in the constraint of \eqref{problem:decoupled} yields
\begin{align*}
&\frac{1}{N_1}\left\{\sum_{i=1}^{\lceil N_1(1-\beta)\rceil} \mathbb{I}\left( \gamma_i \leq \gamma_{OA} \right) + \sum_{i=\lceil N_1(1-\beta)\rceil +1}^{N_1} \mathbb{I}\left( \gamma_i \leq \gamma_{OA} \right) \right\} \\
&=\frac{1}{N_1}\left\{ \lceil N_1(1-\beta) \rceil + 0 \right\} \geq 1-\beta.
\end{align*}
Thus, there exists a bounded $\gamma_{OA}$ which satisfies the constraint of \eqref{problem:decoupled}. This concludes the proof.
\end{proof}

The interpretation of \textit{Lemma \ref{bound_OA_multi}} is that the system operator tolerates a fraction ($\beta \times 100 \%$) of cases where the impact \eqref{s2} is unbounded. Conversely, even-though the impact is unbounded for certain elements of the uncertainty set, the risk assessed by the system operator will still be bounded. This allows the system operator to be less pessimistic.

In this section, we have discussed a scenario-based approach to determine the VaR$_{\beta}$. The risk assessment problem was posed as a convex SDP. However, the number of optimization problems necessary to determine the risk is $N_1$. Next, we discuss an approach to solve the optimization problem \eqref{w1}.
\subsection{Approximating the uncertainty set }
To recall, we are interested in determining the maximum expected loss associated with the impact caused by a rational adversary. Unfortunately, this problem is computationally intensive or in general NP-hard. Thus, as a first step toward solving \eqref{w1}, we approximate the objective function in \textit{Lemma \ref{Expectation_sample_lemma}}.
\begin{lemma}\label{Expectation_sample_lemma}
Let $\delta_i$ be sampled uncertainties from $\Omega$. Let us define
\begin{equation}
    \hat{\mathbb{E}}_{\Omega_{N_2}}(X^A(a,\delta)) \triangleq \frac{1}{N_2}\sum_{i=1}^{N_2} X^A(a_i,\delta_i).
\end{equation}
Then, it holds that $\lim_{N_2 \to \infty} \hat{\mathbb{E}}_{\Omega_{N_2}}(X^A(a_i,\delta_i)) = \mathbb{E}_{\Omega}(X^A(a,\delta)).$
\end{lemma}
\begin{proof}
The proof follows from applying \cite[{Theorem 7.2}]{tempo2012randomized} to approximate the expectation operator in \eqref{w1}.
\end{proof}

\textit{Lemma \ref{Expectation_sample_lemma}} states that the continuous set $\Omega$ can be approximated by a discrete set $\Omega_{N_2}$ of cardinality $N_2$. The approximation becomes more accurate as $N_2 \to \infty$. This approximation introduces a curse of dimensionality since, to obtain a good estimate of the risk and to obtain a well feasible attack vector, the required number of samples is too large for practical implementation. To circumvent this practical issue, we next show that an attack vector obtained by solving \eqref{w1} with a discrete uncertainty set as mentioned in \textit{Lemma \ref{Expectation_sample_lemma}} is partially feasible to the original optimization problem \eqref{w1} with a continuous uncertainty set.

It might not be immediately apparent that the notion of feasibility applies to \eqref{w1} since there are no external constraints present. Thus, we begin by simplifying the optimization problem \eqref{w1}.
\begin{lemma}\label{pp31}
The optimization problem \eqref{w1} is equivalent to the optimization problem
\begin{equation}\label{pp3}
\begin{aligned}
 \sup_{a \in \ell_{2e},\beta \in [0,\;1]} \; & {\mathbb{E}_{\Omega}} \left\{||y_p(\delta)||_{\ell_2}^2\; \vert \; \begin{pmatrix}
||y_r(\delta)||_{\ell_2}^2 \leq1\\
x(\delta)[\infty]=0
\end{pmatrix} \right\} (1-\beta) \\
\textrm{s.t.} \quad & {\mathbb{P}_{\Omega}}
\begin{pmatrix}
||y_r(\delta)||_{\ell_2}^2 \leq1\\
x(\delta)[\infty]=0
\end{pmatrix}\geq 1-\beta.
\end{aligned}
\end{equation}
\end{lemma}
\begin{proof}
Consider the function $X^A(a,\delta)$ in \eqref{newimpact}. By expanding its indicator function, we can write $\mathbb{E}_{\Omega}(X^A(a,\delta))$ as
\begin{align}\label{w3}
\nonumber & {\mathbb{E}_{\Omega}} \left\{ ||y_p(\delta)||_{\ell_2}^2\; \vert \begin{pmatrix}
||y_r(\delta)||_{\ell_2}^2 \leq1\\
x(\delta)[\infty]=0
\end{pmatrix} \right \}{\mathbb{P}_{\Omega}}\begin{pmatrix}
||y_r(\delta)||_{\ell_2}^2 \leq1\\
x(\delta)[\infty]=0
\end{pmatrix}
\end{align}
Using the above definition in \eqref{w1}, we obtain $\gamma_{RA}$ as
\begin{equation}
\sup_{a \in \ell_{2e}} {\mathbb{E}_{\Omega}} \left\{ ||y_p(\delta)||_{\ell_2}^2\; \vert \;\begin{pmatrix}
||y_r(\delta)||_{\ell_2}^2 \leq1\\
x(\delta)[\infty]=0
\end{pmatrix} \right \}{\mathbb{P}_{\Omega}}\begin{pmatrix}
||y_r(\delta)||_{\ell_2}^2 \leq1\\
x(\delta)[\infty]=0
\end{pmatrix}
\end{equation}
which can be rewritten as \eqref{pp3}. This concludes the proof.
\end{proof}

\textit{Lemma \ref{pp31}} uncovers the constraints present in the optimization problem \eqref{w1}. We can now discuss the notion of feasibility in regards to the optimization problem \eqref{pp3}. So, we continue by simplifying \eqref{pp3} as follows. In reality, $\beta$ represents the fraction of the uncertainty set with respect to which the adversary is not stealthy. Let us assume that the adversary is injecting a maximally robust attack, i.e., the attack vector is stealthy w.r.t. all uncertainties. Thus we could set $\beta =0$. Motivated by the above argument, we rewrite \eqref{pp3} as
\begin{equation}\label{fi1}
\begin{aligned}
\sup_{a \in \ell_{2e}} \; & {\mathbb{E}_{\Omega}} \left\{||y_p(\delta)||_{\ell_2}^2\; \vert \;\begin{pmatrix}
||y_r(\delta)||_{\ell_2}^2 \leq1\\
x(\delta)[\infty]=0
\end{pmatrix}\right\} \\
\textrm{s.t.} \quad &
\begin{pmatrix}
||y_r(\delta)||_{\ell_2}^2 \leq1\\
x(\delta)[\infty]=0
\end{pmatrix} \forall \delta \in \Omega
\end{aligned}
\end{equation}
Recalling the approximation result in \textit{Lemma \ref{Expectation_sample_lemma}}, assume that $\Omega$ is replaced with a discrete uncertainty set $\Omega_{N}$, so that \eqref{fi1} is approximated by
\begin{equation}\label{uncertain:adversary}
\begin{aligned}
\gamma_{RA}(N_2) = \sup_{a \in \ell_{2e}} & \hat{\mathbb{E}}_{\Omega_{N_2}} \left\{||y_p(\delta)||_{\ell_2}^2\; \vert \;\begin{pmatrix}
||y_r(\delta)||_{\ell_2}^2 \leq1\\
x(\delta)[\infty]=0
\end{pmatrix}\right\} \\
\textrm{s.t.} &
\begin{pmatrix}
||y_r(\delta)||_{\ell_2}^2 \leq1\\
x(\delta)[\infty]=0
\end{pmatrix} \forall \delta \in \Omega_{N_2}.
\end{aligned}
\end{equation}

Let the resulting optimal attack vector be denoted by $a_{N_2}^*$. Then the following theorem provides \textit{a posteriori} results on the feasibility of the attack vector $a_{N_2}^*$ to \eqref{fi1}.

\begin{theorem}\label{Thm_sampling}
Let the number of samples $N_2$ and the confidence level $\lambda \in (0,1)$ be predefined constants. Define $\epsilon(\cdot)$ such that
\begin{equation}\label{sample:constraint2}
\epsilon(N_2) = 1 \;,\; \sum_{k=0}^{N_2-1} {N_2\choose k}(1 - \epsilon(k))^{N_2-k} = \lambda.
\end{equation}
Let $s_{N_2}^*$ represent the cardinality of the support subsample for $(\delta_1,\dots,\delta_{N_2})$ (see \cite[\textit{Definition 2}]{campi2018general}). Then it holds that
\begin{equation}
\mathbb{P}^{N_2}\{ \mathbb{P}_{\Omega}\{\delta \in \Omega \;\vert\; a_{N_2}^* \notin \Theta\}> \epsilon(s_{N_2}^*)\}\leq \lambda,
\end{equation}
where $a_{N_2}^*$ is the argument of the optimization problem \eqref{uncertain:adversary} and $\Theta$ is defined as
\begin{align}\label{theta}
    \Theta &\triangleq \bigcap_{\delta \in \Omega} \Theta_{\delta}, \text{where} \;\Theta_{\delta} \triangleq \left\{ a \in \ell_{2e} \vert
\begin{pmatrix}
||y_r(\delta)||_{\ell_2}^2 \leq1\\
x(\delta)[\infty]=0
\end{pmatrix}
\right\}
\end{align}
\end{theorem}
\begin{proof}
See Appendix \ref{Main_thm}.
\end{proof}

In words, \textit{Theorem \ref{Thm_sampling}} states that the attack vector obtained by solving \eqref{uncertain:adversary} is stealthy and state-vanishing for all the closed-loop system of the form \eqref{system:CL:uncertain} with uncertainties belonging to the set $\Omega$ except for the fraction $\epsilon(s_{N_2}^*)$ of $\Omega$. It also states that the accuracy ($\epsilon(\cdot)$) and confidence ($\lambda$) of the solution are independent of the distribution. This result is a direct consequence of \cite[\textit{Theorem 1}]{campi2018general}.

In conclusion, it follows that \eqref{fi1} can be solved approximately with a discrete set $\Omega_{N_2}$ of arbitrary but bounded cardinality. Thus, the next section will focus on solving the optimization problem \eqref{uncertain:adversary}.
\subsection{Risk assessment}
The optimization problem \eqref{uncertain:adversary} has two main disadvantages namely $(i)$ it is a non-convex optimization problem, and $(ii)$ its constraints lie in the infinite-dimensional attack space. To resolve these disadvantages, we adopt S-procedure and dissipative system theory and revisit the optimization problem \eqref{uncertain:adversary} in the theorem below. To begin with, given a sampled uncertainty $\delta_i \in \Omega$, we define $\tilde{\Sigma}_{p,i} \triangleq ({A}_{cl,i}, {B}_{cl,i}, {C}_{p,i}, {D}_{p,i})$ and  $\tilde{\Sigma}_{r,i} \triangleq ({A}_{cl,i}, {B}_{cl,i}, {C}_{r,i}, {D}_{r,i})$ with $y_p(\delta_i)=y_{pi},y_r(\delta_i)=y_{ri}$ and $x(\delta_i)=x_{i}$ as the outputs and states of $\tilde{\Sigma}_{p,i}$ and  $\tilde{\Sigma}_{r,i}$ correspondingly. Now we are ready to present the main theorem of this section.

\begin{theorem}\label{Thm2}
Let $N_2$ be a predefined constant. The maximum expected loss \eqref{uncertain:adversary} associated with the impact caused by a rational adversary injecting a maximally robust attack vector on \eqref{system:CL:uncertain} can be obtained by solving the convex SDP \footnote{With an abuse of notation, we denote that every element of the vector $\gamma$ is non-negative by $\gamma \geq 0$}
\begin{equation}\label{problem:coupled:final}
\begin{aligned}
\min_{\gamma \geq 0, P = P^T}  & \mathbf{1}^T\begin{bmatrix}
\gamma_1 & \dots & \gamma_{N_2}
\end{bmatrix}^T \\
\textrm{s.t.} \quad &
\begin{bmatrix}
\Bar{A}^TP\Bar{A}-P &\Bar{A}^TP\Bar{B}\\ \Bar{B}^TP\Bar{A}&\Bar{B}^TP\Bar{B}\end{bmatrix}+ \Psi(\gamma) \preceq 0,
\end{aligned}
\end{equation}
where $\Psi(\gamma) \triangleq \frac{1}{N_2}
\begin{bmatrix}
\bar{C}_{p}^T\\
\bar{D}_{p}^T
\end{bmatrix}
\begin{bmatrix}
\bar{C}_{p} & \bar{D}_{p}
\end{bmatrix}
- \begin{bmatrix}
\bar{C}_{r}^T\\
\bar{D}_{r}^T
\end{bmatrix}\Gamma(\gamma)
\begin{bmatrix}
\bar{C}_{r} & \bar{D}_{r}
\end{bmatrix}$,
\begin{align}\label{balaji}
\left[
\begin{array}{c|c}
\bar{A} & \bar{B}\\
\hline
\\[\dimexpr-\normalbaselineskip+2pt] \bar{C}_p & \bar{D}_p\\
\hline
\\[\dimexpr-\normalbaselineskip+2pt] \bar{C}_r & \bar{D}_r
\end{array}
\right] =
\left[
\begin{array}{c|c}
\begin{matrix}
A_{cl,1} & \dots & 0 \\
\vdots & \ddots & \vdots\\
0 & \dots & A_{cl,N_2}
\end{matrix}     & \begin{matrix}
B_{cl,1} \\ \vdots \\ B_{cl,N_2}
\end{matrix} \\
\hline
\\[\dimexpr-\normalbaselineskip+2pt] \begin{matrix}
C_{p,1} & \dots & 0\\
\vdots  & \ddots & \vdots \\
0 & \dots & C_{p,N_2}
\end{matrix} & \begin{matrix}
D_{p,1}^T \\ \vdots \\ D_{p,N_2}^T
\end{matrix}\\
\hline
\\[\dimexpr-\normalbaselineskip+2pt] \begin{matrix}
C_{r,1} & \dots & 0\\
\vdots  & \ddots & \vdots \\
0 & \dots & C_{r,N_2}
\end{matrix} & \begin{matrix}
D_{r,1}^T \\ \vdots \\ D_{r,N_2}^T
\end{matrix}
\end{array}
\right],
\end{align}
and $\Gamma(\gamma) =I_{n_r} \bigotimes \text{diag}(\gamma_1,\dots.\gamma_{N_2})$. The dimension of each matrix is given in TABLE \ref{matrix_dimention}.
\begin{table}
\centering
\begin{tabular}{|| m{3em} | m{6em}|| m{3em} | m{6em}||}
 \hline
Matrix & Dimension & Matrix & Dimension  \\
\hline
 $\bar{A}$ &  $n_cN_2 \times n_cN_2$ & $\bar{B}$ & $n_cN_2 \times n_a$\\
 $\bar{C}_p$ & $n_pN_2\times n_cN_2$ & $\bar{D}_p$ & $n_pN_2\times n_a$\\
 $\bar{C}_r$ & $n_rN_2\times n_cN_2$ & $\bar{D}_r$ & $n_rN_2\times n_a$\\
 \hline
\end{tabular}
\caption{Dimension of matrices,$\;n_c\triangleq n_x+n_z+n_s$}
\label{matrix_dimention}
\end{table}
\end{theorem}
\begin{proof}
See Appendix \ref{App2}.
\end{proof}

The optimization problem \eqref{uncertain:adversary} is the primal problem with its optimizer being the attack vector $a$. This optimization problem is non-convex. With the help of S-procedure and dissipative system theory, \eqref{uncertain:adversary} is converted to its equivalent dual SDP form \eqref{problem:coupled:final} with its optimizer $\gamma, P$, which is convex. This equivalency also helps us to conclude that the duality gap is zero. The necessary and sufficient conditions for the value of \eqref{problem:coupled:final} to be finite is given \textit{Lemma \ref{Bound_RA}}.

\begin{lemma}[Boundedness]\label{Bound_RA}
Consider $N_2$ i.i.d. realizations of the closed-loop system \eqref{system:CL:uncertain}, each with an uncertainty $\delta_i$. The optimal value of \eqref{problem:coupled:final} with the above mentioned system realizations is bounded if and only if the system with $\bar{\Sigma}_p=(\Bar{A},\Bar{B},\bar{C}_p,\bar{D}_p)$ and $\bar{\Sigma}_r=(\Bar{A},\Bar{B},\bar{C}_r,\bar{D}_r)$ satisfy one of the following:
\begin{enumerate}
    \item The system $\bar{\Sigma}_r$ has no zeros on the unit circle.
    \item The zeros on the unit circle of the system $\bar{\Sigma}_r$ (including multiplicity and input direction) are also zeros of $\bar{\Sigma}_p$.
\end{enumerate}
\end{lemma}
\begin{proof}
See Appendix \ref{bound_2}.
\end{proof}
Let the outputs of $\bar{\Sigma}_p$ and $\bar{\Sigma}_r$ be represented by $\bar{y}_r$ and $\bar{y}_p$ respectively. Then, in words, \textit{Lemma \ref{Bound_RA}} states that the maximum expected loss \eqref{uncertain:adversary} is bounded if either, there does not exist an attack vector which makes the output $\bar{y}_r$ identically zero, or for all attack vectors which yields $\bar{y}_r$ identically $0$, it also yield $\bar{y}_p$ identically zero.

It is important to study the conditions for unbounded risk because, if the conditions of the lemma do not hold, it means that there exists an attack vector that can cause a very huge impact whose upper-bounded cannot be determined and remain stealthy. However, as the conditions of the lemma are necessary and sufficient, the operator can alter the system matrices so that the conditions hold and consequently reduce the vulnerability of the system to such attacks. In the next section, we illustrate the results with a numerical example.

\subsection{Risk assessment for a rational adversary}\label{EXMP-2}

In view of \textit{Lemma \ref{Expectation_sample_lemma}}, we choose $N_2=21$ to approximate the set $\Omega$. With this approximation, by solving the convex SDP \eqref{problem:coupled:final}, we obtain $\gamma_{RA}(N_2) = 617.267$. To recall, $\gamma_{RA}$ represents the maximum expected performance loss of the system operator. In this implementation, the value of risk was obtained for $\beta=0$ since we considered a maximally robust adversary.

Next, we discuss the validity of the approximation $\Omega_{N_2}$. For varying values of $N_2$, the number of non-zero $\gamma_{i}$s obtained while solving \eqref{problem:coupled:final} is shown in the first two columns of TABLE \ref{guarentees_table}. In view of \textit{Theorem \ref{Thm_sampling}}, if we solve the problem \eqref{problem:coupled:final} with an arbitrary $N_2$, we can provide guarantees on the optimization problem \eqref{w1} as shown in column $3$ of TABLE \ref{guarentees_table}. Here, the function $\epsilon(\cdot)$ is evaluated according to \cite[(7)]{campi2018general} where $\lambda=10^{-2}$ and $s_{N_2}^* = |\text{supp}(\gamma)|$\footnote{The confidence is denoted by $\lambda$ here whereas it is denoted by $\beta$ in \cite{campi2018general}}. And $\epsilon(\cdot)$ should be interpreted as follows: The attack vector obtained by solving \eqref{problem:coupled:final} with $N_2$ samples, with probability $1-\lambda$, will be at-most feasible for the fraction $(1-\epsilon(\cdot))\Omega$ of the set $\Omega$.

Since $\epsilon(\cdot)$ represents the fraction of the uncertainty set to which the attack vector is infeasible, it is intuitive to expect this value to be close to zero. Numerically we have observed from TABLE \ref{guarentees_table} that $\text{supp}(\gamma)$ is always one. So, assuming again that $s_{N_2}^*=1$, the number of samples required to guarantee that the attack vector, with a probability $1-\lambda$, is feasible for $1-\epsilon(\cdot)$ of the uncertainty set can be obtained by picking $N_2$ such that \eqref{N2_closed_form} holds \cite[(7)]{campi2018general}
\begin{equation}\label{N2_closed_form}
  \sqrt[\leftroot{-2}\uproot{2}N_2-1]{\frac{\lambda}{N_2^2}} = 1-\epsilon(\cdot).
\end{equation}
Thus \eqref{N2_closed_form} gives an idea of how $N_2$ increases as $\epsilon(\cdot)$ decreases. And consequently gives an idea of the scalability of the proposed approach as the dimension of the matrix inequality \eqref{problem:coupled:final} depends on $N_2$. We also depict the computation complexity by providing the time taken to solve \eqref{problem:coupled:final} in the last column of TABLE \ref{guarentees_table}.

\begin{table}
\centering
\begin{tabular}{|| c | c | c | c | c||}
 \hline
 $N_2$ & $\vert\text{supp}(\gamma)\vert$ & $\epsilon(s_{N_2}^*)$ & $\gamma_{RA}(N_2)$ & Time (seconds)\\
 \hline
 8 & 1 & 0.7141 & 638.04 & 17\\
 \hline
15 & 1 & 0.5112 & 628.60 & 116\\
 \hline
 21 & 1 & 0.4142 & 617.26 & 578\\
 \hline
\end{tabular}
\caption{Rational adversary - a posteriori Guarantees}
\label{guarentees_table}
\end{table}

\section*{Appendix}
\section{Explaining the performance criterion}\label{app_explain}
Consider an LTI system which was described in \eqref{P} without attack. Let us consider that, the operator is evaluating the performance of the system by a linear-quadratic (LQ) cost with cost matrix $Q \succeq 0$ on the states and  cost matrix $R \succ 0$ on the control input. Here, with a slight loss of generality, we consider that the cost on the cross-term is zero. Next, let us define
    \begin{equation}
        y_p[k] = \underbrace{\begin{bmatrix}
        \sqrt{Q} \\ 0
        \end{bmatrix}}_{C_J} x[k] + \underbrace{\begin{bmatrix}
        0 \\ \sqrt{R}
        \end{bmatrix}}_{D_J} u[k] = \begin{bmatrix}
        \sqrt{Q} & 0\\
        0 & \sqrt{R}
        \end{bmatrix} \begin{bmatrix}
        x[k] \\ u[k]
        \end{bmatrix}
    \end{equation}
    Here, since the matrix $Q$ and $R$ are positive definite and positive semi-definite, their square root exists. Then, evaluating this LQ cost is equivalent to evaluating the $\ell_2$ norm of $y_p[k]$. That is
    \begin{equation*}
    ||y_p||_{\ell_2}^2 \triangleq
    \sum_{k=0}^{N} \begin{bmatrix}
    x[k]^T & u[k]^T
    \end{bmatrix}
    \begin{bmatrix}
        \sqrt{Q} & 0\\
        0 & \sqrt{R}
        \end{bmatrix}
        \begin{bmatrix}
        \sqrt{Q} & 0\\
        0 & \sqrt{R}
        \end{bmatrix} \begin{bmatrix}
        x[k] \\ u[k]
        \end{bmatrix}
    \end{equation*}
    Thus we depict $(i)$ how to choose the nominal matrices $C_J,$ and $D_J$, and $(ii)$ the reason to evaluate the attack impact in terms of $||y_p||_{\ell_2}^2$. Similarly, we assume that faults which have a very low effect on the detection output can occur rarely (say faults which change the energy of the detection output to $\epsilon_r$). Thus an attack is considered to be detected when $||y_p||_{\ell_2}^2 > 1$.

\section{Proof of \textit{Theorem \ref{Thm_sampling}}}\label{Main_thm}
Before presenting the proof, an introduction to scenario-based constraint satisfaction is provided.
\subsection*{Scenario-based constraint satisfaction \cite{campi2018general}}
Consider the constrained non-convex optimization problem $\inf_{\theta\in \Theta}  f(\theta,\delta),$ where $\delta \in \Delta$ is the uncertainty and $\theta$ is the infinite dimensional decision parameter which lies in the set $\Theta \triangleq \bigcap_{\delta \in \Delta} \Theta_{\delta}$, where $\Theta_{\delta}$ is the constraint set which includes all the admissible parameters for the isolated uncertainty $\delta$.

\begin{definition}[Violation probability]
Let us define the violation probability as $\mathbb{V}(\theta) \triangleq \mathbb{P}_{\Delta}\{\delta \in \Delta \;|\;\theta \notin \Theta\}.$ $\hfill\triangleleft$
\end{definition}

\begin{definition}[$\epsilon$-level solution]
Let $\epsilon \in (0,1)$. Then, $\theta \in \Theta$ is an $\epsilon$ level solution if $\mathbb{V}(\theta) \leq \epsilon. \hfill\triangleleft$
\end{definition}

\begin{definition}[Confidence level : $\lambda$]
Let $\lambda \in (0,1)$. Then, the confidence level $\lambda$ represents the probability that $\theta$ is not an $\epsilon$ level solution. i.e., $\lambda \triangleq \mathbb{P}\{\mathbb{V}(\theta) > \epsilon\}.\hfill\triangleleft$
\end{definition}

\begin{proof}[Proof of \textit{Theorem \ref{Thm_sampling}}]
The optimization problem \eqref{fi1} can be reformulated as
\begin{equation}\label{w7}
    -\inf_{a \in \ell_{2e}} \left\{ \mathbb{E}_{\Omega} \{-||y_p(\delta)||_{\ell_2}^2 \}\left. \right|\  \begin{pmatrix}
||y_r(\delta)||_{\ell_2}^2 \leq1\\
x(\delta)[\infty]=0
\end{pmatrix} \forall \delta \in \Omega \right\}
\end{equation}
In view of \cite[Theorem 1]{campi2018general}, let us define the objective function as $f(a,\delta) \triangleq \mathbb{E}_{\Omega} \{-||y_p(\delta)||_{\ell_2}^2\},$
where $y_p(\delta)$ is also a function of the attack vector $a$. Let us define the set $\Theta$ as in \eqref{theta}. Let us define a confidence level $\lambda \in (0,1)$, a constant $N_2$ and $\epsilon(\cdot)$ such that \eqref{sample:constraint2} holds. Then applying \cite[Theorem 1]{campi2018general}, we obtain that
\begin{equation}
\mathbb{P}^{N_2}\{\mathbb{P}_{\Omega}\{\delta \in \Omega\;\vert\;a_{N_2}^* \notin \Theta\} > \epsilon(s_{N_2}^*)\}\leq \lambda,
\end{equation}
Thus, with a probability level $1-\lambda$, the solution $a_{N_2}^*$ is $\epsilon(s_{N_2}^*)$ feasible to the optimization problem \eqref{w7}. In our problem setting, $a_{N_2}^*$ is the optimal argument of the optimization problem
\begin{equation}
\begin{aligned}
-arg \inf_{a \in \ell_{2e}} \; & \left \{\frac{-1}{N_2} \sum_{i=1}^{N_2} \{||y_{p,i}||_{\ell_2}^2\; \vert \;
\begin{pmatrix}
||y_{r,i}||_{\ell_2}^2 \leq 1\\
x_{i}[\infty]=0
\end{pmatrix}\} \right\}\\
\textrm{s.t.} \quad & \begin{pmatrix}
||y_{r,i}||_{\ell_2}^2 \leq1\\
x_{i}[\infty]=0
\end{pmatrix}, \forall i \in \{1,\dots,N_2\},
\end{aligned}
\end{equation}
which can be rewritten as \eqref{uncertain:adversary}. This concludes the proof.
\end{proof}
\section{Proof of \textit{Theorem \ref{Thm2}}}\label{App2}
Before presenting the proof, an introduction to S-system and S-procedure is provided as it helps formulating the proof.
\subsection*{S-system \cite[\textit{Definition 4.3.1}]{petersen2012robust}}
Let $\mathcal{L}$ be a real Hilbert space with a well defined inner product denoted by $\langle \cdot,\cdot \rangle$. Let $\mathcal{G}_0(\cdot),\dots,\mathcal{G}_k(\cdot)$ be quadratic functionals mapping $\mathcal{L} \to \mathbb{R}$. Let $\omega$ be a discrete time signal.
\begin{definition}[S-system]\cite[\textit{Definition 4.3.1}]{petersen2012robust}\label{S-system}
The quadratic functionals $\mathcal{G}_0(\cdot),\mathcal{G}_1(\cdot),\dots,\mathcal{G}_k(\cdot)$ form an S-system if there exist a bounded linear operator $\mathbf{T}_i: \mathcal{L} \rightarrow \mathcal{L},i=1,2,\dots,$ such that
\begin{enumerate}
    \item $\langle \mathbf{T}_i\omega_1,\omega_2\rangle \rightarrow 0$ as $i \rightarrow \infty,\;\forall\; \omega_1,\omega_2 \in \mathcal{L}$.
    \item If $\omega \in \mathcal{M},$ then $\mathbf{T}_i\omega \in \mathcal{M},\;\forall i=1,2,\dots$, where $\mathcal{M}$ is a linear subspace of $\mathcal{L}$.
    \item $\mathcal{G}_j(\mathbf{T}_i\omega) \to \mathcal{G}_j(\omega)$ as $i \rightarrow \infty,\;\forall \omega\in \mathcal{L}$ and $j=0,1,\dots,k$.$\hfill\triangleleft$
\end{enumerate}
\end{definition}
We next present a theorem which helps in proving \textit{Theorem \ref{Thm2}}.
\begin{theorem}\label{thm_new}
Let us define a stable discrete time linear time-invariant system of the form
\begin{align}\label{eq1}
    {\eta}[k+1] &= \Phi\eta[k] + \Lambda\mu[k]\\
    \sigma[k] &=\Pi\eta[k] + \Upsilon \mu[k]\qquad \eta[0]=\eta_0, \; \eta[\infty]=0.
\end{align}
Let us define the set $\mathcal{L}$ as
\begin{equation}\label{eq2}
    \mathcal{L} \triangleq \left\{  \omega = \begin{bmatrix}
    \sigma\\\mu
    \end{bmatrix} \vert \begin{bmatrix}
    \sigma, \mu \text{\;are related by\;} \eqref{eq1}\\
    \mu \in \ell_{2e}, \; \eta[0]=\eta_0, \; \eta[\infty]=0.
    \end{bmatrix} \right\}
\end{equation}
Let us also define the functionals $ \mathcal{G}_0(\omega) \triangleq \sum_{k=0}^{\infty} \omega[k]^TM_0\omega[k] + \zeta_0, \dots,   \mathcal{G}_k(\omega) \triangleq \sum_{k=0}^{\infty} \omega[k]^TM_k\omega[k] + \zeta_k$.
where $M_0,\dots,M_k$ are given matrices and $\zeta_0,\dots,\zeta_k$ are given constants. Then, the functionals $-\mathcal{G}_0(\cdot),\dots,\mathcal{G}_k(\cdot)$ form an S-system.
\end{theorem}
\begin{proof}[Proof of \textit{Theorem \ref{thm_new}}]
In view of \textit{Definition \ref{S-system}}, let us define the operator $\mathbf{T}_i$ as
\[\mathbf{T}_i\omega[k]=
\begin{cases}
    0, & \text{if}\;\; 0\leq k \leq i\\
    \omega[k-i],  & \text{if}\;\; k > i
\end{cases}.\]
For any $\omega_1,\omega_2 \in \mathcal{L}$, $\langle \mathbf{T}_i\omega_1,\omega_2 \rangle =  |\sum_{k=0}^\infty \langle \mathbf{T}_i\omega_1[k],\omega_2[k]\rangle|^2 = $
\begin{equation}
|\sum_{k=i}^\infty \langle \omega_1[k-i],\omega_2[k]\rangle|^2 \leq_{1} \sum_{k=i}^{\infty} ||\omega_1[k-i]||^2 \sum_{k=i}^{\infty} ||\omega_2[k]||^2
\end{equation}
where the inequality $1$ holds due to the Cauchy Schwartz inequality. Since the theorem states that $\lim_{k\to \infty}\eta[k] =0$, it holds from \textit{Lemma \ref{vanish_lemma}} that $\lim_{k\to \infty}\mu[k] =0$. Following which, it immediately holds that $\lim_{k\to \infty}\sigma[k] =0$ which implies that $\lim_{k\to \infty}\omega[k] =0$. Due to the above reasoning, it holds that $\lim_{i\to \infty} \sum_{k=i}^{\infty}||\omega_2[k]||^2 \to 0 \; \forall \omega \in  \mathcal{L}$ . Thus condition $1)$ of \textit{Definition \ref{S-system}} holds.

Let us consider a set $\mathcal{M} = \mathcal{L}|_{\eta_0=0}$. Then, if $\omega \in \mathcal{M}$, due to the time-invariant property of \eqref{eq1}, $\mathbf{T}_i\omega \in \mathcal{M}$. This proves that $\exists \;\mathcal{M} \subset \mathcal{L}$ which is invariant under the operator $\mathbf{T}_i$. Thus condition $2)$ of \textit{Definition \ref{S-system}} holds. Let us consider $\mathcal{G}_0(\mathbf{T}_i\omega) = $
\begin{align*}
\sum_{k=0}^{\infty} (\mathbf{T}_i\omega[k])^TM_0\mathbf{T}_i\omega[k] + \zeta_0 = \sum_{k=i}^{\infty} \omega[k-i]^TM_0\omega[k-i] + \zeta_0
\end{align*}
$= \mathcal{G}_0(\omega)$. Similarly, it can be observed that $\mathcal{G}_j(\mathbf{T}_i\omega) = \mathcal{G}_j(\omega)\;\forall j=\{0,\dots,k\}$. Thus condition $3)$ of \textit{Definition \ref{S-system}} holds. Since we have shown that the functionals $-\mathcal{G}_0(\omega(\cdot)),\mathcal{G}_1(\omega(\cdot)),\dots,\mathcal{G}_k(\omega(\cdot))$ satisfy conditions $1),2)$ and $3)$ of \textit{Definition \ref{S-system}}, they form an S-system. This concludes the proof.
\end{proof}
Now we are ready to present the proof of \textit{Theorem \ref{Thm2}}.

\begin{proof}
\subsubsection*{Step 1}[Problem reformulation]
Using the hypograph formulation, \eqref{uncertain:adversary} can be rewritten as
\begin{equation}\label{eqq1}
\sup_{\upsilon, a \in \ell_{2e}} \left\{  \upsilon \; \Bigg|\; \begin{aligned}
& \frac{1}{N_2} \sum_{i=1}^{N_2} ||y_{p,i}||_{\ell_2}^2 \geq \upsilon\\
& ||y_{r,i}||_{\ell_2}^2 \leq 1,\;x_{i}[\infty]=0,\; \forall i \in \{1,\dots,N_2\}
\end{aligned}\right\}
\end{equation}
From now, in the next two steps of the proof, we focus on the optimization problem \eqref{eqq1} without the state constraints
\begin{equation}\label{eqq2}
\sup_{\upsilon, a \in \ell_{2e}} \left\{  \upsilon \; \Bigg|\; \begin{aligned}
& \frac{1}{N_2} \sum_{i=1}^{N_2} ||y_{p,i}||_{\ell_2}^2 \geq \upsilon\\
& ||y_{r,i}||_{\ell_2}^2 \leq 1,\; \forall i \in \{1,\dots,N_2\}
\end{aligned}\right\}
\end{equation}

The reason to focus on \eqref{eqq2} rather than \eqref{eqq1} is that S-procedure becomes convenient to apply when there are no equality constraints. Thus, we drop the state constraints now and introduce them back at the end of \textit{Step 3}. Equivalently \eqref{eqq2} be reformulated as
\begin{equation}\label{eq4}
\inf_{\upsilon}\left\{\upsilon \Bigg|\left\{a\in\ell_{2e}\Big|\;\begin{aligned}
&  \frac{1}{N_2} \sum_{i=1}^{N_2} ||y_{p,i}||_{\ell_2}^2 - \upsilon > 0 \\
& 1-||y_{r,i}||_{\ell_2}^2 \geq 0, \forall i
\end{aligned}\right\}=\emptyset\right\}
\end{equation}
\subsubsection*{Step 2}
Recall that the system matrices for the system with the isolated uncertainty $\delta_i$, with attack vector as input and performance and detection outputs as outputs are
$\tilde{\Sigma}_{p,i} \triangleq ({A}_{cl,i}, {B}_{cl,i}, {C}_{p,i}, {D}_{p,i})$ and  $\tilde{\Sigma}_{r,i} \triangleq ({A}_{cl}, {B}_{cl,i}, {C}_{r,i}, {D}_{r,i})$. Let us consider a linear time-invariant system of the form \eqref{eq1} with the attack vector $a$ as input $\mu(\cdot)$ and the vector $\begin{bmatrix}
y_{p,1}&y_{r,1},\dots,y_{pN_2}&y_{r,N_2}\end{bmatrix}^T$ as output $\sigma(\cdot)$. This system will be stable due to \textit{Assumption \ref{assume_stable}} and the system matrices of \eqref{eq1} would read
\begin{align}
\left[
\begin{array}{c|c}
\Phi & \Lambda\\
\hline
\\[\dimexpr-\normalbaselineskip+2pt] \Pi & \Upsilon
\end{array}
\right] =
\left[
\begin{array}{c|c}
\begin{matrix}
A_{cl,1}& 0 & 0 \\
\; & \ddots & \;\\
0 & 0 & A_{cl,N_2}
\end{matrix}     & \begin{matrix}
B_{cl,1} \\ \vdots \\ B_{cl,N_2}
\end{matrix} \\
\hline
\\[\dimexpr-\normalbaselineskip+2pt] \begin{matrix}
C_{p,1} & 0 & 0\\
C_{r,1} & 0 & 0\\
\; & \ddots & \; \\
0 & 0 & C_{p,N_2}\\
0 & 0 & C_{r,N_2}
\end{matrix} & \begin{matrix}
D_{p,1}\\
D_{r,1}\\
\vdots\\
D_{p,N_2}\\
D_{r,N_2}
\end{matrix}
\end{array}
\right],
\end{align}
For this system, let us define the set $\mathcal{L}$ as in \eqref{eq2} where $\omega = \begin{bmatrix} a^T& \sigma \end{bmatrix}^T \in \mathbb{R}^{n_a+N_2(n_p+n_r)}$. In view of \eqref{eq4}, let us also define the functionals $\mathcal{G}_0(\omega) \triangleq $
\begin{equation}\label{G_0}
\frac{1}{N_2} \sum_{i=1}^{N_2} ||y_{p,i}||_{\ell_2}^2 - \upsilon= \sum_{k=0}^{\infty} \omega[k]^TM_0\omega[k] + \zeta_0
\end{equation}
where $M_0 \in \mathbb{R}^{(n_a+N_2(n_p+n_r)) \times (n_a+N_2(n_p+n_r))}, \zeta_0=-\upsilon$ and $M_0(i,j)= 1, \text{if}\; i=j, n_a+i$  is odd and $0$ elsewhere. Similarly let us define
\begin{equation}\label{G_1}
    \mathcal{G}_k(\omega) \triangleq  -||y_{r,k}||_{\ell_2}^2+1 \;=\;\sum_{k=0}^{\infty} \omega[k]^TM_k\omega[k] + \zeta_k,
\end{equation}
$\forall k =\{1,2,\dots,N_2\}$, where $\zeta_k = -1$ and $M_k(i,j)= -1, \text{if}\; i=j, i=n_a+2k-1$ and $0$ elsewhere. Here, $M_k$ has same dimension as $M_0\;\forall k =\{1,2,\dots,N_2\}$. Therefore, we have shown that the constraints of \eqref{eq4} can be rewritten as functionals of the set $\mathcal{L}$. We can now see that the functionals $-\mathcal{G}_0(\cdot),\mathcal{G}_1(\cdot),\dots,\mathcal{G}_k(\cdot)$ along with \textit{Lemma \ref{vanish_lemma}} satisfy the conditions under which \textit{Theorem \ref{thm_new}} holds. Thus, by applying \textit{Theorem \ref{thm_new}}, if follow that the functionals $-\mathcal{G}_0(\cdot),\mathcal{G}_1(\cdot),\dots,\mathcal{G}_k(\cdot)$ form an S-system. Let this be argument 1.

In the case the adversary chooses not to attack the system, i.e., $a=0 \in \ell_{2e}$, it follows that $||y_{r,i}||_{\ell_2}^2 \approx 0 \;\forall \delta_i$. The residual energy $||y_{r,i}||_{\ell_2}^2$ is not strictly zero since there might be residual outputs due to difference in initial condition between the system and the detector. The threshold ($\epsilon_r = 1$) is chosen in such a way that $||y_{r,i}||_{\ell_2}^2 \ll 1\; \forall \delta_i$ when $a=0$. Thus, it holds that $\exists \; \omega_0 = [a,\;\omega] = [0,\; 0_{+}] $ s.t. $ -||y_{r,k}||_{\ell_2}^2+1 = \mathcal{G}_k(\omega_0) >0\; \forall k=\{1,\dots,N_2\}$. Here $0_{+}$ represents a real number close to zero. Let this be argument 2. And from \eqref{eq4}, we know that the system $\mathcal{G}_0(\omega) > 0 , \mathcal{G}_i(\omega) \geq 0 \;i=\{1,\dots,k\}$ is not solvable. Let this be argument 3. Using the above arguments (1-3) and \cite[Theorem 4.3.1]{petersen2012robust}, we can conclude that when we consider the functionals defined in \eqref{G_0} and \eqref{G_1}, $\exists \; \gamma_1 \geq 0,\gamma_1 \geq 0,\dots,\gamma_{N_2}\geq 0$ such that the following inequality holds
\begin{equation}\label{eq5}
 \mathcal{G}_0(\omega) + \sum_{i=1}^{N_2}\gamma_i\mathcal{G}_i(\omega) \leq 0,\; \forall \; \omega \in \mathcal{L}.
\end{equation}
To conclude, in this step we have shown that the constraint of \eqref{eq4} holds only if \eqref{eq5} is true. And the converse is generally true.
\subsubsection*{Step 3}
We have shown that the constraint of \eqref{eq4} holds iff \eqref{eq5} is true. Then, we reformulate \eqref{eq4} as
\begin{equation}
\inf_{\upsilon,\gamma_1\geq 0,\dots,\gamma_{N_2}\geq 0} \left\{ \upsilon \;\Big|\; \mathcal{G}_0(\omega) + \sum_{i=1}^{N_2}\gamma_i\mathcal{G}_i(\omega) \leq 0,\; \forall \; \omega(\cdot) \in \mathcal{L}\right\}.
\end{equation}
Substituting the definition of $\mathcal{G}_0(\omega)$, we obtain
\begin{equation}\label{eq6}
\inf_{\gamma\geq 0}\left\{ \inf_{\upsilon} \left\{\upsilon\Big|\sum_{i=1}^{N_2} \left\{ \frac{1}{N_2}||y_{p,i}||_{\ell_2}^2\; + \gamma_i\mathcal{G}_i(\omega)\right\} \leq \upsilon,\forall\omega \right\}\right\}
\end{equation}
where $\gamma=[\gamma_1,\dots,\gamma_{N_2}]^T$. The inner optimization problem of \eqref{eq6} resembles an epigraph formulation which can be rewritten as
\begin{equation}
\inf_{\gamma\geq 0}\left\{
\begin{aligned}
\sup_{\omega} \; & \left\{\sum_{i=1}^{N_2} \left\{ \frac{1}{N_2}||y_{p,i}||_{\ell_2}^2\;
 + \gamma_i\mathcal{G}_i(\omega)\right\}\right\}
\end{aligned}\right\}.
\end{equation}
Substituting the definition of $\mathcal{G}_i(\omega),\;\forall i=\{1,\dots,N_2\}$, we obtain
\begin{equation}\label{eq7}
\inf_{\gamma\geq 0}\underbrace{
\sup_{\omega} \left\{\sum_{i=1}^{N_2} \left\{ \frac{1}{N_2}||y_{p,i}||_{\ell_2}^2-\gamma_i||y_{r,i}||_{\ell_2}^2\right\}+\sum_{i=1}^{N_2} \gamma_i \right\}
}_{\kappa}
\end{equation}
Observe that $\kappa$ is a maximization problem with a quadratic term in its objective. Thus, it holds that
\[
\kappa =
\begin{cases}
    \mathbf{1}^T \gamma, & \text{if}\;\;\sum_{i=1}^{N_2} \left\{ \frac{1}{N_2}||y_{p,i}||_{\ell_2}^2\;
 - \gamma_i||y_{r,i}||_{\ell_2}^2\right\} \leq 0\\
    +\infty,  & \text{otherwise}
\end{cases}.
\]
Using the above result in \eqref{eq7}, we obtain
\begin{equation*}
\inf_{\gamma \geq 0}\; \left\{ \mathbf{1}^T \gamma \Bigg| \sum_{i=1}^{N_2} \frac{1}{N_2}||y_{p,i}||_{\ell_2}^2 - \gamma_i ||y_{r,i}||_{\ell_2}^2 \leq 0, \forall a \in \ell_{2e}\right\}
\end{equation*}
where $\omega$ is replaced by $a$ since it is only the control variable. Finally, we add the state constraints that were removed while formulating the optimization problem \eqref{eqq2}
\begin{equation}\label{eqq3}
\inf_{\gamma \geq 0} \left\{ \mathbf{1}^T \gamma \;\Bigg| \begin{aligned}
&\sum_{i=1}^{N_2} \frac{1}{N_2}||y_{p,i}||_{\ell_2}^2 - \gamma_i ||y_{r,i}||_{\ell_2}^2 \leq 0, \forall a \in \ell_{2e}\\
&x_{i}[\infty]=0 \;\forall i \in \{1,\dots,N_2\} \end{aligned}\right\}
\end{equation}
Thus in this step, we have shown using S-procedure that the optimization problem \eqref{eqq1} and \eqref{eqq3} are equivalent.
\subsubsection*{Step 4}
Define $\bar{x}[\infty] =$
$ \begin{bmatrix} x_{1}[\infty]^T & \dots & x_{N_2}[\infty]^T \end{bmatrix}^T, \bar{x}[0] = \begin{bmatrix} x_{1}[0]^T & \dots & x_{N_2}[0]^T \end{bmatrix}^T$, \\
$\bar{y}_p = \begin{bmatrix} y_{p1}^T & \dots & y_{pN_2}^T
\end{bmatrix}^T$, and $\bar{y}_r = \begin{bmatrix} y_{r1}^T & \dots & y_{rN_2}^T \end{bmatrix}^T$. Using these definitions, the constraint of \eqref{eqq3} can be rewritten as
\begin{align}\label{con1}
-\frac{1}{N_2}||\bar{y}_p||_{\ell_2}^2 + ||\sqrt{\Gamma(\gamma)}\bar{y}_r||_{\ell_2}^2 \geq 0,  \forall a \in \ell_{2e},\;\; \bar{x}[\infty]=0
\end{align}
where $\Gamma(\gamma)$ is defined in the theorem statement. Additionally due to \textit{Assumption \ref{equil_begin}}, we have $\bar{x}[0] =0$. Next let us define $y_1 = \sqrt{\Gamma(\gamma)}\bar{y}_r, y_2=\frac{1}{\sqrt{N_2}}\bar{y}_p$ and the supply rate $s(\cdot) \triangleq  ||{y_1}||_{2}^2 -||{y_2}||_{2}^2$. Then, we have shown from \eqref{con1} that \cite[\textit{Proposition 2, 2)}]{truncated} holds. Equivalently, using \cite[\textit{Proposition 2, 3)}]{truncated}, we replace \eqref{con1} by the constraint of \eqref{problem:coupled:final} where $\bar{\Sigma}_p = (\bar{A},\bar{B},\bar{C}_p,\bar{D}_p)$ and $\bar{\Sigma}_r = (\bar{A},\bar{B},\bar{C}_r,\bar{D}_r)$ represent the system with the attack as input and $\bar{y}_p$ and $\bar{y}_r$ as system outputs respectively. Constructing these system matrices, as we did in \textit{Step 2} of this proof concludes the proof.
\end{proof}

\section{Proof of \textit{Lemma \ref{Bound_RA}}}\label{bound_2}

\begin{proof}
To recall, the optimization problem \eqref{problem:coupled:final} was formulated using \cite[\textit{Proposition 2, 3)}]{truncated} where $y_1=\sqrt{\Gamma(\gamma)} \bar{y}_r$ and $y_2= \frac{1}{\sqrt{N_2}}\bar{y}_p$. Here $\bar{y}_p$ and $\bar{y}_r$ represents the outputs of $\bar{\Sigma}_p$ and $\bar{\Sigma}_r$ respectively. Due to the equivalency between $3)$ and $4)$ of \cite[\textit{Proposition 2, 3)}]{truncated}, the FDI \cite[\textit{Proposition 2, 4)}]{truncated} should hold $\;\forall \;|z|=1$. Since we know that $y_1 = \sqrt{\Gamma(\gamma)} \bar{y}_r$ and $y_2 = \bar{y}_p$, we can deduce that $G_1(z)$ corresponds to $\sqrt{\Gamma(\gamma)}\bar{G}_{r}(z)$ and $G_2(z)$ to $\frac{1}{N_2}\bar{G}_{r}(z)$ in \cite[\textit{Proposition 2, 4)}]{truncated}, where $\bar{G}_{r}(z)= \bar{C}_{r}(z_1I-\bar{A})^{-1}\bar{B}+\bar{D}_{r}$ and $\bar{G}_{p}(z) \triangleq \bar{C}_{p}(z_1I-\bar{A})^{-1}\bar{B}+\bar{D}_{p}$. Thus, \eqref{problem:coupled:final} can be rewritten as
\begin{equation}\label{lem_s2}
\inf\Big\{ \textbf{1}^T\gamma\Big|\bar{G}_{r}(\bar{z})^T\Gamma(\gamma)\bar{G}_{r}({z})-\bar{G}_{p}^T(\bar{z})\bar{G}_{p}({z})\succeq 0,\forall|z|=1\Big\}
\end{equation}
Let us define the following sets such that $\mathbb{C}^{n_a} = \mathcal{Z}_{pr} \cup \mathcal{Z} \cup \mathcal{Z}_{r} \cup \mathcal{Z}_{p}$.
\begin{align}
    \mathcal{Z}_{pr} & \triangleq \{x \in \mathbb{C}^{n_a} \;|\; \bar{G}_{r}({z})x = 0 , \bar{G}_{p}({z})x = 0\},\\
    \mathcal{Z} & \triangleq \{x \in \mathbb{C}^{n_a} \;|\; \bar{G}_{r}({z})x \neq 0 , \bar{G}_{p}({z})x \neq 0\},\\
    \mathcal{Z}_r & \triangleq \{x \in \mathbb{C}^{n_a} \;|\; \bar{G}_{r}({z})x = 0 , \bar{G}_{p}({z})x \neq 0\},\\
    \mathcal{Z}_p & \triangleq \{x \in \mathbb{C}^{n_a} \;|\; \bar{G}_{r}({z})x \neq 0 , \bar{G}_{p}({z})x = 0\}.
\end{align}
\textit{Sufficiency:} For any given $z$ such that $|z|=1$ , if $ x \in \mathcal{Z}_p$ or $ x \in \mathcal{Z}_{pr}$, choosing $\Gamma(\gamma) =0$ satisfies the constraint of \eqref{lem_s2}. Similarly, if $ x \in \mathcal{Z}$, let us pick $\Gamma(\gamma) = cI_{n_{r}}$ where c is a constant. Then, the value of \eqref{lem_s2} is bounded if there exists a bounded $c$ that makes $ \Xi \triangleq \sup_{|z|=1, x \in \mathcal{Z}} \frac{x^H\big\{\bar{G}_{r}^T(\bar{z})\bar{G}_{r}({z})\big\}x}{x^H\big\{\bar{G}_{p}^T(\bar{z})\Gamma(\gamma)\bar{G}_{p,i}({z})\big\}x}$ bounded. $\Xi$ is bounded  since the denominator cannot become zero (since $x \in \mathcal{Z}$ and $\Gamma(\gamma)$ is full rank), and we have assumed that the $\bar{G}_{r}({z})$ and $\bar{G}_{p}({z})$ are stable (\textit{Assumption \ref{assume_stable}}). Next we prove sufficiency when $x \in \mathcal{Z}_r$.

When condition $1)$ of the lemma statement holds, by definition of a zero $ \forall |z|=1, \nexists s \neq 0 \in \mathbb{C}^{n_a}$ such that $\bar{G}_{r}({z})s = 0$. Thus it follows that $\mathcal{Z}_r = \mathcal{Z}_{pr} = \emptyset$. When condition $2)$ of the lemma statement holds, by definition of a zero $ \forall |z|=1, \nexists s \neq 0$ such that $\bar{G}_{r}({z})s = 0$ and $\bar{G}_{p}({z})s \neq 0$. Thus it follows that $\mathcal{Z}_r = \emptyset$. \\
\textit{Necessity:} Assume that there exists a bounded $\Gamma(\gamma)$ that solves \eqref{lem_s2}. We also assume that there exists a complex number $z_1$ on the unit circle which is a zero of the system $\bar{\Sigma}_{r}$ (including multiplicity and input direction) but are not zeros of $\bar{\Sigma}_{p}$. By definition of a zero, it holds that $\exists s \neq 0$ such that $\bar{G}_{r}(z_1)s=0$ and $\bar{G}_{p}(z_1)s \neq 0$. Under these assumptions, when $z=z_1$ and $x = s$, the constraint of \eqref{lem_s2} can be rewritten as $-s^H\bar{G}_{p}^T(\bar{z}_1)\bar{G}_{p,i}({z_1})s \geq 0$ which cannot hold since $\bar{G}_{p}({z_1})s \neq 0$. That is, the feasibility set of \eqref{lem_s2} is empty which contradicts our assumption. This concludes the proof.
\end{proof}
\section{Introduction}\label{Intro}

\input{1_Introduction}

\section{Problem background}\label{description}

\input{2A_Description}

\section{Problem formulation}\label{PF}

\input{2B_Problem_Formulation}

\section{Risk assessment for a bounded-rational adversary}\label{rational}

\input{4_Rational}

\section{Numerical example}\label{Example}

\input{5_Example}

\section{Conclusion and Future work}\label{Conclusion}

\input{6_Conclusion}

\bibliographystyle{ieeetr}


\section{Introduction}
\label{sec:intro}

Hyperspectral (HS) cameras sense the electromagnetic spectrum to produce images
that depict a scene across several contiguous bands. Such technology has led to
improved performance in various tasks, including anomaly detection~\cite{Xu2016Anomaly} and remote sensing~\cite{Ojha2015spectral}.

HS cameras, however, suffer from low spatial and temporal resolution
\cite{Schechner2002generalized}. RGB cameras, on the other hand, acquire images
with high spatial resolution but low spectral resolution, since they split the
spectrum into three broad bands. This tradeoff between spatial and spectral
resolution stems from the physics of the acquisition process. To overcome it,
one can fuse a low (spatial) resolution hyperspectral (LRHS) image with a high resolution (spatial) RGB/multispectral (HRMS) image to obtain the high resolution hyperspectral (HRHS) image, a method known as fusion based HS image super-resolution.

\mypar{Acquisition model}
Let $\boldsymbol{X} \in \mathbb{R}^{M_0 \cdot N_0 \times S_0}$ denote the unknown HRHS that we wish to reconstruct. Each column represents a vectorized HRHS of size $M_0 \times N_0$ on a given spectral band, and there are $S_0$ spectral bands. We assume access to a low-resolution HS image acquired by a HS camera, $\boldsymbol{Z} \in \mathbb{R}^{M\cdot N \times S_0}$, where $M < M_0$ and $N < N_0$, and to a high-resolution RGB image acquired by a conventional camera, $\boldsymbol{Y} \in \mathbb{R}^{M_0 \cdot N_0 \times S}$, where $S<S_0$. Typically, $S = 3$ and $S_0$ is greater than 20. All these quantities are related as

\begin{equation}
\boldsymbol{Z} = \boldsymbol{AX}, \quad \boldsymbol{Y}=\boldsymbol{XR},
\label{eq:cond1}
\end{equation}
where $\boldsymbol{A} \in \mathbb{R}^{M \cdot N \times M_0 \cdot N_0}$ is the downscaling operator, and $\boldsymbol{R} \in \mathbb{R}^{S_0 \times S}$ is the camera spectral response (CSR) function of the RGB camera i.e., it integrates all bands of the spectrum into the RGB image.

\begin{figure}[t]
	\includegraphics[width=8.5cm, height=4.3cm]{figs/frame_new.png}
	\caption{Our framework. A high-resolution RGB image $\boldsymbol{Y}$ (HRMS) and a low-resolution HS image $\boldsymbol{Z}$ (LRHS) are first super-resolved into $\boldsymbol{W}$ with a fusion-based hyperspectral imaging method (typically, a convolutional neural network). We then create a measurement-consistent high-resolution HS image $\boldsymbol{\widehat{X}}$ (HRHS) by solving TV-TV minimization using $\boldsymbol{Z}$, $\boldsymbol{Y}$, and $\boldsymbol{W}$ as inputs.}
	\label{fig:frameworkV2}
	\vspace{-0.4cm}
\end{figure}

\mypar{Overview of existing approaches} Reconstructing $\boldsymbol{X}$ from $\boldsymbol{Z}$ and $\boldsymbol{Y}$ in (\ref{eq:cond1}) is an ill-posed problem, as there are generally more unknowns than equations. This task thus requires assumptions about the structure of $\boldsymbol{X}$ (prior knowledge). Examples include
sparsity \cite{Akhtar2014sparsespatiospectral} and spectral correlation \cite{nezhad2016fusion}. Other approaches to recover $\boldsymbol{X}$ encode prior knowledge about the HS observations via matrix \cite{Kawakami2011tfm, Yokoya2012cnmf} and tensor \cite{dian2017nlstf} factorization. We refer to this class of methods, in which prior knowledge is explicitly encoded, as \textit{model-based} methods.

Deep learning (DL) methods have also been applied to reconstruct $\boldsymbol{X}$ from $\boldsymbol{Z}$ and $\boldsymbol{Y}$, e.g., \cite{xie2019fusion, qu2018unsupervised}. Although they require extensive training datasets, they bypass the need to encode assumptions explicitly and are computationally efficient once trained. Currently, they are the state-of-the-art for this task. DL methods, however, suffer from two important drawbacks: they generally lack interpretability, and during testing, they are unable to enforce measurement consistency between the inputs ($\boldsymbol{Z}, \boldsymbol{Y}$) and their output, i.e., enforce (\ref{eq:cond1}) after they have been trained. Conventional model-based algorithms do not suffer from these drawbacks, but produce outputs of lower quality.

\mypar{Problem statement} Our aim is then to design a method for fusion based HS image super-resolution that not only produces high quality outputs, but also guarantees measurement consistency for every test image. By guaranteeing such consistency, we aim to preserve fine details in the measurements, which can be critical in applications like airport security and space surveillance.

We propose to combine DL (or any other RGB-guided HS image super-resolution method) with model-based methods. Fig. \ref{fig:frameworkV2} illustrates our scheme, described in Section \ref{ourframework}. In brief, the inputs $\boldsymbol{Y}$ and $\boldsymbol{Z}$ are passed through a fusion based super-resolution method to obtain $\boldsymbol{W}$. We then solve an
optimization problem termed as TV-TV minimization that takes $\boldsymbol{W}, \boldsymbol{Z}$ and $\boldsymbol{Y}$ as inputs to generate an improved version of $\boldsymbol{W}$ denoted by $\boldsymbol{\widehat{X}}$. This approach extends our original work \cite{vella2019tvtvminimization, Vella20-Robust}, where we considered single image super-resolution. In this paper, we generalize our original framework to fusion based HS images by considering the additional constraint $\boldsymbol{Y} = \boldsymbol{XR}$, which relates data from different modalities. Our extensive experiments show that the proposed method surpasses the current state-of-the-art for fusion HS super-resolution in terms of several reconstruction metrics. To the best of our knowledge, this is the first work that tackles the data consistency problem characterizing DL networks in the field of fusion based HS imaging super-resolution.

\section{Related Work}

In recent years, single-image super-resolution methods that process conventional RGB images have been extended to HS imaging. Traditional HS image super-resolution methods take a LRHS image as an input and output a HRHS image, without the aid of any other modality. Because of this, such methods have two shortcomings: the scaling factor is usually limited to four, and they increase the spectral resolution only. To overcome these shortcomings, several techniques fuse a LRHS image with a conventional HRMS image. Most methods are model-based and require explicit assumptions about the HRHS image, but more recent algorithms build on DL networks. The latter is still an emerging field with limited literature. We give a brief overview of both techniques since we utilize both in our work.

\mypar{Model-based methods} Classic model-based methods recover $\boldsymbol{X}$ by making explicit assumptions about it. For example, in \cite{Simoes2015tv}, it is assumed that $\boldsymbol{X}$ has a small TV norm. The work in \cite{Akhtar2014sparsespatiospectral} assumes that the spatial information of the HRHS image has a sparse representation that can be learned from a dictionary.

A different line of work is based on matrix factorization methods \cite{Kawakami2011tfm, Yokoya2012cnmf}. These assume that the HS observations can be factorized into two matrices, one representing the spectral response of the material, and the other containing the proportion of materials at each pixel location. For example, \cite{Yokoya2012cnmf} proposed an algorithm that uses coupled non-negative matrix factorization (CNMF) and assumes no sparsity constraints. Even though matrix factorization approaches are popular, they ignore the spectral-spatial correlations. Motivated by this, more recent methods use tensor factorization. For example, \cite{dian2017nlstf} uses non-local sparse tensor factorization to decompose each cube of the HS image as a sparse tensor and dictionaries of three different modes. In addition, similar cubes are processed together to leverage the non-local self-similarities of the HS image. Such assumptions are common in other HS-related tasks, for example, tensor denoising \cite{gong2020}.

\mypar{Deep learning methods} DL has been widely applied to computer vision and image processing, e.g., \cite{MSLapSRN, Kim2016VDSR}, but its application to HS image super-resolution is still relatively recent \cite{Palsson2017fusion,Han2018ssfcnn, xie2019fusion, Sidrov2019Deephsi, Zhang2019hsir}. DL methods are computationally efficient during testing and offer superior performance. In reconstruction tasks, trained DL networks, however, cannot ensure that the output is consistent with the input \cite{Vella20-overcoming}. This is because the weights learned during training are applied to all test images, which makes it difficult for such networks to ensure that the observation model holds in the testing stage. The work in \cite{xie2019fusion} attempts to solve this problem by unfolding an iterative algorithm into a deep network that takes into consideration the low rankness and sparsity of the HRHS image. Also, \cite{qu2018unsupervised} proposed an unsupervised method consisting of two networks with an encoder-decoder structure used for representation learning. The first network takes the HRMS image as an input while the other network processes the LRHS image. To exploit spectral and spatial information, the learned weights are then shared before reconstructing the HRHS image. The proposed network also forces the representations of both inputs to follow a sparse Dirichlet distribution. Unsupervised networks solve a fundamentally harder task, thus their performance is naturally worse than supervised methods \cite{l2019handsonul}.

 \begin{table}[t]
   \centering
   \renewcommand{\arraystretch}{1}
   \caption{Average PSNR, SSIM, SAM, ERGAS and RMSE on all HS images from the Harvard dataset.}

   \begin{tikzpicture}
     \filldraw[color=black!10] (-2.4,-1.71) rectangle +(2.9,3.4);
     \filldraw[color=black!10] (0.73,-1.71) rectangle +(3.0,3.4);
     \draw node {
         \begin{tabular}{@{}llllll@{}}\toprule
           \textcolor{thisblue}{\textbf{Metric}} &
           \textcolor{darkgreen}{\sf \textit{CNMF}} \cite{Yokoya2012cnmf}  & 		\textcolor{darkredc} {\sf \textit{\textbf{Ours}}}  & 			\textcolor{darkgreen}{\sf \textit{uSDN}} \cite{qu2018unsupervised}  &
           \textcolor{darkredc} {\sf \textit{\textbf{Ours}}} \\
           \midrule
           {{PSNR} } &41.381 & \textbf{41.574}  & 38.564 & \textbf{41.120}\\
           \midrule
           {SSIM} & \phantom{0}0.985 & \phantom{0}\textbf{0.986} & \phantom{0}{0.990}&  \phantom{0}\textbf{0.993}\\
           \midrule
           {SAM} & \phantom{0}3.951 &\phantom{0}\textbf{3.831} & \phantom{0}4.503&  \phantom{0}\textbf{3.441} \\
           \midrule
           {ERGAS} & \phantom{0}0.314 &\phantom{0}\textbf{0.311} & \phantom{0}0.563 & \phantom{0}\textbf{0.431}\\
           \midrule
           {RMSE} & \phantom{0}2.456&  \phantom{0}\textbf{2.407}& \phantom{0}2.753 & \phantom{0}\textbf{1.817}\\
           \bottomrule
         \end{tabular}
       };
   \end{tikzpicture}

   \label{tab:harvard}
   \vspace{-0.4cm}
 \end{table}

  \begin{table*}[t]
	\sffamily\selectfont
	\centering
	\renewcommand{\arraystretch}{1}
	\caption{Average PSNR, SSIM, SAM, ERGAS and RMSE on 12 test HS images used by MHF from the CAVE dataset.}
    \begin{tikzpicture}
    \filldraw[color=black!10] (-5.8,-1.71) rectangle +(3.05,3.4);
     \filldraw[color=black!10] (-2.5,-1.71) rectangle +(3.05,3.4);
     \filldraw[color=black!10] (0.79,-1.71) rectangle +(3.05,3.4);
     \filldraw[color=black!10] (4.1,-1.71) rectangle +(3.15,3.4);
     \draw node {
		\begin{tabular}{@{}lllllllll@{}}

			\toprule[1.1pt]
			\textcolor{thisblue}{\sf \textbf{Metric}} &
			\textcolor{darkgreen}{\sf \textit{CNMF}} \cite{Yokoya2012cnmf} &
			\textcolor{darkredc} {\sf \textit{\textbf{Ours}}} &
			\textcolor{darkgreen}{\sf \textit{NLSTF}} \cite{dian2017nlstf} &
			\textcolor{darkredc}{\sf \textit{\textbf{Ours}}}   &
			\textcolor{darkgreen}{\sf \text{uSDN}} \cite{qu2018unsupervised}  &
			\textcolor{darkredc} {\sf \textit{\textbf{Ours}}} &
			\textcolor{darkgreen}{\sf \textit{MHF-net}} \cite{xie2019fusion}   &
			\textcolor{darkredc} {\sf \textit{\textbf{Ours}}}
			\\
			\midrule
			{PSNR} & 37.149 & \textbf{ 37.881} & 40.859 & \textbf{40.942} & 36.365 & \textbf{39.712 }& 37.507 & \textbf{38.151}\\
			\midrule
			{SSIM} &  $\phantom{0}$0.982 & $\phantom{0}$ \textbf{0.984 }&$\phantom{0}$ 0.991 &$\phantom{0}$\textbf{0.991}  &$\phantom{0}$0.979 &$\phantom{0}$\textbf{0.987} & $\phantom{0}$0.977 &$\phantom{0}$\textbf{0.982} \\
			\midrule
			{SAM} &$\phantom{0}$7.482 &$\phantom{0}$ \textbf{7.263}  &$\phantom{0}$ 4.723 &$\phantom{0}$\textbf{4.719}& $\phantom{0}$ 6.823&$\phantom{0}$\textbf{6.051} &$\phantom{0}$8.221  &$\phantom{0}$\textbf{7.946}\\
			\midrule
			{ERGAS} &$\phantom{0}$0.629 &$\phantom{0}$ \textbf{0.568}  & $\phantom{0}$ 0.384  &$\phantom{0}$\textbf{0.379}  & $\phantom{0}$ 0.716 & $\phantom{0}$\textbf{0.533}  &$\phantom{0}$0.590 & $\phantom{0}$\textbf{0.541}\\
			\midrule
			{RMSE} &$\phantom{0}$3.805 &$\phantom{0}$ \textbf{3.487}&$\phantom{0}$ 2.598  &$\phantom{0}$\textbf{2.568}   & $\phantom{0}$ 4.235 &$\phantom{0}$\textbf{2.877} & $\phantom{0}$3.845 &$\phantom{0}$\textbf{3.565}\\
			\bottomrule
		\end{tabular}
	};
	\end{tikzpicture}
	\label{tab:cave}
	\vspace{-0.4cm}
\end{table*}

\section{Proposed Framework}
\label{ourframework}

Our goal is to recover the HRHS image $\boldsymbol{X} \in \mathbb{R}^{M_0.N_0 \times S_0}$ from a LRHS image $\boldsymbol{Z} \in \mathbb{R}^{M.N \times S_0}$ and an RGB image $\boldsymbol{Y} \in \mathbb{R}^{M_0.N_0 \times S}$. To achieve this, we use the scheme depicted in Fig. \ref{fig:frameworkV2}. First, the inputs $\boldsymbol{Z}$ and $\boldsymbol{Y}$ are super-resolved into $\boldsymbol{W} \in \mathbb{R}^{M_0.N_0 \times S_0}$ via a fusion based HS image super-resolution method. We use a combination of both model and DL based methods and then post-process $\boldsymbol{W}$ via an additional block, TV-TV minimization, which enforces consistency for both the spatial and spectral measurements.

\mypar{TV-TV minimization} The images $\boldsymbol{W}$, $\boldsymbol{Y}$ and $\boldsymbol{Z}$ are then simultaneously processed via TV-TV minimization, which generates an estimate $\boldsymbol{\widehat{X}}$ that satisfies the model in (\ref{eq:cond1}), is not too different from $\boldsymbol{W}$, and has a small total variation (TV) norm. More specifically, we solve

 \begin{equation}
\begin{array}{cl}{\underset{\boldsymbol{X} \in \mathbb{R}^{M_0.N_0 \times S_0}}{\operatorname{minimize}}} & {\|\boldsymbol{X}\|_{\mathrm{\text{TV}}}+\beta\|\boldsymbol{X} - \boldsymbol{W}\|_{\text{TV}}} \\  {\text {subject to}} & {\boldsymbol{AX}=\boldsymbol{Z}}\\ & {\boldsymbol{XR}=\boldsymbol{Y},}\end{array}
\label{eq:tv-tv}
\end{equation}
where $\boldsymbol{X} \in \mathbb{R}^{M_0.N_0 \times S_0}$ is the optimization variable and $\beta \geq 0$ is a trade-off parameter. We chose the TV norm since it is a widely used prior for image processing tasks. Our framework can be easily adapted to accommodate other priors. The TV-norm $\| \boldsymbol{X} \|_{\text{TV}}$ in (\ref{eq:tv-tv}) for 3D tensors is defined as the sum of the 2D TV-norms of each spectral band i.e., $\| \boldsymbol{X} \|_{\text{TV}} = \sum_{i=1}^{S_0}\left\|X_{s}\right\|_{\mathrm{TV}}$ where $X_s \in \mathbb{R}^{M_0 \times N_0}$ represents the band $s$ of $X$. The 2D TV norm of a vectorized image $x \in \mathbb{R}^{M_0. N_0}$ is defined as
$$
\left\|x\right\|_{\mathrm{TV}}:=\sum_{i=1}^{M} \sum_{j=1}^{N}\big|D_{i j}^{v} x\big|+\big|D_{i j}^{h} x\big| = \|\boldsymbol{D}x\|_1,
$$
where $D^v_{ij}$ (resp. $D^h_{ij}$) is a row-vector that extracts the vertical (resp. horizontal) difference at pixel $(i, j)$ of $\boldsymbol{X}$, and $\boldsymbol{D}$ is the vertical concatenation of all the $D^v_{ij}$ and $D_{ij}^h$.

To solve (\ref{eq:tv-tv}), we use the alternating direction method of multipliers (ADMM) \cite{boyd_distributed_2011}. First we introduce a set of auxiliary variables, $u_s$ and $v_s$, which we constrain as  $\boldsymbol{D}v_s = u_s$ and $x_s=v_s$ where $u_s$ (resp. $v_s$ and $x_s$) represents the vectorized $s$th spectral band of $u$ (resp. $v$ and $x$), as suggested in \cite{liu_depth_2015}, and dualize both constraints. Specifically, we rewrite (\ref{eq:tv-tv}) as
\begin{equation}
  \begin{array}[t]{ll}
    \underset{(u,\boldsymbol{X}),v}{\text{minimize}} &  \sum^{S_0}_{s=1} \big( \|u_s\|_1 + \beta \|u_s-\overline{w}_s\|_1 + \\ &
    \text{i}_{\{x_s:Ax_s=z_s\}}(x_s) \big)+ \text{i}_{\{\boldsymbol{X}:\boldsymbol{XR}=\boldsymbol{Y}\}}(\boldsymbol{X})\\ \\
        \text{subject to} & \boldsymbol{D}v_s=u_s, s = 1, \cdots , S_0 \\
        & v_s=x_s, s = 1, \cdots , S_0 \, ,
  \end{array}
  \label{eq:tvvec}
\end{equation}
where $\overline{w}_s := \boldsymbol{D}w_s \in \mathbb{R}^{M_0.N_0}$ and $\text{i}_F(u)$ is
an indicator function, i.e., $\text{i}_F(u) = 0$ if $u \in F$, and
$\text{i}_F(u) = +\infty$ otherwise. Once we obtain a solution ($u^\star, X^\star, v^\star$) of (\ref{eq:tvvec}), a solution of (\ref{eq:tv-tv}) is given by $X^\star = \begin{bmatrix} x^{\star}_1 & \cdots & x^{\star}_{S_0}\end{bmatrix}$.

\mypar{ADMM iterates} We use two dual variables $\lambda \in \mathbb{R}^{2 \cdot M \cdot N}$ and $\mu \in \mathbb{R}^{M \cdot N} $ and solve (\ref{eq:tvvec}) via the following ADMM iterates:
\begin{align}
  \boldsymbol{U}^{k+1}
  &=
  \underset{u_1 \cdots u_{S_0}}{\arg\min}\,\,\,
  \sum_{s=1}^{S_0} \|u_s\|_1 + \beta \|u_s-w_s\|_1 +
  \notag
  \\
  &\hspace{3.1cm}
  \big(\lambda_s^{k} - \rho  \boldsymbol{D}v_s^k\big)^T u_s
  + \frac{\rho_s}{2} \|u_s\|^2_2
  \label{Eq:ADMMSubProbU}
  \\
  \boldsymbol{X}^{k+1}
  &=
  \underset{\boldsymbol{X}}{\arg\min}\,\,\,
  \mu^{k^T}\boldsymbol{X} + \frac{\rho}{2}
  \|\boldsymbol{X}-\boldsymbol{V}^k\|^2_2
  +
  \text{i}_{\{\boldsymbol{X}:\boldsymbol{AX}=\boldsymbol{Z}\}}(\boldsymbol{X})
  \notag
  \\
  &\hspace{4.3cm}
  +
  \text{i}_{\{\boldsymbol{X}:\boldsymbol{XR}=\boldsymbol{Y}\}}(\boldsymbol{X})
  \label{Eq:ADMMSubProbX}
  \\
  \boldsymbol{V}^{k+1}
  &=
  \underset{v_1 \cdots v_{S_0}}{\arg\min}\,\,\,
  \sum_{s=1}^{S_0} -\lambda_s^{k^T} \boldsymbol{D}v_s - \mu^{k_s^T}v_s
  + \frac{\rho}{2} \|x_s^{k+1}-v_s\|^2_2
  \notag
  \\
  &\hspace{1.7cm}
  + \frac{\rho}{2} \|x_s^{k+1}-v_s\|^2_2 + \frac{\rho}{2} \|u_s^{k+1}- \boldsymbol{D}v_s\|^2_2
  \label{Eq:ADMMSubProbV}
  \\
  \lambda^{k+1}
  &=
  \sum^{S_0}_{s=1} \lambda_s^k + \rho(u_s^{k+1}- \boldsymbol{D}v_s^{k+1})
  \\
  \mu^{k+1}
  &=
  \sum_{s=1}^{S_0} \mu_s^k + \rho(x_s^{k+1}- v^{k+1}_s)\,,
  \label{Eq:lmult}
\end{align}
where $\rho > 0$ is the augmented Lagrangian parameter. Recall that the subscript $s$ represents the $s$th band of each of the hyperspectral image. Then the concatenation of each band gives the corresponding hyperspectral cube. Problem (\ref{Eq:ADMMSubProbU}) is solved by equating the element-wise derivative with respect to $u_s$ to zero, (\ref{Eq:ADMMSubProbX}) is the projection of a point onto a linear subspace, and (\ref{Eq:ADMMSubProbV}) can be solved by equating the first order derivative with respect to $v_s$ to zero. Thus, each subproblem has a closed form solution.
Note that problems (\ref{Eq:ADMMSubProbU}) and (\ref{Eq:ADMMSubProbV}) can each be decomposed into $S_0$ independent problems that can be solved in parallel, which considerably speeds up the ADMM iterates.

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

We now describe our testing procedure and present visual and quantitative results. Code to replicate our experiments is available online\footnote{https://github.com/marijavella/hs-sr-tvtv}.

\mypar{Datasets} In this paper we run experiments on two popular HS imaging datasets, CAVE \cite{yasuma2010cave} and Harvard \cite{chakrabarti2011statistics}. The CAVE dataset consists of 32 HRHS images, each with dimensions $512 \times 512$ and 31 spectral bands. The spectral images are taken with a wavelength ranging between $400 \sim 700$nm at an interval of 10nm. The Harvard dataset consists of 50 HRHS images with dimensions $1392 \times 1040$ and 31 spectral bands each. Each spectral image is taken at 10nm intervals and the wavelength ranges between $420 \sim 720$nm. For this dataset, the top left 1024 $\times$ 1024 pixels were considered.

\mypar{Base methods} We use and compare against four state-of-the-art base methods: CNMF \cite{Yokoya2012cnmf}, NLSTF \cite{dian2017nlstf}, MHF \cite{xie2019fusion} and uSDN \cite{qu2018unsupervised}. The first and second methods are based on matrix and tensor factorization, respectively, while the last two use DL networks. The outputs from each base method were obtained using the default settings and parameters provided by the respective authors. We assume noiseless images and thus remove any added noise in any of the comparison methods.

\mypar{Experimental Settings} All images were first normalized between 0 and 1. For the CAVE dataset, we used the same 12 HS images/cubes tested with MHF while for the Harvard dataset, we use all the 50 available HS images/cubes. The HSLR images $\boldsymbol{Z}$ were obtained according to (\ref{eq:cond1}) with $\boldsymbol{A}$ representing the averaging operation over non-overlapping blocks of size $32 \times 32$, that is, as a downscaling operation with a factor of 32. The HRMS images $\boldsymbol{Y}$ were generated according to (\ref{eq:cond1}), with $\boldsymbol{R}$ representing either the camera spectral response of Nikon D700 \cite{qu2018unsupervised, dian2017nlstf} or an estimate as in \cite{Yokoya2012cnmf}. Details about the CSR function used in MHF \cite{xie2019fusion} are omitted in that paper but we used the matrix provided in the official repository\footnote{https://github.com/XieQi2015/MHF-net/tree/master/CMHF-net}. Since the CSR function varies in different methods, our aim is not to compare base methods against each other, but rather how our framework improves each of their outputs.

For all base methods, we set $\beta = 1$ as the trade-off parameter in (\ref{eq:tv-tv}), and $\rho = 0.2$ as the augmented Lagragian parameter in (\ref{Eq:ADMMSubProbU})-(\ref{Eq:lmult}). The algorithm was terminated when the primal or dual residual was smaller than 0.001, or the number of iterations exceeded 120.

The performance of each method was assessed by the following quantitative measures: peak-signal-to-noise ratio (PSNR), structural similarity (SSIM) \cite{Wang2004ssim}, spectral angle mapper (SAM) \cite{Yuhas1992DiscriminationAS}, relative dimensionless global error in synthesis (ERGAS) \cite{Wald2002DataFD} and root mean squared error (RMSE).

    \begin{figure}[t]
 	\centering
 	\begin{subfigure}[b]{0.25\textwidth}
 		\centering
 		\includegraphics[width=2.5cm, height=4cm]{figs/Haravrd_GTV3}
 		\caption{GT}
 	\end{subfigure}
 	\begin{subfigure}[b]{2.5cm}
 		\centering
 		\includegraphics[width=2cm, height=4cm]{figs/rgblr_harvardv2}
 		\caption{RGB/LR}
 	\end{subfigure}
 	\begin{subfigure}[b]{0.11\textwidth}
 		\centering
 		\includegraphics[width=2cm, height=4cm]{figs/cnmf_harvard}
 		\caption{CNMF}
 	\end{subfigure}
 	\begin{subfigure}[b]{0.11\textwidth}
 		\centering
 		\includegraphics[width=2cm, height=4cm]{figs/cnmf_harvard_ours}
 		\caption{Ours+CNMF}
 	\end{subfigure}
 	\begin{subfigure}[b]{0.11\textwidth}
 		\centering
 		\includegraphics[width=2cm, height=4cm]{figs/usdn_harvard}
 		\caption{uSDN}
 	\end{subfigure}
 	\begin{subfigure}[b]{0.11\textwidth}
 		\centering
 		\includegraphics[width=2cm, height=4cm]{figs/usdn_harvard_ours}
 		\caption{Ours+uSDN}
 	\end{subfigure}
 	\caption{Reconstruction of band 15 from a sample image from the Harvard dataset using CNMF and uSDN as base methods.}
 	\label{fig:harvard_results}
 	\vspace{-0.5cm}
 \end{figure}

\mypar{Experiments on the CAVE and Harvard datasets} Table \ref{tab:harvard} (resp. \ref{tab:cave}) shows the average PSNR, SSIM, SAM, ERGAS and RMSE of \cite{Yokoya2012cnmf, dian2017nlstf, xie2019fusion, qu2018unsupervised} and of our algorithm using those base methods as input, on the Harvard and CAVE datasets. Larger PSNR and SSIM and smaller ERGAS, SAM and RMSE indicate images of better quality. The results clearly show that our framework consistently achieves better metric scores on all the datasets and methods considered, confirming that we are able to obtain better spatial and spectral resolution by enforcing measurement consistency. For example, for the CAVE dataset, the SAM gains range between 0.004 and 0.772, while the improvements of the ERGAS range between 0.005 and 0.183. Also, note that our method provides larger gains for DL based methods, which cannot easily enforce consistency during testing.

 We also present visual results, using the jet colormap for easier visualizations, in Figs. \ref{fig:harvard_results} and \ref{fig:caveresults} on a sample image from each dataset. The first subfigure a) of each figure contains the ground truth (GT) image and a patch extracted from it, while the second subfigure b) presents the RGB and LR image, respectively. Note that the RGB image corresponds to the one used in MHF. Other networks have a different RGB image due to a different CSR function. The rest of the subfigures contain the output from a base method followed by the output from our method. This is repeated for each base method.

 The superiority of our method can be observed from the extracted patches. The base methods tend to blur out details in the image, while our method preserves the original measurements in the outputs. Although the proposed algorithm improves the quality of the outputs of both model-based and learning-based methods in several metrics, one shortcoming is the additional computation it requires. For example, to upsample a $16 \times 16 \times 31$ HS image with a $32 \times$ factor and a $512 \times 512 \times 3$ RGB image, MHF (supervised method) requires 25 seconds, uSDN (unsupervised method) requires 660 seconds while our method requires 840 seconds.

 \begin{figure}[t]
 	\centering
 	\begin{subfigure}[b]{0.25\textwidth}
 		\centering
 		\includegraphics[width=2.4cm, height=4.3cm]{figs/GTV2}
 		\caption{GT}
 	\end{subfigure}
 	\begin{subfigure}[b]{2.5cm}
 		\centering
 		\includegraphics[width=2cm, height=4.2cm]{figs/rgblrv2}
 		\caption{RGB/LR}
 	\end{subfigure}
 	\begin{subfigure}[b]{0.11\textwidth}
 		\centering
 		\includegraphics[width=2cm, height=4cm]{figs/mhf_cave}
 		\caption{MHF}
 	\end{subfigure}
 	\begin{subfigure}[b]{0.11\textwidth}
 		\centering
 		\includegraphics[width=2cm, height=4cm]{figs/Ours_mhf_cave}
 		\caption{Ours+MHF}
 	\end{subfigure}
 	\begin{subfigure}[b]{0.11\textwidth}
 		\centering
 		\includegraphics[width=2cm, height=4cm]{figs/usdn_cave}
 		\caption{uSDN}
 	\end{subfigure}
 	\begin{subfigure}[b]{0.11\textwidth}
 		\centering
 		\includegraphics[width=2cm, height=4cm]{figs/ours_usdn_cave}
 		\caption{Ours+uSDN}
 	\end{subfigure}
 	\begin{subfigure}[b]{0.11\textwidth}
 		\centering
 		\includegraphics[width=2cm, height=4cm]{figs/cnmf_cavev2}
 		\caption{CNMF}
 	\end{subfigure}
 	\begin{subfigure}[b]{0.11\textwidth}
 		\centering
 		\includegraphics[width=2cm, height=4cm]{figs/ours_cnmf_cavev2}
 		\caption{Ours+CNMF}
 	\end{subfigure}
 	\begin{subfigure}[b]{0.11\textwidth}
 		\centering
 		\includegraphics[width=2cm, height=4cm]{figs/nlstf_cave}
 		\caption{NLSTF}
 	\end{subfigure}
 	\begin{subfigure}[b]{0.11\textwidth}
 		\centering
 		\includegraphics[width=2cm, height=4cm]{figs/nlstf_cave_ours}
 		\caption{Ours+NLSTF}
 	\end{subfigure}
 	\caption{Reconstruction of band 15 of the \textit{flowers} image from the CAVE dataset for all the base methods considered.}
 	\label{fig:caveresults}
 	\vspace{-0.5cm}
 \end{figure}

\section{Conclusions}

We proposed a post-processing step for the fusion of RGB and HS images that can work with either classic or DL based methods. Results show that by guaranteeing consistency in the spatial and spectral measurements, the image quality is improved. Moreover, our method can easily work with different downscaling operators and camera spectral response functions, unlike DL networks that are trained for specific operators. The basis of our framework is TV-TV minimization, which we solve with an ADMM-based algorithm. Future lines of work include using different regularizers and unrolling our algorithm in a neural network to improve reconstruction time.

\section{Acknowledgements}
Work supported by Royal Society and NSFC (IEC/NSFC/181255), EPSRC
(EP/T026111/1), National Key R\&D Program of China (Grant
No.2018YFE0207600), the Natural Science Foundation of China
(61911530216) and the Beijing Natural Science Foundation (L202019).

\bibliographystyle{IEEEbib}