\section{Introduction}
\label{Intro}
We define a class of problems $H$ that are constrained discrete optimization problems with the objective function (OF) satisfying the principle of optimality (PO) without considering the constraints \cite{Bellman}. The constraint functions are assumed to be neither convex nor linear in their decision variables. Nor are the constraints required to satisfy the linear independence constraint qualification (LICQ) \cite{Nocedal}. These problems, in general, are NP-Hard \cite{NPHard}. Solving $H$ to optimality with reduced computational resources is highly desirable. Popular methods like Dynamic Programming (DP), Branch and Bound (BB), and Lagrangian relaxation methods can solve $H$ to optimality \cite{DynP1,DynP2,Bbound, Kolman, Fisher} as long as we place either linearity or quadratic conditions on the OF or constraints \cite{Bbound,LinShu}, while many others do not ensure optimality \cite{Plancher}.\\
\indent Examples of $H$ include resource allocation (power, ADC bits) in massive Multiple-Input Multiple-Output (MaMIMO) systems with constraints \cite{Zakir2}, constrained network optimization problems used in routing protocols \cite{NetApp}, and Reinforcement Learning (RL) applied to Markov decision processes with constraints \cite{RLbook}, to name a few.\\

\noindent The contributions of this paper are as follows
\begin{itemize}
\item we propose a novel multi-survivor Dynamic Programming (msDP) algorithm to solve $H$ optimally at significantly reduced computational cost, and,
\item we show that Constrained Markov Decision Processes (CMDP) with finite horizons belong to $H$.
\end{itemize}

\textit{Notations:}
The column vectors are represented as boldface small letters. The superscripts $T$ denote transpose. The terms $\mathbb{I}$ and $\mathbb{R}$, indicate space of integers and real numbers respectively. We represent a vector $\bold{x} = [x_1, x_2,\cdots,x_N]^T$ whose entries are constants up to $x_i$ for $i \le m$ as $\bold{x}^{m}$ and call such vectors partial, i.e., $\bold{x}^{m} = [a_1, a_2,\cdots, a_m, x_{m+1},\cdots, x_N]^T$ where $a_i$'s are constants. The $x_i$'s belong to a finite set $\mathcal{X}$ with cardinality $M$. If there are multiple of such partial vectors under consideration, we refer to the $k^{\text{th}}$ such vector as $\bold{x}_{k}^{m}$ and its elements as $\bold{x}_{k}^{m} = \big[a_{1}^{k},a_{2}^k,\cdots,a_{m}^k, x_{m+1}^k,x_{m+2}^k,\cdots,x_N^k \big]^T$.  The $j^{\text{th}}$ element of the set $\mathcal{X}$ is represented as $\{ \mathcal{X} \}_j$. If the $m^{\text{th}}$ element of $\bold{x}_{k}^{m}$ takes a value $\{ \mathcal{X} \}_j$, we represent such partial vector as $\bold{x}_{k}^{mj}$.\\

\indent The rest of this paper is organized as follows. We define the problem $H$ in Section \ref{Prob}. We first show in Section \ref{Cmdp} that the CMDP framework typically used in RL belongs to $H$. In Section \ref{TF} we describe the theoretical foundations that lead to the proposed msDP algorithm. We describe the proposed algorithm and evaluate its worst-case computational complexity in Section \ref{Algo}. We compare the performance and computational complexity of the proposed method with exhaustive search (ES) and the heuristic algorithm in Section \ref{Sim}, followed by the conclusions in Section \ref{Conc}.

\section{Problem Setup}
\label{Prob}
We define $H$ as stated below, where $\bold{x}^*$ is the optimal solution to $\eqref{cdo_1}$ if it exists.
\begin{equation}\label{cdo_1}
\begin{split}
&\underbrace{\text{max}}_{\substack{\bold{x}}} f(\bold{x}),\\
\text{ such that }&c_i(\bold{x}) \le \alpha_i;\text{ for }1 \le i \le Q_{I},\\
&h_j(\bold{x}) = \beta_j;\text{ for }1 \le j \le Q_E,
\end{split}
\end{equation}
where $\bold{x} = [x_1, x_2,\cdots,x_N]^T \in \mathbb{R}^N$. The OF $f(\bold{x})$ satisfies the PO without the constraint functions $c_i(\bold{x})$ and $h_j(\bold{x})$, where $\alpha_i, \beta_j \in \mathbb{R}, \forall i,j$ \cite{Bellman}. Thus, it can be written in the form $f(\bold{x}) = \sum_{i=1}^{N} b_i \phi_i(x_i)$ where $b_i \in \mathbb{R}$ are constants, $x_i \in \mathcal{X}$ can only take values from the set $\mathcal{X}$ whose cardinality is $M$, and, $\phi_i:\mathcal{X} \longrightarrow \mathcal{Y},\text{ for } 1\le i \le N$. The mapping $\phi_i$ need not be known in closed form. The constraint functions $c_i(\bold{x})$ and $h_j(\bold{x})$ are not limited to linear mappings in $x_i$, nor convex and need not satisfy LICQ \cite{Nocedal}. The terms $Q_{I}$ and $Q_E$ represent the number of inequality and equality constraints, respectively.

\section{Finite Horizon CMDP as problem $H$}
\label{Cmdp}
A finite-state Markov decision process is a 5-tuple $(X,A,P,r,\mu)$ where $X$ is the finite set of states, $A$ is the finite set of actions, $P:X \times A \times X^{'} \rightarrow [0,1]$ are the state transition probabilities $p_{x,a}(x^{'})$ that a state $x^{'}$ is attained when an action $a \in A$ is taken in state $x$. Both $x, x^{'} \in X$. A reward $r:X \times A \rightarrow \mathbb{R}$ is associated with an $a \in A$ from a state $x \in X$. The map $\mu : X \rightarrow [0,1]$ is the starting state distribution \cite{RL2}.\\
\indent A policy $\pi = (D_1, D_2,\cdots,D_{N-1})$ is defined as a set of decisions taken by an agent over states that maximize a certain cumulative reward. A stationary deterministic policy has a decision rule $D_i:X \rightarrow A$ that maps states to actions. In case of randomized policies, the decision rule $D_i$ is a random variable and described using state dependent distributions on possible actions. Hence, $D_i$ is described as $D_i:X \rightarrow \mathcal{P}(A)$ where $\mathcal{P}(A)$ is the probability distribution over actions defined by $q_i(a|x) = \mbox{Prob}[D_i=a|x_i]$ for every action $a \in A$.\\
\indent The cumulative reward is defined as $R = \sum_{i=1}^N \gamma^i r_i$ where $r_i$ is a reward accrued at $i^{\text{th}}$ time instance when the agent uses a policy $\pi$. The discount factor $\gamma \in [0,1]$. We can write the expected sum of discounted rewards as a function of a given policy and an initial state $x$ as $V^{\pi}(x) = \mathbb{E}[R \mid x_1=x] $\cite{RL2}, which can be rewritten as
\begin{equation}\label{cdo_3}
\begin{split}
&V^{\pi}(x) = \mathbb{E}\bigg[ \sum_{i=1}^N \gamma^i r_i\big(x_i,D_i(x_i)\big) \mid x_1=x\bigg],\\
&= \sum_{i=1}^N \gamma^i \hat{r_i}(x_i),\text{ where }\hat{r_i}(x_i) =  \sum_{a \in A} q_i(a|x_i)r_i(x_i,a),
\end{split}
\end{equation}
assuming a random stationary policy. Considering an initial distribution $\mu$ on all the possible starting states on set $X$, we can define the value of the policy as \cite{RL2}
\begin{equation}\label{cdo_5}
\begin{split}
J(\pi) &= \mathbb{E}_{x \sim \mu}\bigg[ V^{\pi}(x)\bigg] = \sum_{x \in X} \mu(x)V^{\pi}(x).\\
&= \sum_{x \in X} \mu(x) \sum_{i=1}^N \gamma^i \hat{r_i}(x_i) = \sum_{i=1}^N \gamma^i \phi_i(x_i),\\
\end{split}
\end{equation}
where $\phi_i(x_i) = \sum_{x \in X} \mu(x)\hat{r_i}(x_i)$, which is obtained by substituting $V^{\pi}(x)$ from \eqref{cdo_3}. We see that with $b_i = \gamma^i$ the OF $J(\pi)$ of \eqref{cdo_5} is similar to the OF $f(\bold{x}) = \sum_{i=1}^{N} b_i \phi_i(x_i)$ of $H$ in \eqref{cdo_1}.\\
\indent For the CMDP case, we have constraint functions that could be either secondary reward functions or penalty functions. The constraint functions $c$ defined as $c:X \times A \rightarrow \mathbb{R}$ such that $c^{\pi}(x) \le d$. The whole optimization problem is written as \cite{RL2}
\begin{equation}\label{cdo_6}
\underbrace{\text{max}}_{\substack{\pi}} J(\pi),\text{ such that }c(\pi) \le d.
\end{equation}
In \eqref{cdo_6}, $c(\pi) = \mathbb{E}[c^{\pi}(x)] = \mathbb{E} \big[ \sum_{i=1}^N \gamma^i c_i(x_i,a_i) \big]$.

\section{Theoretical Foundations of msDP}\label{TF}
In this section, we construct the theoretical foundations for the proposed msDP. We first define the Constraint Satisfaction Function (CSF).
We define a boolean function $A:\bold{x} \rightarrow \{0,1\}$ such that
\begin{equation}\label{cdo_8}
    A(\bold{x})=
\begin{cases}
    1,& \text{if } \bold{x} \text{ satisfies all the constraints }\\
    &c_i(\bold{x}), h_j(\bold{x}) \text{ of H; for all } i,j.\\
    0,              & \text{otherwise.}
\end{cases}
\end{equation}
We also define $A(\cdot)$ on a partially assigned vector $\bold{x}_{k}^{m}$ as $A(\bold{x}_{k}^{m}) = 1$ for $1 \le m < N$ if there exists at least one vector ${\bold{b}_k} = [b_{m+1}^k,b_{m+2}^k,\cdots,b_{N}^k]^T$ such that $b_i \in \mathcal{X}$ for $i > m$ that satisfies all the constraints. This is represented as
\begin{equation}\label{cdo_9}
    A(\bold{x}_{k}^{m})=
\begin{cases}
    1,& \text{if there exists at least one } \bold{b}_k \text{ defined}\\
    &\text{ above, that satisfies all the constraints }\\
    &c_i(\bold{x}_{k}^{m}), h_j(\bold{x}_{k}^{m}) \text{ for all } i,j.\\
    0,              & \text{otherwise.}
\end{cases}
\end{equation}

\subsection{Principle of Optimality}
We define a partially computed OF of $H$ as $\lambda(\bold{x}_{k}^{m}) = \sum_{i=1}^{m} b_i \phi_i(a_i)$ for $1 \le m < N$. It is to be noted that $A(\bold{x}) = A(\bold{x}^N)$ and $\lambda(\bold{x}^{N}) = f(\bold{x})$. Here $\bold{x}^N$ is a completely known vector.\\
\indent Now we show that $H$ does not satisfy the PO. We then introduce a multi-valued function $\mbox{max}_l\{\cdot\}$ that returns $l$ maximum values from its domain. Adopting $\mbox{max}_l\{\cdot\}$, we reformulate $H$ to satisfy the PO.\\

\newtheorem{thrm}{Theorem}
\begin{thrm}\label{thrm1}
If $\bold{x}^*$ is the optimal solution to $H$ described in \eqref{cdo_1} and if there exists at least one infeasible solution $\bold{x}_1$ to $H$ such that $f(\bold{x}_1) > f(\bold{x}^*)$ with constraint satisfaction function being true at some intermediate stage $l$ in the recursion $A(\bold{x}_1^l) > 0$, then it can be shown that $J = \underbrace{\max}_{\bold{x},A(\bold{x}) > 0}f(\bold{x})$ does not obey the PO.
\end{thrm}
\begin{proof}
We say that \eqref{cdo_1} satisfies the PO if we can express it as a value function \cite{Bellman}. That is, given
\begin{equation}\label{cdo_10}
\begin{split}
J &= f(\bold{x}^*) = J(x_1) = \underbrace{\text{max}}_{\substack{\{x_i\}_{i=1}^N,\\A(\bold{x}) > 0}}\Big\{\sum_{i=1}^{N} b_i \phi_i(x_i)\Big\},\text{ then}\\
J(x_i) &= \underbrace{\text{max}}_{\substack{x_i;\\A(\bold{x}^{i})>0}}\Bigg\{ b_i\phi_i(x_i) +  \underbrace{\text{max}}_{\substack{\{x_j\}_{j={i+1}}^N;\\{A(\bold{x}^{i+1}) > 0}}} \sum_{j={i+1}}^{N} b_j \phi_i(x_j) \Bigg\},\\
J(x_i) &= \underbrace{\text{max}}_{\substack x_i;A(\bold{x}^{i})>0}\Bigg\{ b_i\phi_i(x_i) +  J(x_{i+1}) \Bigg\}.
\end{split}
\end{equation}
\indent Let $\bold{x}_1$ be an infeasible solution to $H$ such that
\begin{equation}\label{cdo_10a}
\begin{split}
f(\bold{x}_1) &> f(\bold{x}^{*}),\text{ where}\\
A(\bold{x}_1) &= 0\text{ and }A(\bold{x}^*) > 0.
\end{split}
\end{equation}
However, at some intermediate stage $l < N$ in the recursion, the condition depicted below would arise as a consequence of no conditions placed on the constraints \cite{Nocedal, Boyd}.\\
As $\lambda(\bold{x}_1^l) > \lambda(\bold{x}^{*l}),\text{ such that }A(\bold{x}_1^l) > 0\text{ and }A(\bold{x}^{*l}) > 0: \lambda(\bold{x}_1^l) = \underbrace{\text{max}}_{\substack x_l;A(\bold{x}^{l})>0}\bigg\{ \lambda(\bold{x}_1^l), \lambda(\bold{x}^{*l}), \cdots \bigg\}$.
Thus, we can write the recursion at stage $l$ in \eqref{cdo_10} as $J = \lambda(\bold{x}_1^l) +  J(x_{l+1})$.
As a result, we see that the optimal solution $\bold{x}^{*l}$ is dropped against $\bold{x}_1^{l}$ at stage $l$. However, in the subsequent stage $r$, where $l < r \le N$, another candidate solution $\bold{x}_2^r$ is picked against $\bold{x}_1^r$ because the infeasibility of $\bold{x}_1^{r}$ shows up as $A(\bold{x}_2^r) > 0,\text{ and }A(\bold{x}_1^r) = 0$.
Consequently, $\bold{x}_1^{r}$ is dropped and we have the recursion equation at stage $r \le N$ as
\begin{equation}\label{cdo_14}
\begin{split}
J &= \lambda(\bold{x}_2^r) +  J(x_{r+1}).
\end{split}
\end{equation}
Extending the result to stage $N$, we can see that
\begin{equation}\label{cdo_14}
\begin{split}
J &= f(\bold{x}_2) = \lambda(\bold{x}_2^N) = \underbrace{\text{max}}_{\substack{\bold{x};A(\bold{x}) > 0}}\sum_{i=1}^{N} b_i\phi_i(x_i).
\end{split}
\end{equation}
However, by definition we have
\begin{equation}\label{cdo_14}
J = f(\bold{x}^*) = \lambda(\bold{x}^{*N}) = \underbrace{\text{max}}_{\substack{\bold{x};A(\bold{x}) > 0}}\sum_{i=1}^{N} b_i\phi_i(x_i),
\end{equation}
where $\bold{x}^{*} \ne \bold{x}_2$, hence, $H$ defined in \eqref{cdo_1} does not satisfy the PO.
\end{proof}
We now define a multi-valued function $\mbox{max}_{N_e}:\mathbb{R}^m \rightarrow \mathbb{R}^{N_e}$ such that it returns at-most $N_e$ maximum values over a combination of $m$ variables in its domain. We now rewrite $H$ in \eqref{cdo_1} using $\Gamma$ as defined below
\begin{equation}\label{cdo_15}
\begin{split}
&\underbrace{\text{max}}_{\substack {\bold{x}};A(\bold{x})>0} f(\bold{x}) = J = \text{max}\Big\{ \Gamma(x_1) \Big\},\text{ where}\\
&\Gamma(x_1) = \boldsymbol{\lambda}(\bold{\vec{x}}) = \underbrace{\mbox{max}_{N_e}}_{\substack \{x_j\}_{j=1}^{N};A(\bold{\vec{x}})>0}\Big\{\sum_{i=1}^{N} b_i \phi_i(x_i)\Big\},
\end{split}
\end{equation}
where $\boldsymbol{\lambda}(\bold{\vec{x}})$ is a vector of $N_e$ maximum values. That is, $\boldsymbol{\lambda}(\bold{\vec{x}}) = [\lambda(\bold{x}_1), \lambda(\bold{x}_2),\cdots,\lambda(\bold{x}_{N_e})]^T$ that corresponds to the $N_e$ solutions $\bold{\vec{x}} = [\bold{x}_1,\bold{x}_2,\cdots,\bold{x}_{N_e}]^T$. It is also to be noted that the CSF $A(\bold{\vec{x}})>0$ is true iff $A(\bold{x}_i)>0$ is satisfied for all  $\bold{x}_i \in \bold{\vec{x}}$.\\
\indent We also define the partially evaluated objective at stage $m \le N$ as $\boldsymbol{\lambda}(\bold{\vec{x}}^m) = [\lambda(\bold{x}_1^m), \lambda(\bold{x}_2^m),\cdots,\lambda(\bold{x}_{N_e}^m)]^T = \underbrace{\text{max}_{N_e}}_{\substack \{x_j\}_{j=1}^{m};A(\bold{\vec{x}}^m)>0}\Big\{\sum_{i=1}^{m} b_i \phi_i(x_i)\Big\}$. As before, the CSF  $A(\bold{\vec{x}}^m)>0$ is true iff $A(\bold{x}_i^m)>0$ is satisfied for all  $\bold{x}_i^m \in \bold{\vec{x}}^m$.\\
Given the optimal solution $\bold{x}^*$, we say that the value function $\Gamma(\cdot)$ in \eqref{cdo_15} satisfies the PO if we can write
\begin{equation}\label{cdo_16}
\begin{split}
\Gamma(x_i) &= \underbrace{\mbox{max}_{N_e}}_{\substack x_i;A(\bold{\vec{x}}^{i})>0}\Big\{b_i\phi_i(x_i) + \Gamma(x_{i+1})\Big\},\text{ where}\\
J &= f(\bold{x}^*) = \text{max}\Big\{ \Gamma(x_1) \Big\}.
\end{split}
\end{equation}
That is, the optimal solution is maintained at every stage of the recursion in $\Gamma(x_i)$.
\begin{thrm}\label{thrm2}
If $\bold{x}^*$ is the optimal solution to $H$ described in \eqref{cdo_1}, then with an appropriate choice of $N_e$ that is dependent on the constraint satisfaction function $A(\cdot)$, it can be shown that $\Gamma(x_1) = \underbrace{\max_{N_e}}_{\substack \{x_j\}_{j=1}^{N};A(\bold{\vec{x}})>0}\Big\{\sum_{i=1}^{N} b_i \phi_i(x_i)\Big\}$ always satisfies the PO.
\end{thrm}
\begin{proof}
Let $\bold{x}_1$ be an infeasible solution to $H$ as described in \eqref{cdo_10a} in Theorem \ref{thrm1}. Let the PO \eqref{cdo_16} hold true; we can then write the same at recursion $l \le N$ as
\begin{equation}\label{cdo_17}
\begin{split}
\Gamma(x_1) &= \underbrace{\mbox{max}_{N_e}}_{\substack x_l;A(\bold{\vec{x}}^l)>0}\Big\{b_l \phi_l(x_l)  + \Gamma(x_{l+1})\Big\},\\
&= \boldsymbol{\lambda}(\bold{\vec{x}}^l) + \Gamma(x_{l+1}).
\end{split}
\end{equation}
We now let $N_e$ be large enough such that $\boldsymbol{\lambda}(\bold{\vec{x}}^l) = [\lambda(\bold{x}_1^l), \lambda(\bold{x}^{*l}),\cdots]^T$. That is $(\bold{x}_1^l, \bold{x}^{*l}) \in \bold{\vec{x}}^l$.\\
\indent However, at some later stage $r$ where $l < r \le N$, we can write the recursion as
\begin{equation}\label{cdo_18}
\begin{split}
\Gamma(x_1) &= \underbrace{\mbox{max}_{N_e}}_{\substack x_r;A(\bold{\vec{x}}^r)>0}\Big\{b_r \phi_l(x_r)  + \Gamma(x_{r+1})\Big\},\\
&= \boldsymbol{\lambda}(\bold{\vec{x}}^r) + \Gamma(x_{r+1}).
\end{split}
\end{equation}
Now we note that $\lambda(\bold{x}_2^r) \in \boldsymbol{\lambda}(\bold{\vec{x}}^r)$ and $\lambda(\bold{x}_1^r) \not\in \boldsymbol{\lambda}(\bold{\vec{x}}^l)$. This is because, at this stage in the recursion, the infeasibility of the candidate solution $\bold{x}_1^r$ becomes known as indicated by $A(\bold{x}_1^r) = 0$. Thus the solution $\bold{x}_1^r$ is not part of $\bold{\vec{x}}^r$.\\
With the appropriate choice of $N_e$, it can be ensured that $(\bold{x}_2^r, \bold{x}^{*r}) \in \bold{\vec{x}}^r$. The determination of $N_e$ is discussed in Lemmas \ref{lemma3} and \ref{lemma4} in section \ref{TF}. Extending $r$ to $N$, we have $(\bold{x}_2^N, \bold{x}^{*N}) \in \bold{\vec{x}}^N$. Hence, we can write
\begin{equation}\label{cdo_19}
\begin{split}
\Gamma(x_1) &= \boldsymbol{\lambda}(\bold{\vec{x}}^N) =  \underbrace{\mbox{max}_{N_e}}_{\substack \{x_j\}_{j=1}^{N};A(\bold{\vec{x}}^N)>0}\Big\{\sum_{i=1}^{N} b_i \phi_i(x_i)\Big\},\\
&\text{with }\bold{x}^{*} \in \bold{\vec{x}}^N;\text{ Thus, }f(\bold{x}^*) = \text{max}\Big\{ \Gamma(x_1) \Big\}.
\end{split}
\end{equation}
Thus, with an appropriate choice of $N_e$, we see that the PO for \eqref{cdo_15} is always satisfied.
\end{proof}

\subsection{Trellis for H}\label{graph}
We define a graph $G(V,E,A)$ with nodes $V$, edges $E$, and edge label $A$ as
\begin{equation}\label{mlba_4}
\begin{split}
V_i &= \Big\{ x_{ij} \mid x_{ij} = \{ \mathcal{X} \}_j; 1 \le i \le N; 1 \le j \le M \Big\} \cup S \cup F,\\
E &=\Big\{(x_{ij},x_{{i+1}m}) \mid x_{ij} = \{ \mathcal{X} \}_j; x_{{i+1}m} = \{ \mathcal{X} \}_m;\\
&1 \le j,m \le M; 1 \le i \le N-1; (S,x_{1j}),  (x_{Nm},F) \Big\},\\
A &= \Big\{ \psi_i(x_{i+1} \mid x_{i}) =  \phi_i(x_i) \leftarrow (x_{i},x_{i+1}) \mid \forall x_i \in \mathcal{X};\\
&1 \le i \le N; \phi_N(x_N) \leftarrow (x_N,F) \mid \forall x_N \in \mathcal{X}; \\
& 0  \leftarrow (S,x_1) \mid \forall x_1 \in \mathcal{X} \Big\}.
\end{split}
\end{equation}
The nodes $S$ and $F$ are the starting and terminal nodes, respectively. The optimal solution $\bold{x}^* = [x_{1j},x_{2j},\cdots,x_{Nj}]^T$ can be thought of as a walk through this graph starting at node $S$ and ending at the destination node $F$ passing through the nodes $x_{ij} \in V_i$ from stage $i=1$ to $i=N$, along the edges $(x_{ij},x_{{i+1}m}) \in E$ accruing a reward $\psi_i(x_{{i+1}m} \mid x_{ij}) = \phi_i(x_{ij}) \in A$ that results in the maximum total accumulated reward $f(\bold{x}^*)$ such that all the constraints are satisfied. That is  $A(\bold{x}^*) > 0$. It can also be noted that this graph $G(V,E,A)$ can be qualified as a trellis  $T(V,E,A)$ \cite{Kolman}. An example $T(V,E,A)$ with $M=4$ is shown in Fig \ref{Fig2_mlba}.\\
\begin{figure*}[http]
\centering
\resizebox{.7\textwidth}{!}{
\input{Fig1_ICML2020_t.eps}
}
\caption{An example trellis $T(V,E,A)$ for $M=4$.}
\label{Fig2_mlba}
\end{figure*}
\indent A DP approach can not be used to obtain $\bold{x}^*$ as shown using Theorem \ref{thrm1}. We note from Theorem \ref{thrm2} that the PO is satisfied by using the multi-valued function $\mbox{max}_{N_e}$. Thus, maintaining multiple survivor paths $(N_e)$ at each stage of the trellis ensures optimality if $N_e$ is picked properly. An appropriate choice of  $N_e$ is discussed in the next subsection. Thus the DP algorithm is modified to perform Add, Compare, and Multiple Select (ACMS) as basic operation at each node in the trellis as opposed to Add, Compare, and Select (ACS). The algorithm is described in section \ref{Algo}.

\subsection{Determination of $N_e$}\label{SecNe}
In the Lemmas below, we establish a relationship between the number of survivors $N_e$ and CSF $A(\cdot)$ to ensure optimality. Due to the page limitation, the proofs are omitted.
\newtheorem{lemm}{Lemma}
\begin{lemm}\label{lemma3}
The number of survivor paths (or partial solutions) that need to be maintained at a given node $x_{ij}$ in trellis $T(V,A,E)$ to obtain optimal solution to $H$ is dependent on the constraint satisfaction function $A(\bold{x}^{ij})$ at the node $x_{ij}$ in the trellis. It can be shown that
\end{lemm}
\begin{equation}\label{cdo_lem1}
B(x_{ij}) \triangleq \sum_{k=1}^{L_{ij}} A(\bold{x}^{ij}_k),
\end{equation}
where $L_{ij}$ is the number of paths incident on the node $x_{ij}$, out of which $B(x_{ij})$ paths satisfy all the constraints as per the definition of CSF in \eqref{cdo_9}.

\begin{lemm}\label{lemma4}
The minimum number of maximum paths $N_e$ that need to be maintained in \eqref{cdo_16} at each stage in the recursion for optimality of $H$ can be expressed as a lower bound:
\end{lemm}
\begin{equation}\label{cdo_lem4_3}
N_e \ge \underbrace{\text{max}}_{\substack i} \sum_{j=1}^M \sum_{k=1}^{L_{ij}} A(\bold{x}^{ij}_k).
\end{equation}

\section{Multi-survivor Dynamic Programming Algorithm}
\label{Algo}
\begin{algorithm}[h!]
  \caption{Multi-survivor Dynamic Programming.}\label{AlgoMVA}
  \begin{algorithmic}
    \small
    \State {\bfseries Input:}{$M$,$N$,$A(\cdot)$,$\mathcal{X}$}
      \State $M\gets \text{Cardinality of sets }\mathcal{X}$
      \State $N\gets \text{Number of stages or variables of }f(\cdot)\text{ of problem H}$
      \State $A(\cdot) \gets \text{Constraint satisfaction function of H}$
      \State $\mathcal{X} \gets \text{Set of discrete values that }x_i\text{'s can take }$
      \State \text{Determine the path valency $E(x_{ir})$ for all nodes}
      \State \text{Initialize paths and its associated costs}
      \State $\scriptstyle \bold{x}^{1,1} = \{\mathcal{X}\}_1, \bold{x}^{1,2} = \{\mathcal{X}\}_2, \cdots, \bold{x}^{1,M} = \{\mathcal{X}\}_M$
      \State $\scriptstyle \lambda(\bold{x}^{1,1})=b_1\phi_1(\{\mathcal{X}\}_1), \lambda(\bold{x}^{1,2})=b_1\phi_1(\{\mathcal{X}\}_2), \cdots,$
      \State $\scriptstyle \lambda(\bold{x}^{1,M})=b_1\phi_1(\{\mathcal{X}\}_M)$
      \For{$i=2:N$}
           \For{$r=1:M$}
           \State{$m \gets 1$}
           	\For{$b=1:M$}
		   \If {$(x_{{i-1}b},x_{ir}) \in T(A)$}
			\For{$k=1:E(x_{ir})$}
			   \If {$A(\bold{x}_k^{ir}) > 0$}
			       \State{$x_{ir} =  \{\mathcal{X}\}_r$}
			       \State{$\bold{x}^{i,r}_k = [(\bold{x}^{i-1,b}_k)^T, x_{ir}]^T$}
    			       \State{$\lambda(\bold{x}^{i,r}_k) = \lambda(\bold{x}^{i-1,b}_k) + b_i\phi_i(x_{ir})$}
			   \EndIf
			\EndFor
		   \EndIf
		\EndFor
		\State{\text{Keep only $\{\bold{x}_m^{i,r}\}_{m=1,\cdots,E(x_{i,r})}$ paths, that }
		\State{$\text{ correspond to }E(x_{ir})$ maximum values of}
		\State{$\scriptstyle \Big\{\lambda(\bold{x}_m^{i,r})\Big\}_{m=1,\cdots,E(x_{ir})}$ and ignore the rest!}}
           \EndFor
     \EndFor
     \State{Process the termial node "F"}
     \State{$m \gets 1$}
     \For{$b=1:M$}
		\For{$k=1:L_{T}\text{ All paths incident on terminal node F}$}
			\State{$\scriptstyle \lambda(\bold{x}_m^F) = \lambda(\bold{x}^{N-1,b}_k)+b_N\phi_N(x_{Nb})$}
			\State{$\scriptstyle m \gets m+1$}
		\EndFor
     \EndFor
     \State{\text{Pick $\bold{x}^*$ that has maximum of $\Big\{\lambda(\bold{x}^F_k)\Big\}_{k=1,\cdots,m-1}$ }}
     \State{$\textbf{return }\bold{x}^* \text{ Optimal Solution}$}
  \end{algorithmic}
\end{algorithm}

\subsection{Evaluation of CSF}\label{const}
It can be shown that the $A(\bold{x}^m)$ can be evaluated at any intermediate stage $m \le N$ by posing it as an unconstrained discrete optimization (UDO) problem. The UDO problems can be solved in polynomial time using graph structures like Depth-First-Search (DFS), Breadth-First-Search (BFS), or BB algorithms \cite{Bbound, Ndale}.

\subsection{Computational Complexity Analysis}
\label{CCA}
The computational blocks of the msDP involve (i) evaluating the CSF $A(x^{ij}_k)$, and (ii) performing the ACMS operation on the $E(x_{ij})$ paths incident on a given node $x_{ij}$ of the trellis. For (i), a BFS or DFS can be used in such cases which has complexity of $O(E + V)$, where $E$ and $V$ are the numbers of edges and vertices of the trellis \cite{Ndale}. In the worst case of a fully connected trellis, $E = (N-1)M^2,$ and $V = NM$. Also, in the worst case scenario of $N_e^{'} = \frac{N_e}{M}$ incident on all the nodes in the given trellis, the net complexity of CSF evaluation $T_{CSF}$ can be expressed as
\begin{equation}\label{cdo_cca1}
T_{CSF} = N_e\Big((N-1)NM^2+N^2M\Big)\text{ for all the }NM\text{ nodes}.
\end{equation}
In case of (ii), the total number of ACMS evaluations for the trellis in the worst case is $T_{ACMS} = \frac{N_e}{M}NM = N_eN$. Thus, the complexity of msDP is $O(N_eN^2M^2)$ or in general $O(N_eN^{l})$, where $l \in \mathbb{R}$. On the other hand, the ES requires an evaluation of the CSF and OF for $M^N$ solutions. Hence its complexity is $O(M^N)$. In the case of learning techniques like RL, to guarantee an optimal solution, the complexity required is $\approx O(M^N)$ \cite{Watkins}.

\section{Simulations}
\label{Sim}
In this section, we apply the proposed msDP to solve (i) Analog to Digital Converter (ADC) Bit-Allocation (BA) problem for MaMIMO receiver, and (ii) DNA fragment assembly (DFA) problem that is a part of DNA sequencing \cite{Maria,DFA2}. We evaluate the performance and compare the number of computations\footnote{Number of computations = Actual number of CSF evaluations + ACMS operations.} against the ES and the heuristic suboptimal simulated annealing (SA) algorithms\footnote{The optimality is compromised against complexity with heuristic algorithms like SA.} \cite{SimAn2}.

\subsection{ADC Bit Allocation for MaMIMO}
The ADC BA problem is to assign the number of bits to be used by Variable-Resolution ADCs on different Radio Frequency (RF) paths of the MaMIMO receivers. An optimal BA ensures that the performance, of the receiver is maximized under a non-linear power constraint. It is to be noted that the OF is non-linear. In \cite{Zakir2}, the authors reduce this to a problem in $H$, which is described as
\begin{equation}\label{cdo_sim_ba1}
J = \underbrace{\text{max}}_{\substack{\{x_i\}_{i=1}^N;A(\bold{x}) > 0}}\Big\{\sum_{i=1}^{N} \frac{a_i^2}{b_i^2+d_i2^{x_i}}\Big\},
\end{equation}
where $a_i$, $b_i$, and $d_i$ are constants $\in \mathbb{R}$ that represent channel singular value, noise power, and coefficient of quantization noise due to bit allocation $x_i$ on the $i^{th}$ RF path, respectively. Here $N$ is the number of RF paths in the receiver. The CSF $A(\bold{x}) > 0$ iff the power constraint $\sum_{i=1}^{N} 2^{x_i} \le P_t$, and bit-ordering  constraint $x_1 \ge x_2  \ge \cdots  \ge x_N$ are satisfied. The total ADC power budget is $P_t$. Hence we have
\begin{equation}\label{cdo_sim3}
    A(\bold{x}) =
\begin{cases}
    1,\text{ if } &\sum_{i=1}^{N} 2^{x_i} \le P_t,\\
       &x_1 \ge x_2  \ge \cdots  \ge x_N.\\
    0, &\text{ Otherwise}.\\
\end{cases}
\end{equation}
We perform simulations considering the parameters used by the authors in \cite{Zakir2}. They consider a maximum bit resolution of $M=4$, such that $x_i$'s (bits) take values from the set $\mathcal{X} = \{1,2,3,4\}$ for $1 \le i \le N$. We assume the number of RF paths in the receiver as $N=12$ and the power budget of $P_t=48$. For this case, we achieve $N_e \ge 13$. The number of computations with msDP, ES, and SA algorithms are shown in Table \ref{tab2}. Both the brute-force and msDP methods obtain the optimal solution of $\bold{x}^* = [4,2,2,2,2,2,1,1,1,1,1,1]^T$ \cite{Zakir5}.

\subsection{DNA Fragment Assembly}
We apply our proposed formalism to the DFA problem, which has equivalence to the TSP \cite{Maria, boundTSP}. DFA is the challenging process of the DNA sequencing \cite{Jay}. DNA sequencing's main problem is that the current technology can not read an entire genome in one shot, sometimes not even more than 1000 bases. Even the simplest organisms, like bacteria and viruses, have much longer genome lengths. Consequently, the genomes are broken down into smaller readable fragments and sequenced \cite{Sanger}. In this step, $N$ copies of DNA are created.  A short fragment is derived from each of the replicated DNA at some random location. These short fragments are sequenced. The final and challenging step is to assemble these sequenced fragments to obtain the original DNA sequence. This step is called the DFA, which is illustrated through an example below \cite{Maria}.\\
\indent We assume the DNA sequence to be $TTACCGTGC$, and the fragments sequenced using 4 DNA copies being $F_1 = ACCGT$, $F_2 = CGTGC$, $F_3 = TTAC$, and $F_4 = TACCGT$. The overlap of each fragment with the other three fragments are computed using similarity measure. Based on this similarity measure, the order of fragments are determined which in the case of this example is $F_3F_4F_1F_2$.\\
\indent The DFA problem is posed as a maximization problem where the sum of the similarity measures\footnote{The similarity scores between different pairs of fragments are derived using the Smith-Waterman algorithm \cite{Maria, Smith}.} between two adjacent DNA fragments is maximized. This is subject to the constraint that there is no repetition of the fragments in the sequence. Formally, this problem is defined as \cite{Maria}
\begin{equation}\label{cdo_sim1}
J = \underbrace{\text{max}}_{\substack{\{F_{\sigma_i}\}_{i=1}^N;A(\mathcal{F}) > 0}}\Big\{\sum_{i=1}^{N-1} \phi(F_{\sigma_i},F_{\sigma_{i+1}})\Big\},
\end{equation}
where $\mathcal{F} = \{ F_{\sigma_1}, F_{\sigma_2},\cdots,F_{\sigma_N} \}$ is the optimal solution indicating the original DNA sequence, $\sigma_i$ is the fragment index, and $\phi(F_{\sigma_i},F_{\sigma_{i+1}})$ is the similarity measure between the fragments $F_{\sigma_i}$ and $F_{\sigma_{i+1}}$. Here, the set $\mathcal{X}$ is collection of DNA fragments $\{F_{\sigma_j}\}_{j=1}^{N}$, where $N$ is the number of fragments. In this example we have $M=N$. For this problem the CSF $A:\mathcal{F} \rightarrow \{0,1\}$ is defined on the set $\mathcal{F}$ as $A(\mathcal{F}) > 0$ iff all the fragments in $\mathcal{F}$ are unique. The CSF denoted as
\begin{equation}\label{cdo_sim3}
    A(\mathcal{F}) =
\begin{cases}
    1,&\text{ if }\bigcap\limits_{i=1}^{N} F_{\sigma_i} = \varnothing,\\
    0, &\text{ Otherwise}.\\
\end{cases}
\end{equation}
We perform the DFA of a small section of the DNA sequence of bacterium \textit{Escherihia Coli (E. coli)} \cite{Maria}. The original section of the DNA is represented as  $TACTAGCAATACGCTTGCGTTCGGT$. We consider $N=10$ fragments each with 8 bases, as follows: $F_1 = \textit{ACGCTTGC}$, $F_2 = \textit{TTGCGTTC}$, $F_3 = \textit{ACTAGCAA}$, $F_4 = \textit{CGTTCGGT}$, $F_5 = \textit{AGCAATAC}$, $F_6 = \textit{TACTAGCA}$, $F_7 = \textit{AATACGCT}$, $F_8 = \textit{CTTGCGTT}$, $F_9 = \textit{ATACGCTT}$, and $F_{10} = \textit{CTAGCAAT}$.
Using the proposed method, we obtain the optimal solution $\mathcal{F} = \{F_6 F_3 F_{10} F_5 F_7 F_9 F_1 F_8 F_2 F_4 \}$. It can be shown that by imposing a second constraint using an upper bound \cite{boundTSP}, we achieve $N_e \ge 454$. Solving the DFA problem to optimality with the imposition of these additional constraints is not possible using the traditional DP. The results are shown in Table \ref{tab3}.

\begin{table}[t]
\caption{Comparison of the total actual number of computations $^\mathsection$ required for msDP, ES, and SA algorithms.}
\label{tab2}
\begin{center}
\begin{sc}
\resizebox{\columnwidth}{!}{%
\begin{tabular}{ |c|c|c|c|}
\hline
& \small Proposed & \small Exhaustive & \small Simulated \\
& \small Method  & \small Search & \small annealing \\
& \tiny $O(N_eN^2M^2)$  & \tiny $O(M^N),O(N!)$ & \tiny $\approx O(M^N),O(N!)^{\ddagger}$ \\
\hline
\small ADC BA  & \small 6912 $(N_e = 13)$ & \small $16.7 \times 10^6$ & \small 5220 \\
\hline
\small DFA-1{$^{\dagger}$} & \small 83890 $(N_e = 434)$ & \small $3.62 \times10^6$ & \small 380 \\
\hline
\end{tabular}%
}
\end{sc}
\tiny{$^{\dagger}$} DFA with a second constraint. \tiny{$\mathsection$ Number of computations = Actual number of CSF evaluations + ACMS operations}. \tiny{$^{\ddagger}$} To guarantee optimal solution.
\end{center}
\vskip -0.1in
\end{table}

\begin{table}[h!]
\caption{Comparison of performance of the msDP, ES, and SA algorithms.}
\label{tab3}
\begin{center}
\begin{sc}
\resizebox{\columnwidth}{!}{%
\begin{tabular}{ |c|c|c|c|}
\hline
& \small Proposed & \small Exhaustive & \small Simulated \\
& \small Algorithm  & \small Search & \small annealing\\
& \small \textbf{\textcolor{blue}{Optimal}}  & \small \textbf{\textcolor{blue}{Optimal}} & \small \textbf{\textcolor{red}{Suboptimal}} \\
\hline
\small ADC BA  & \small \textcolor{blue}{$[4,2,2,2,2,2,1,1,1,1,1,1]^T$}  & \small \textcolor{blue}{$[4,2,2,2,2,2,1,1,1,1,1,1]^T$} & \small \textcolor{red}{$[4,4,2,2,1,1,1,1,1,1,1,1]^T$} \\
\hline
\small DFA-1 \small{$^{\dagger}$} & \small \textcolor{blue}{$[6,3,10,5,7,9,1,8,2,4]^T$} & \small \textcolor{blue}{$[6,3,10,5,7,9,1,8,2,4]^T$} & \small \textcolor{red}{$[6,3,10,5,2,8,4,7,9,1]^T$} \\
\hline
\end{tabular}%
}
\end{sc}
\end{center}
\vskip -0.1in
\end{table}

\section{Conclusion}
\label{Conc}
Constrained discrete optimization problems with the objective function that satisfies the optimality principle and place no conditions on the constraints are NP-Hard. In general, such problems arise in many areas of science, engineering, finance, and management. We propose a multi-survivor dynamic programming framework that solves these problems to optimality with worst-case complexity reduction from either $O(K^N)$ or $O(N!)$ to $O(N_eN^l)$, where $N_e$ and $l$ are constants that depend on the constraints.

\section*{Acknowledgment}
The authors would like to thank National Instruments for the financial support extended for this work.
\bibliographystyle{IEEEtran}


\section{Introduction}
Cyber-Physical System (CPS) and Internet of Things (IoT) are playing an increasingly important role in critical infrastructures and everyday life.
The related area of CPS and IoT continues to emerge and expand as costs drop and the confluence of sensors, platforms and networks increases\cite{2018DHS_report}.
Simultaneously, the cyber-security risks and attack surfaces are also increasing \cite{cardenas2008research} since CPS relies on remote sensing devices, communication channels, and spatially distributed processors, which are prone to failures under cyber attacks on the data acquisition and communication channels.
The unintentional faults or malicious attacks could cause severe system damage, economic loss, and environmental degradation, e.g., the Stuxnet launched on Iran’s nuclear facilities\cite{STUXNET}, power blackouts in Brazil \cite{samba_stop} and Ukraine \cite{Ukraine_Blackout}, etc. The research community has recognized the importance of CPS security, especially the design of secure detection, estimation, and control strategy\cite{cardenas2009challenges}.

Recently, substantial research efforts have been devoted to secure state estimation against malicious sensors.
One of the research paths is the finite horizon approach, where only measurements in a finite time-window are considered when recovering the system state.
The problem is combinatorial in nature since searching for the corrupted sensors involves minimizing the $\ell_0$ norm \cite{FawziTAC2014}, which is an NP problem.
Fawzi et al.~\cite{FawziTAC2014} propose an algorithm using its convex relaxation, i.e., $\ell_1$ optimization, to solve the problem in the absence of noise.
Similarly, Shoukry and Tabuada \cite{ShoukryTAC2016} adopt a $2$-norm batch optimization approach to solve the state estimation problem.
Shoukry et al.~\cite{Shoukry2017} propose an algorithm in the presence of noise by searching for reliable sensors using consistency check, and the searching complexity is reduced based on Satisfiability Modulo Theory.
However, in these works, the sensory data out of the time window are discarded, which may cause performance degradation and estimation delay.

Another solution is the switch estimator\cite{Mishra2017TCNS}\cite{yorie}\cite{luTAC2019} where multiple estimates are maintained based on measurements from different subsets of sensors, and the system operator switches between these estimates based on the evaluation of their reliability by consistency checking or malicious detection algorithms. However, the combinatorial nature of estimation candidates or sensor subsets may incur heavy computation and storage burden of the devices.
In view of these problems, Liu et al. \cite{liuxinghua-TAC2020} design local estimators whose weighted sum coincides with the Kalman estimate with a certain probability in the absence of attack. The local estimates are fused securely by a quadratic programming problem with an $\ell_1$ term to handle the sparse outliers.

Even though there are efficient algorithms computing the estimation by introducing $\ell_1$ relaxation, the design of the secure state estimation is still computational challenging. For secure state recovery problem in the presence of $p$ malicious sensors, it is required that the system needs to be $2p$-sparse observable\cite{ShoukryTAC2016}, and calculating the sparse observability index for general systems is proved to be NP-hard\cite{mao2021computational}. Besides sparse observability, other results \cite{FawziTAC2014}\cite{liuxinghua-TAC2020}\cite{sandberg_TAC2014} impose stronger conditions for system resiliency, whose validations are also NP problems.
However, it has been observed by Mao et al. \cite{mao2021computational} that when the eigenvalues of the system matrix $A$ have unitary multiplicity ($A$ is non-derogatory), the computational complexity of secure state reconstruction without noise can be significantly reduced.
Leveraging upon this assumption, we propose our secure estimation scheme in the formulation of LTI system with Gaussian noise, and it has the following merits:
\begin{itemize}[left=0pt]
	\item In the presence of $p$ compromised sensors, the proposed estimation is secure if the system is $2p$-sparse observable. In the absence of attack, the proposed estimation coincides with Kalman estimation with certain probability, which can be adjusted.
	\item During the designing phase, the sparse observability index can be computed with low complexity. During the algorithm operating phase, the proposed estimation is formulated as the solution of a convex optimization problem which can be computed efficiently.
\end{itemize}

\textit{Organization:} We introduce the problem formulation and preliminary results in Section \ref{sec:problem}. The main results are provided in Section \ref{sec:main_result} and collaborated by numerical simulation in Section \ref{sec:sim}. Section \ref{sec:conclusion} finally concludes the paper.

\textit{Notations:}
Cardinality of a set ${\mathcal{S}}$ is denoted as $|{\mathcal{S}}|$. Conjugate transpose of a matrix $A$ is denoted as $A'$.
The determinant of a matrix is represented by $\det(\cdot)$.
Diagonal matrix with diagonal elements $A_1,\cdots,A_k$ is denoted as $\text{diag}(A_1,\cdots,A_k)$.
Denote the span of row vectors of matrix $A$ as $\rs(A)$.
All-one vector with size $m\times 1$ is denoted as $\mathbf{1}_{m}$.
${\mathbb{R}}^n$ is the space of real vectors with $n$ entries.
The set of complex matrices with $m$ rows and $n$ columns is denoted as ${\mathbb{C}}^{m\times n}$.

\section{Problem Formulation and Preliminary Results}\label{sec:problem}
\subsection{Secure dynamic state estimation}
In this paper, we consider the linear time-invariant system with Gaussian noise:
\begin{align}
x(k+1)&=A x(k)+w(k) , \label{eq:system} \\
y(k)&=C x(k)+v(k)+a(k) ,\label{eq:y_i_def}
\end{align}
where $x(k) \in \mathbb{R}^{n}$ is the system state, $w(k) \sim {N}(0, Q)$ and $v(k) \sim {N}(0, R)$ are i.i.d. Gaussian process noise and measurement noise with zero mean and covariance matrix $Q$ and $R$.
The vector $y(k)\in \mathbb{R}^{m}$ is the collection of measurements from all $m$ sensors, and $i$-th entry $y_i(k)$ is the measurement from sensor $i$.
The vector $a(k)$ denotes the bias injected by an adversary and $a_i(k)$ is the attack on sensor $i$. Define $$z(k)=C x(k)+v(k)$$ as the true measurements without the attack.
The initial state $x(0) \sim {N}(0, \Sigma)$ is assumed to be zero mean Gaussian random variable and is independent from the process noise $\{w(k)\}$.

The secure dynamic estimation problem aims at recovering system state $x(k)$ at every time $k$ based on all history observations $\{y(t),0\leq t\leq k\}$ which have been partly manipulated by the malicious attacker.
It is conventional in the literature \cite{FawziTAC2014}\cite{Shoukry2017} that the attacker can only compromise a fixed subset of sensors with known maximum cardinality.
Denote the index set of all sensors as ${\mathcal{R}} \triangleq\{1,2, \ldots, m\}$.
Define the support of vector $a\in{\mathbb{R}}^{n}$ as $\supp(a)\triangleq \left\{i| 1\leq i\leq n , a_i\neq0 \right\}$ where $a_i$ is the $i$-th entry of vector $a$.
We have the following assumptions on the malicious adversary.

\begin{definition}[Sparse Attack]\label{def:attack}
	The attack is called a $p$-sparse attack if the vector sequence $a(k)$ satisfies that,
	there exists a time invariant index set ${\mathcal{I}}\subseteq {\mathcal{R}} $ with $|{\mathcal{I}}| = p$ such that $\bigcup_{k=1}^{\infty} \supp\left\{a(k)\right\} = {\mathcal{I}}$.
\end{definition}

Closely related to the sparse attack, we introduce the notion of sparse observability that characterizes the system observability in the presence of attack.

\begin{definition}[Sparse observability]\label{df:sparse_obs}
	The sparse observability index of system \eqref{eq:system}-\eqref{eq:y_i_def} is the largest integer $s$ such that system\footnote{The matrix $C_{{\mathcal{R}}\setminus{\mathcal{I}}}$ represents the matrix composed of rows of $C$ with row index in ${\mathcal{R}}\setminus{\mathcal{I}}$.} $(A,C_{{\mathcal{R}}\setminus{\mathcal{I}}})$ is observable for any set of sensors ${\mathcal{I}}\subset{\mathcal{R}}$ with cardinality $|{\mathcal{I}}| = s$. When the sparse observability index is $s$, we say that the system with pair $(A,C)$ is $s$-sparse observable.
\end{definition}
Define $y(k_1:k_2)$ as the sequence $\{y(k_1),y(k_1+1),\cdots,y(k_2)\}$. Similar notation is also applied on $z(k)$.
For linear Gaussian noise system, the estimation is secure if the estimation error is bounded by a constant term irrelevant to the attack.

\begin{definition}[Secure estimator]\label{def:resi}
	An estimator is an infinite sequence of mappings $g=\{g_k\}_{k=1}^{{\infty}}$ where $g_k$ is a mapping from all the history observations to an estimation at time $k$:
	$$g_k\left(y(0:k)\right)=\hat{x}(k).$$
	Define the estimation difference introduced by attack as
	$$d_k\triangleq \left\|g_k\left(z(0:k)\right)-g_k\left(y(0:k)\right) \right\|_2 .
	$$
	The estimator is said to be secure against the $p$-sparse attack if the following holds:

	$$\sup_{k\in{\mathbb{Z}}^+} {\mathbb{E}} \left[d_k^2\right] < {\infty} ,$$
	where ${\mathbb{E}}$ is the expectation with respect to the probability measure generated by the Gaussian noise $\{w(k)\}$ and $\{v(k)\}$.
\end{definition}
If all sensors are benign, i.e., $a(k)=\mathbf{0}$ for all $k$, the optimal state estimator is the classical Kalman filter:
\begin{align*}
\hat{x}(k)&=\hat{x}(k | k-1)+K(k)\left[y(k)-C \hat{x}(k  | k-1)\right] ,\\
P(k)&=P(k  | k-1)-K(k) C P(k  | k-1),
\end{align*}
where
\begin{align*}
&\hat{x}(k  | k-1)=A \hat{x}(k-1), P(k  | k-1)=A P(k-1) A{'}+Q ,\\
&K(k)=P(k  | k-1) C{'}\left(C P(k  | k-1) C{'}+R\right)^{-1},
\end{align*}
with initial condition $\hat{x}(0  |-1)=0,\ P(0  |-1)=\Sigma $.
It is well-known that for observable system, the estimation error covariance matrices $P(k)$ and the gain $K(k)$ will converge to the following matrices $P$ and $K$:
\begin{align*}
&P \triangleq \lim _{k \rightarrow \infty} P(k),\ P_{+}=A P A{'}+Q ,\\
&K \triangleq P_{+} C{'}\left(C P_{+} C{'}+R\right)^{-1}.
\end{align*}
Since typically the control system will be running for an extended period of time, we focus on the case where the Kalman filter is in steady state, and thus the Kalman filter reduces to the following fixed-gain linear estimator:
\begin{equation}\label{eq:fix_gain_kalman}
\hat{x}(k+1)=(A-K C A) \hat{x}(k)+K y(k+1) .
\end{equation}
Before introducing our work, we first recall some results in the previous work that decomposes the fix gain Kalman filter to local estimates and recovers it by an optimization problem.

\subsection{Preliminary Results}\label{sec:preli}
The following preliminary results in this subsection are from \cite{liuxinghua-TAC2020}.
We introduce the following assumption:
\begin{assumption}\label{as:distinct_eigvalue}
	$A-K C A$ has $n$ distinct eigenvalues. Moreover, $A-K C A$ and $A$ do not share any eigenvalue.
\end{assumption}
Since $A-K C A$ has distinct eigenvalues, it can be diagonalized as:
\begin{equation}\label{eq:VLambda}
A-K C A=V \Pi V^{-1},
\end{equation}
where $\Pi$ is a diagonal matrix with the eigenvalues of $A-KCA$ as diagonal entries. Denote these $n$ eigenvalues (diagonal entries) as $\pi_{1},\cdots,\pi_{n}$.
We design local estimation $\zeta_{i}(k)$ with the following dynamic:
\begin{equation}\label{eq:def_zeta}
\zeta_{i}(k+1)=\Pi \zeta_{i}(k)+\mathbf{1}_{n} y_{i}(k+1) .
\end{equation}
In other words, $\zeta_{i}(k)$ is the state of a estimator with dynamics similar as $A-KCA$ while only takes observations from sensor $i$.

Define $G_i$ as
\begin{equation}\label{eq:def_Gi}
G_{i} \triangleq\left[\begin{array}{c}
C_{i} A\left(A-\pi_{1} I\right)^{-1} \\
\vdots \\
C_{i} A\left(A-\pi_{n} I\right)^{-1}
\end{array}\right].
\end{equation}
In the following, we claim that local estimate $\zeta_i(k)$ is actually a stable estimate of state $x(k)$ transformed by the constant matrix $G_i$, i.e., their difference $\epsilon_i(k)\triangleq\zeta_{i}(k)-G_ix(k)$ is a stationary Gaussian process.
The following result has been proved in \cite{liuxinghua-TAC2020}.
\begin{lemma}[\hspace{-0.0001pt}\cite{liuxinghua-TAC2020}]\label{lm:epsilon}
	Let $\epsilon_i(k)\triangleq\zeta_{i}(k)-G_ix(k)$, then
	\begin{align}
	\epsilon_{i}(k+1)= \Pi \epsilon_{i}(k)+\left(G_{i}-\mathbf{1}_{n} C_{i}\right) w(k) \notag \\
	-\mathbf{1}_{n} v_{i}(k+1)&-\mathbf{1}_{n} a_{i}(k+1) . \label{eq:epsilon}
	\end{align}
\end{lemma}
Noticing that $\Pi$ is a strictly stable matrix (all eigenvalues are within the open unit disk) and $w(k),v(k)$ are Guaissain zero-mean noise, when the attack is absence, i.e., $a_i(k)=0$, the residue $\epsilon_i(k)$ is a stationary Gaussian process.
Define $\tilde{Q} \in {\mathbb{C}}^{m n \times m n}$ as the covariance of noise term $\left(G_{i}-\mathbf{1}_{n} C_{i}\right) w(k) -\mathbf{1}_{n} v_{i}(k+1)$ for all $i$, i.e.,
\begin{align*}
\tilde{Q} \triangleq&
\cov\left(	\begin{bmatrix}
G_1-\mathbf{1}_{n} C_1 \\
\vdots \\
G_m-\mathbf{1}_{n} C_m
\end{bmatrix} w(k)\right)+
\cov\left(	\begin{bmatrix}
\mathbf{1}_{n} v_{1}(k+1)\\
\vdots \\
\mathbf{1}_{n} v_{m}(k+1)
\end{bmatrix} \right)
\\
=&
\begin{bmatrix}
G_1-\mathbf{1}_{n} C_1 \\
\vdots \\
G_m-\mathbf{1}_{n} C_m
\end{bmatrix}
Q\begin{bmatrix}
G_1-\mathbf{1}_{n} C_1 \\
\vdots \\
G_m-\mathbf{1}_{n} C_m
\end{bmatrix}^{'}
+ R\otimes \mathbf{1}_{n\times n},
\end{align*}
where $\otimes$ is the Kronecker product.
Define $\tilde{\Pi}\in {\mathbb{C}}^{m n \times m n}$ as
$$
\tilde{\Pi} \triangleq\left[\begin{array}{ccc}
\Pi & & \\
& \ddots & \\
& & \Pi
\end{array}\right].
$$
As $k{\rightarrow}{\infty}$, the stable covariance of $\epsilon(k)\triangleq\left[\epsilon_1(k){'},\cdots,\epsilon_m(k){'}\right]{'}$ is the solution $\tilde{W}$ of the following Lyapunov equation:
$$
\tilde{W}=\tilde{\Pi} \tilde{W} \tilde{\Pi}{'}+\tilde{Q}.
$$
The matrix $\tilde{W}$ is well-defined since $\Pi$ is strictly stable. As a result, the secure estimation can be recovered by the solution of the following optimization problem where $\zeta(k)\triangleq\left[\zeta_1(k){'},\cdots,\zeta_m(k){'}\right]{'} $ and $G\triangleq\left[G{'}_1,\cdots,G{'}_m\right]{'} $.
\begin{subequations}\label{pb:old_lasso}
	\begin{align}
	&\underset{\check{x}(k), \mu(k), \nu(k)}{\operatorname{minimize}}\quad \frac{1}{2} \mu(k){'} \tilde{W}^{-1} \mu(k) + \gamma \|\nu(k)\|_1 \label{min:old_lasso} \\
	& \text{ subject to}\quad
	\zeta(k)=
	G \check{x}(k)+\mu(k)+\nu(k). \label{eq:old_lasso}
	\end{align}
\end{subequations}
The parameter $\gamma$ is a non-negative constant chosen by the system operator.
According to \cite{liuxinghua-TAC2020}, the solution $\check{x}(k)$ to problem (\ref{pb:old_lasso}) is a secure estimation and has the following properties.
\begin{theorem}[\hspace{-0.01pt}\cite{liuxinghua-TAC2020} Theorem 4]\label{th:TAC}
	In the presence of $p$-sparse attack, the state estimation $\check{x}(k)$ solved from problem \eqref{pb:old_lasso} is secure if the following inequality holds for all $u \neq \mathbf{0}$, $u\in{\mathbb{R}}^n$:
	\begin{equation}\label{eq:cond}
	\sum_{i \in \mathcal{I}}\left\|G_{i} u\right\|_{1}<\sum_{i \in \mathcal{I}^{c}}\left\|G_{i} u\right\|_{1}, \quad \forall\ {\mathcal{I}}\subset {\mathcal{R}}, |\mathcal{I}|\leq p .
	\end{equation}
\end{theorem}
Even though Theorem \ref{th:TAC} establishes the sufficient condition of the estimation to be secure, validating \eqref{eq:cond} is NP-hard.
In the following section, we reduce the complexity of checking condition \eqref{eq:cond} by transforming $G_i$ to its canonical form under the assumption that the geometric multiplicities of all the eigenvalues of $A$ are $1$ ($A$ is non-derogatory).

\section{Main Results}\label{sec:main_result}
In this section, under the assumption on $A$, we first prove that the span of rows of $G_i$ coincides with the observability space of $(A,C_i)$, which implies that the matrix $G_i$ has a canonical form under row operations. Based on the canonical form, a new secure estimation scheme is proposed.

\subsection{Canonical form of $G_i$}\label{subsec:transform}
In order to prevent degeneration problem, we introduce the following assumption.
\begin{assumption}\label{as:geo_mul_1}
	System matrix $A$ is non-singular, and all the eigenvalues of $A$ have geometric multiplicity $1$.
\end{assumption}
Without loss of generality, based on Assumption \ref{as:geo_mul_1}, we assume that $A$ is in the following Jordan canonical form:
\begin{align*}
A=\begin{bmatrix}
J_{1} & \mathbf{0} & \cdots & \mathbf{0} \\
\mathbf{0} & J_{2} & \cdots & \mathbf{0} \\
\vdots & \vdots & \ddots & \vdots \\
\mathbf{0} & \mathbf{0} & \cdots & J_{l}
\end{bmatrix},\
J_k=
\begin{bmatrix}
\lambda_{k} & 1 & 0 &  \cdots & {0} \\
{0} & \lambda_{k} & 1 & \cdots & {0} \\
\vdots & \ddots & \ddots   & \ddots & \vdots \\
\vdots & \ddots & \ddots  & \ddots & 1 \\
{0} & \cdots & \cdots & 0 & \lambda_{k}
\end{bmatrix},
\end{align*}
and the (complex) eigenvalues with respect to distinct blocks are different, i.e., $k_1\neq k_2$ implies $\lambda_{k_1}\neq\lambda_{k_2}$.
Define the observability matrix of system $(A,C_i)$ as
\begin{equation}\label{eq:def_O}
O_{i} \triangleq\left[\begin{array}{c|c|c|c}
C_{i}{'} &
\left(C_{i} A\right){'} &
\cdots &
\left(C_{i} A^{n-1}\right){'}
\end{array}\right]{'}.
\end{equation}
Before continuing on, we need the following notation of state observability.
Define ${\mathcal{S}}_j$ as the index set of sensors that can observe state $j$, i.e.
\begin{equation}\label{eq:def_Oc}
{\mathcal{S}}_j\triangleq \{i\in{\mathcal{R}}\ |\ O_i{'} e_j\neq \mathbf{0} \},
\end{equation}
where ${\mathcal{R}}\triangleq \{1,2,\cdots,m\}$ is the index set of all sensors and $e_j$ is the $n$-dimensional canonical basis vector with 1 on the $j$-th entry and 0 on the other entries.
We have the following theorem characterizing the structure of $G_i$.
\begin{theorem}\label{th:span}
	Assume system matrix $A$ satisfies Assumption \ref{as:geo_mul_1}, then the following equation holds:
	\begin{align}
	\rs(G_i)=\rs(O_i)=\rs(H_i) ,
	\end{align}
	where $H_i$ is the following diagonal matrix
	\begin{equation*}
	H_i\triangleq \begin{bmatrix}
	{\mathbb{I}}_{i\in{\mathcal{S}}_1} & & & \\
	&{\mathbb{I}}_{i\in{\mathcal{S}}_2} & &  \\
	& & \ddots &  \\
	& & & {\mathbb{I}}_{i\in{\mathcal{S}}_n}
	\end{bmatrix} ,
	\end{equation*}
	and $\mathbb{I}_{i\in{\mathcal{S}}}$ is the indicator function that takes the value 1 when ${i\in{\mathcal{S}}}$ and value 0 when ${i\notin{\mathcal{S}}}$.
\end{theorem}

\begin{proof}
	Define the characteristic polynomial of $A$ as $p(x)=a_n x^n +\cdots+a_1 x +a_0$.
	Define polynomial fraction  $q_\pi(x)$ with respect to constant $\pi$ as
	$q_\pi(x)\triangleq\frac{p(x)-p(\pi)}{x-\pi}$ where $x\neq \pi$.
	For polynomial fraction $q_\pi(x)$, we have
	\begin{align*}
	&\frac{p(x)-p(\pi)}{x-\pi}\\
	=&\frac{1}{x-\pi}\left(a_n(x^n-\pi^n)+a_{n-1}(x^{n-1}-\pi^{n-1})+\cdots+a_1(x-\pi)\right)\\
	=&a_n\left(x^{n-1}+\pi x^{n-2}+\pi^2x^{n-3}+\cdots+\pi^{n-1}\right)\\
	&+a_{n-1}\left(x^{n-2}+\pi x^{n-3}+\pi^2x^{n-4}+\cdots+\pi^{n-2}\right)\\
	&+\cdots \\
	&+a_{2}\left(x+\pi\right)\\
	&+a_1,
	\end{align*}
	where the last inequality comes from the fact that
	$$\frac{x^n-\pi^n}{x-\pi}=\pi^0x^{n-1}+\pi^1 x^{n-2}+\pi^2x^{n-3}+\cdots+\pi^{n-1}x^0.$$
	The above interpretation of $q_\pi(x)$ is also valid when $x$ is a square matrix.
	As a result, by rearranging the terms of $A$ with the same power, one obtains
	\begin{align*}
	q_{\pi_j}(A)=&a_n\left(A^{n-1}+\pi_j A^{n-2}+\cdots+\pi_j^{n-1}I\right)\\
	&+a_{n-1}\left(A^{n-2}+\pi_j A^{n-3}+\cdots+\pi_j^{n-2}I\right)\\
	&+\cdots+a_{2}\left(A+\pi_j I\right)+\cdots+a_1 I\\
	=&a_nA^{n-1}+\left(a_n\pi_j+a_{n-1}\right)A^{n-2}+\cdots\\
	&+\left(a_n\pi_j^{n-1}+a_{n-1}\pi_j^{n-2}+\cdots+a_2\pi_j+a_1\right)I.
	\end{align*}
	Define $b_{j,k}$ as the constant scalar coefficient of $A^k$ in $q_{\pi_j}(A)$:
	\begin{equation}\label{eq:bjk}
	b_{j,k}\triangleq\sum_{i=0}^{n-k-1} a_{i+k+1} \pi_j^i.
	\end{equation}
	Thus, $q_{\pi_j}(A)=\sum_{k=0}^{n-1}b_{j,k}A^{k}$.
	The $j$-th row of matrix $G_i$ can be reformulated as
	\begin{align*}
	&C_{i} A\left(A-\pi_{j} I\right)^{-1}\notag \\
	=&-\frac{1}{p(\pi_j)} C_i A q_{\pi_j}(A)\notag \\
	=&-\frac{1}{p(\pi_j)} C_i \left(\sum_{k=0}^{n-1}b_{j,k}A^{k+1}\right)\notag \\
	=&-\frac{1}{p(\pi_j)}\begin{bmatrix}
	b_{j,0} & b_{j,1} & \cdots  & b_{j,n-1}
	\end{bmatrix}
	\begin{bmatrix}
	C_i A \\
	C_i A^2 \\
	\vdots \\
	C_i A^n
	\end{bmatrix}\notag \\
	=&-\frac{1}{p(\pi_j)}\begin{bmatrix}
	b_{j,0} & b_{j,1} & \cdots  & b_{j,n-1}
	\end{bmatrix} O_i A.
	\end{align*}

	Therefore, $G_i$ can be interpreted as follows
	\begin{align}\label{eq:Gi=DOA}
	G_i = \mathcal{D}_1\begin{bmatrix}
	b_{1,0} & b_{1,1} & \cdots  & b_{1,n-1} \\
	b_{2,0} & b_{2,1} & \cdots  & b_{2,n-1} \\
	\vdots & \vdots & \ddots  & \vdots \\
	b_{n,0} & b_{n,1} & \cdots  & b_{n,n-1}
	\end{bmatrix}
	O_i A
	= & \mathcal{D}_1\mathcal{D}_2\mathcal{D}_3 O_i A ,
	\end{align}
	where $\mathcal{D}_1\triangleq\text{diag}\left(-\frac{1}{p(\pi_1)},-\frac{1}{p(\pi_2)},\cdots,-\frac{1}{p(\pi_n)}\right)$,
	\begin{align*}
	\mathcal{D}_2\triangleq\hspace{-2pt}
	\begin{bmatrix}
	\pi_1^{n-1} & \pi_1^{n-2} & \cdots  & 1 \\
	\pi_2^{n-1} & \pi_2^{n-2} & \cdots  & 1 \\
	\vdots & \vdots & \cdots  & \vdots \\
	\pi_n^{n-1} & \pi_n^{n-2} & \cdots  & 1
	\end{bmatrix},
	\mathcal{D}_3\triangleq\hspace{-2pt}
	\begin{bmatrix}
	a_n & 0 & \cdots &   0 \\
	a_{n-1} & a_n & \cdots &   0 \\
	\vdots & \vdots & \ddots  & \vdots \\
	a_1 & a_2 & \cdots  & a_n
	\end{bmatrix}.
	\end{align*}
	According to Assumption \ref{as:distinct_eigvalue}, all $\pi_j$ are distinct eigenvalues and they are not the eigenvalues of $A$, i.e. the diagonal matrix $\mathcal{D}_1$ and the Vandermonde matrix $\mathcal{D}_2$ are invertible. Moreover, $a_n=1$. Therefore, the lower triangular Toeplitz matrix $\mathcal{D}_3$ is invertible and thus $\rs(G_i)=\rs(O_i A)$ from equation \eqref{eq:Gi=DOA}.
	We continue to prove $\rs(O_i)=\rs(O_i A)$. Considering that $A^n=-a_{n-1}A^{n-1}-\cdots-a_0 I$, one obtains the following equation \eqref{eq:O_OA}.
	\begin{equation}
	\label{eq:O_OA}
	O_i A=
	\begin{bmatrix}
	0 & 1 & 0 &  \cdots & 0 \\
	0 & 0 & 1 &  \cdots & 0 \\
	\vdots & \vdots & \vdots & \ddots & \vdots \\
	0 & 0 & 0 &  \cdots & 1 \\
	-a_0 & -a_1 & -a_2 & \cdots &  -a_{n-1}
	\end{bmatrix}
	O_i .
	\end{equation}
	According to Assumption \ref{as:geo_mul_1}, $A$ is invertible and $a_0=(-1)^n\det(A)\neq 0$, which leads to the equation that $\rs(O_i)=\rs(O_i A)$.
	Since $A$ is assumed to be in the Jordan canonical form and all eigenvalues have geometric multiplicity 1, one can verify that nonzero columns of $O_i$ are linear independent, and equivalently nonzero columns of $G_i$ are linear independent. Therefore, $i\in{\mathcal{S}}_j$ is equivalent to that $j$-th column of $G_i$ is non-zero, i.e., $G_i$ has the same row-span with the canonical form $H_i$.
\end{proof}

Recall that $\rs(\cdot)$ represents the span of row vectors of a matrix.
According to Theorem \ref{th:span} and definition of $\rs(\cdot)$, one directly obtains the following corollary:
\begin{corollary}\label{co:P_i}
	For every sensor index $i\in{\mathcal{R}}$, there exists an invertible matrix $P_i$ such that $P_iG_i=H_i$.
\end{corollary}

After transformation $P_i$, matrix $G_i$ is transformed into canonical form $H_{i}$ whose rows are either canonical basis vectors or zero vectors.
The non-zero entries of $H_i$ records the state observability of sensor $i$. Therefore, the sparse observability index can be directly obtained from diagonal matrices $H_i,i\in\{1,2,\cdots,m\}$ or, equivalently, from sets ${\mathcal{S}}_j,j\in\{1,2,\cdots,n\}$.
\begin{corollary}\label{co:sparse_obs}
	Denote $s$ as the sparse observability index of system \eqref{eq:system}-\eqref{eq:y_i_def}. Then
	\begin{equation*}
	s=\min_{j\in\{1,2,\cdots,n\}} |{\mathcal{S}}_{j}| - 1 .
	\end{equation*}
\end{corollary}
\begin{proof}
	According to the definition of $s$, for arbitrary $\overline{s}$ that satisfies $\overline{s}\geq s+1$, there exists a state index $j^*$ and a sensor index set ${\mathcal{I}}^*$ with $|{\mathcal{I}}^*|=\overline{s}$ such that ${\mathcal{S}}_{j^*} \cap \left({\mathcal{R}}\setminus{\mathcal{I}}^*\right)=\varnothing$.
	As a result, state $j^*$ can not be observed by any sensor in ${\mathcal{R}}\setminus {\mathcal{I}}^*$, i.e.,
	\begin{equation*}
	e_{j^*}\notin \rs(O_i),\ \forall i\in{\mathcal{R}}\setminus {\mathcal{I}}^*.
	\end{equation*}
	and thus system $(A,C_{{\mathcal{R}}\setminus{\mathcal{I}}^*})$ is not observable.
	For arbitrary $\underline{s}$ that satisfies $\underline{s}\leq s$, arbitrary $j$ and arbitrary ${\mathcal{I}}$ with $|{\mathcal{I}}|=\underline{s}$, one obtains ${\mathcal{S}}_{j} \cap \left({\mathcal{R}}\setminus{\mathcal{I}}\right)\neq\varnothing$, which means for all $j$, there exists $i^*\in{\mathcal{R}}\setminus {\mathcal{I}}$ such that: $e_{j}\in \rs(O_{i^*})$. Therefore, system $(A,C_{{\mathcal{R}}\setminus{\mathcal{I}}})$ is observable. According to Definition \ref{df:sparse_obs}, the system is $s$-sparse observable.
\end{proof}

In view of Corollary \ref{co:sparse_obs}, the system sparse observability index can be obtained with low computational complexity since calculating cardinality of ${\mathcal{S}}_j$ can be done in \emph{polynomial time} with respect to $m$ and $n$.
In the following subsection, we will propose an information fusion scheme that is secure in the presence of $p$-sparse attack as long as the system is $2p$-sparse observable.

\subsection{Secure Information Fusion}
Recalling the transformation $P_i$ introduced in Corollary \ref{co:P_i}, define $\tilde{P} \triangleq \text{diag}\left(P_1,\cdots,P_m\right),\ \tilde{M}\triangleq\tilde{P}\tilde{W}\tilde{P}{'}$ and
\begin{align}
{\mathcal{Y}} (k)\triangleq
\begin{bmatrix}
P_1\zeta_{1}(k) \\
\vdots \\
P_m\zeta_{m}(k)
\end{bmatrix}, \
H\triangleq\begin{bmatrix}
H_{1} \\
\vdots \\
H_{m}
\end{bmatrix} .
\end{align}
We present the following optimization problem whose solution is our proposed estimation. The constant scalar $\gamma$ is an adjustable parameter.
\begin{subequations}\label{pb:lasso}
	\begin{align}
	&\underset{\tilde{x}(k), \mu(k),\nu(k)}{\operatorname{minimize}}\quad \frac{1}{2} \mu(k){'} \tilde{M}^{-1} \mu(k) + \gamma\|\nu(k)\|_1  \\
	&\ \text{subject to }\quad
	{\mathcal{Y}}(k)=H \tilde{x}(k) +\mu(k)+\nu(k).
	\end{align}
\end{subequations}

We first prove that, in the absence of attack, the estimation  $\tilde{x}(k)$ is optimal for certain probability that can be adjusted.
We need to define the following matrices to facilitate the understanding of the result.
Define
$$F_i\triangleq V\diag(V^{-1}K_i),\ F=\begin{bmatrix} F_1&F_2&\cdots & F_m \end{bmatrix},$$
where $V$ is defined in \eqref{eq:VLambda} and $K_i$ is the $i$-th column of Kalman gain $K$. Moreover, $\diag(V^{-1}K_i)$ is an $n\times n$ diagonal matrix with the $j$-th diagonal entry equals to $j$-th element of vector $V^{-1}K_i$.
Recalling that $\epsilon(k)\triangleq\left[\epsilon_1(k)',\cdots,\epsilon_m(k)'\right]'$ and $\epsilon_i(k)=\zeta_i(k)-G_ix(k)$ from Lemma \ref{lm:epsilon},
the following theorem characterizes the estimation $\tilde{x}(k)$ in the absence of attack.

\begin{theorem}\label{th:no_attack}
	In the absence of attack, if the parameter $\gamma$ in problem (\ref{pb:lasso}) satisfies
	\begin{align}\label{eq:kalman_cond}
	\left\|\tilde{M}^{-1}\left(I-GF\right)\epsilon(k)\right\|_\infty \leq \gamma,
	\end{align}
	then our proposed estimation $\tilde{x}(k)$ solved from problem \eqref{pb:lasso} is equivalent to the Kalman estimation $\hat{x}(k)$ defined in (\ref{eq:fix_gain_kalman}), i.e.,
	\begin{equation}\label{eq:eq_to_kalman}
	\tilde{x}(k)=\hat{x}(k).
	\end{equation}
\end{theorem}
\begin{remark}
	Noticing that $\epsilon(k)$ is a stationary Gaussian process from Lemma \ref{lm:epsilon}, the probability that \eqref{eq:kalman_cond} holds can be calculated given the parameters $A,C,Q,R$ of the linear Gaussian system \eqref{eq:system}-\eqref{eq:y_i_def}.
	By tuning design parameter $\gamma$, the probability of recovering the Kalman estimation can be adjusted.
\end{remark}
\begin{proof}
	Consider the following least square problem by suppressing $\nu(k)=\mathbf{0}$ in problem \eqref{pb:lasso}.
	\begin{subequations}\label{pb:least_square}
		\begin{align}
		&\underset{\tilde{x}_\text{ls} (k), \varphi(k)}{\operatorname{minimize}}\quad\ \frac{1}{2} \varphi(k)' \tilde{M}^{-1} \varphi(k)   \\
		&\text{subject to }\quad
		{\mathcal{Y}}(k)=H \tilde{x}_\text{ls}(k) +\varphi(k).
		\end{align}
	\end{subequations}
	It can be verified that the solution is equivalent to problem (30) in \cite{liuxinghua-TAC2020}.
	Thus, the solution $\varphi(k)$ of problem \eqref{pb:least_square} satisfies the following equation due to Theorem 1 in \cite{liuxinghua-TAC2020}:
	\begin{equation}\label{eq:ls_recover}
	\tilde{x}_\text{ls}(k)=\hat{x}(k),\ \varphi(k)=(I-GF)\epsilon(k),
	\end{equation}
	where $\hat{x}(k)$ is the fixed gain Kalman estimation defined in \eqref{eq:fix_gain_kalman}.
	Considering the KKT condition of problem \eqref{pb:lasso}, one obtains that if $\left\|\tilde{M}^{-1} \varphi(k)\right\|_{\infty}\leq\gamma$, then the solution $\tilde{x}(k),\mu(k),\nu(k)$ satisfy
	\begin{equation}\label{eq:lasso_recover}
	\tilde{x}(k)=\tilde{x}_\text{ls}(k)=\hat{x}(k),\ \mu(k)=\varphi(k),\ \nu(k)=\mathbf{0}.
	\end{equation}
	Combining \eqref{eq:ls_recover} and \eqref{eq:lasso_recover}, result in Theorem \ref{th:no_attack} is obtained.
\end{proof}

In the presence of $p$-sparse attack,
we have the following theorem demonstrating that the estimation $\tilde{x}(k)$ solved from problem \eqref{pb:lasso} is secure.

\begin{theorem}\label{th:main}
	In the presence of arbitrary $p$-sparse attack, if the system $(A,C)$ is $2p$-sparse observable, the proposed estimation $\tilde{x}(k)$ solved from (\ref{pb:lasso}) is secure.
\end{theorem}

The proof of Theorem \ref{th:main} is provided in Appendix \ref{ap:th_main} for legibility.
Since sparse observability index only requires simple computation according to Corollary \ref{co:sparse_obs}, this work reduces the complexity of evaluating system vulnerability significantly under the assumption of geometric multiplicity.
For general $A$ that has eigenvalues with geometric multiplicity larger than 1 ($A$ is derogatory), computing sparse observability index is an NP-hard problem \cite{ShoukryTAC2016}\cite{mao2021computational}\cite{yanwen_CDC19}, and there is no computational efficient solution unless P$=$NP. Simultaneously, for algorithm online operation, the computing of our proposed estimation involves solving a convex optimization problem based on LASSO \cite{LASSOTibshirani}, which can be done efficiently.

Moreover, the condition of $2p$-sparse observable is necessary for secure state recovering in the presence of $p$ compromised sensors \cite{ShoukryTAC2016}, which establishes the fundamental limit. Our proposed estimator achieves this limit and also holds optimality in the absence of attack for certain probability, which can be adjusted by parameter $\gamma$, according to Theorem \ref{th:no_attack}.
By properly tuning parameter $\gamma$, the system operator can achieve the trade-off between the performance in normal operation and performance under attack.

\section{Simulation}\label{sec:sim}
We use an inverted pendulum for the numerical simulation. The physical parameters are illustrated in Fig. \ref{fig:invpen}.
The control input $u(k)$ is the force applied on the cart, and we assume that there are no frictions of any form.
The states $x_1,x_2,x_3,x_4$ represent cart position coordinate, cart velocity, pendulum angle from vertical and pendulum angle velocity respectively.
\begin{figure}[htpb]
	\centering\vspace{-10pt}
	\input{inv_pen.tex}\vspace{-10pt}
	\caption{Illustration of the inverted pendulum.}\label{fig:invpen}
\end{figure}

Consider the system linearized at $x_3=x_4=0$, and we sample the continuous-time linear system periodically every $0.1$ seconds. The system matrix $A$ can be transformed into the following Jordan canonical form $$
\begin{bmatrix}
1  & 1 &  0 & 0\\
0  & 1 & 0 &  0\\
0  & 0 & 0.642 & 0  \\
0  & 0 & 0 &  1.557
\end{bmatrix},
$$
and there is a Jordan block with size $2\times 2$ and the geometric multiplicity of all the eigenvalues are 1.

The output matrix $C$ is defined as $C_1=C_2=C_3=[1\ 0\ 0\ 0]$,  $C_4=[0\ 0\ 1\ 0]$, i.e., first three sensors are monitoring the cart position $x_1$ and the fourth sensor is monitoring pendulum angle $x_3$.
In this case, the first 3 sensors can observe all the states, and the fourth sensor can only observe state $x_3$ and $x_4$.
According to Corollary \ref{co:sparse_obs}, the sparse observability index is 2, and the system can at most withstand 1 corrupted sensor.

The noise covariances of the system satisfy $w(k) \sim {N}(0, 0.001\cdot I_4)$ and $v(k) \sim {N}(0, 0.001\cdot I_4)$, where $I_n$ is the identity matrix with size $n\times n$.
The controller of the system is designed as a Linear-Quadratic Regulator (LQR), and the feedback matrix is chosen as
$$K_{\rm lqr}=\begin{bmatrix}
-0.604 & -1.678 & -39.514 & -9.721
\end{bmatrix}.
$$
We demonstrate the performance of different estimators of close-loop system where $u(k)=-K_{\rm lqr} x(k)$.

\begin{figure}[htpb]
	\input{states_under_attack_closeloop.tex}
	\caption{Value of system states and their estimation in the presence of attack. The black solid line represents the real state. The red dashed line and blue dotted line represent Kalman estimation and our proposed secure estimation respectively. The initial state is $x(0)=[0,\ 1,\ 0,\ 1]{'}$.} \label{fig:close_loop}

\end{figure}

Fig. \ref{fig:close_loop} presents the Kalman estimation and our proposed secure estimation in the presence of attack on sensor 4. The attack $a_4(k)$ is a time-independent random value uniformly distributed on open interval $\left(-\frac{\pi}{2},\frac{\pi}{2}\right)$. In other words, the measurement of the pendulum angle is randomly shifted by the malicious adversary in magnitude of $\pm 90 ^{\circ}$. The attack is launched since the beginning time. As shown in Fig \ref{fig:close_loop}, Kalman estimation (denoted as red dashed line) has large estimation error under the attack while our estimation (denoted as blue dotted line) converges to the real state quickly.

\begin{figure}[htpb]
	\centering
	\input{MSE.tex}
	\caption{Estimation mean square error (MSE) in the absence and in the presence of attack with varying tuning parameter $\gamma$.} \label{fig:MSE}

\end{figure}

Fig. \ref{fig:MSE} illustrates the estimation mean square error (MSE) of our proposed estimator with varying $\gamma$ in the absence and in the presence of attack on sensor 4.
The attack launched on the system is the same as in Fig. \ref{fig:close_loop}.
The MSE of the oracle Kalman estimation, .i.e., Kalman estimation not affected by the attack, is illustrated by the red dashed line in Fig. \ref{fig:MSE}.
As shown in the figure, by properly choosing $\gamma$, the MSE of our proposed estimator is smaller than that of Kalman estimation (the blue line is below the red horizontal line), with the cost that MSE without attack is slightly larger.
Moreover, the MSE of our proposed estimator without attack coincides with that of Oracle Kalman filter when $\gamma$ is large (e.g., larger than 50), which validates the performance of our estimator in normal operation.

\begin{figure}[htpb]
	\centering
	\input{MSE_mag.tex}
	\caption{Estimation mean square error (MSE) in the absence and in the presence of attack with varying attack magnitude $\|a\|_{\infty}$. The parameter $\gamma$ is set as 5.} \label{fig:MSE_mag}

\end{figure}

Fig. \ref{fig:MSE_mag} illustrates the estimation mean square error (MSE) of our proposed estimator and Kalman filter with varying attack magnitude $\|a\|_{\infty}$. The
attack signal $a_4(k)$ is uniformly distributed in interval $(-\|a\|_{\infty}, \|a\|_{\infty})$. As the magnitude of attack signal increases, MSE of Kalman estimator increases significantly while MSE of our proposed secure estimator remains low. This result demonstrates that our proposed estimator holds low estimation error in the presence of malicious bias data in large magnitude. At the same time, in the absence of attack (when $\|a\|_{\infty}=0$), the estimation error is only slightly larger than that of the optimal Kalman filter.

\section{Conclusion}\label{sec:conclusion}
This paper considers LTI system with Gaussian noise against malicious attack on a subset of sensors.
We improve upon the previous work by reducing the computational complexity in the designing phase while remaining low computational complexity for online estimation, under the assumption that all the geometric multiplicities of eigenvalues of $A$ are $1$ ($A$ is non-derogatory).
To achieve this, we prove that the span of the rows of $G_i$ is equivalent to the observable space associated to sensor $i$, based on which the canonical form $H_i$ is designed. The proposed estimator is formulated as a convex optimization problem based on $H_i$. We further prove that the proposed estimation is secure if the system is $2p$-sparse observable, which is easily checked due to the simple form of $H_i$. Moreover, in the absence of attack, the proposed estimation coincides with Kalman estimation for certain probability, which can be adjusted by tuning parameter $\gamma$ to balance between the performance with and without attack.