text:
\section{Introduction}

Water distribution systems aim to securely provide drinking water to consumers and for firefighting and hence are key infrastructures in our society.
They are subjected to drastic changes in demands on the daily and seasonal timescales \cite{Butler2006book} and to various threats such as contamination \cite{Rasekh2014EnvirModelSoft} and power outage \cite{Gullick2005JWaterSupplyResTechAqua}.
For safe and secure water supply, it is important that a water distribution system is well designed and operated against potential failures.

In risk management in water engineering, the risk is estimated to be equal to the product of the probability of failures and the consequence of the failure
\cite{Jun2007JEnvironEng,Giustolisi2010UrbanWaterJ}. In contrast to risk management,
resilience management is a relatively new framework
to aid us to prepare for potential failures in water distribution systems.
Resilience is a broad concept and its definition varies according to studies in engineering \cite{Hosseini2016ReliabilityEngSystSafety} and in other fields \cite{XueWangYang2018NatHazards}. The concept of resilience is equally wide for water distribution systems; it is construed as the general capacity of a system to resist, absorb, withstand and recover from stresses \cite{Hosseini2016ReliabilityEngSystSafety,Butler2017GlobalChallenges}.
It is a useful complement, as many failures that happen in real-life are unforeseeable and are not targeted in the traditional risk management due to the low estimated probability.

Water distribution systems are composed of pipes combined with other functional structures
such as reservoirs, pumps and valves. Network disfunctions such as pipe failure are an obvious threat to the functionality of the water distribution system. Network science has accumulated knowledge on how network structure shapes the network's capacity to withhold random failures and intentional attacks \cite{Newman2010book,Barabasi2016book}. Therefore, there has been a surge of interest in water engineering community to deploy network analysis to assess the degree of resilience of water distribution networks.
However, to the best of our knowledge, most of the existing proposals of resilience measures for water distribution networks are based on the network structure
\cite{Yazdani2011Chaos,Yazdani2011EnvironModelSoftware,Yazdani2012WaterResourcesRes,ShuangZhangYuan2014PlosOne,Herrera2016WaterResourcesMan,Pundit2016IntJCritInfrastruct,Hwang2017JWaterResourPlannManag,Meng2018WaterRes}, or
steady-state or quasi-steady-state flow rates
\cite{Todini2000UrbanWater,Zhuang2013JHydraulEng,Yannopoulos2013WaterResourManag,Zhao2015MathProbEng,Diao2016WaterRes,Liu2017JWaterResourPlannManag}; as such, the dynamic, transient water flows in the pipe network cannot be represented, which is in fact the key to the resilience performance of a water distribution system. Furthermore, one needs to apply stresses to the system to measure resilience \cite{Meng2018WaterRes}, which cannot be exhaustive of what failure might occur in a system. Therefore, it is necessary to examine transient dynamics to reveal much of resilience features of the system, in addition to static network structure or flow.
In fact, transient dynamics of water flow in pipe networks are not a new issue in water engineering.
Examination of transient flows is also a practical requirement because live water distribution networks are incessantly undergoing changes due to, for example, routine operational adjustments, breakdowns of their elements and human error. In turn, transient water flows are known to cause various problems in water distribution systems such as pump failure, collapse of vapor cavities and compromised water quality. Motivated by these needs, various methods to simulate and understand transient flows have been developed \cite{Islam1998JHydraulEng,Boulos2006book,Todini2011JHydroinfo,Nault2016JHydraulEng,JungKarney2017UrbanWaterJ}. While useful, these developments are mathematical modelling of the transient water dynamics and their numerical simulation methods. These methods do not themselves provide a measure of resilience.

In the present study, we propose an index of resilience via linear stability analysis of a standard set of dynamical equations modelling the steady-state flow in the given water distribution network. The equations are equivalent to a type of a nonlinear electric circuit system. The proposed resilience index is defined in terms of the dominant eigenvalue of the Jacobian matrix around the steady state. Therefore, the index corresponds to the speed at which the transient dynamics decay towards the steady state in response to a small perturbation in the flow rate.
The index is calculated from the structure of the network, the diameter and length of the individual pipes, and other constraints such as the energy level at the reservoir nodes and water demand at individual consumer nodes. We test the proposed index against different resilience measures that we previously proposed based on realistic numerical simulations of water distribution systems.
MATLAB code for calculating the so-called local stability index is available at Github
(\url{https://github.com/naokimas/LSI_water_distribution_net}).

\section{Model}

Consider a connected and undirected network of pipes. By definition, a pipe connects two nodes, where nodes are junctions, reservoirs and households. A dead-end node that is only connected to a single pipe (hence, the node's degree is equal to one) corresponds to a consumer node such as a household. It should be noted that nodes whose degree is larger than one are junction nodes but can also demand water if they are connected to households. Let us denote by $N$ and $M$ the number of nodes and edges, respectively. Each edge has its own diameter and length, assuming a circular shape of the pipe. These properties affect how much water flows through the pipe under a given condition, thus effectively specifying the weight (i.e., capacity) of the edge. However, we consider the network as unweighted network and explicitly model the effect of the pipe diameter and length.

We consider a rigid water column model on pipe networks
under the assumption that flows in a pipe vary slowly in time
\cite{Wylie1983JHydraulEng,Islam1998JHydraulEng,Ghidaoui2005ApplMechRev,Chaudhry2014book}. The dynamical equation representing transient flow through the $m$th pipe ($1\le m\le M$), which connects the $i$th and $j$th nodes, is given by
\begin{equation}
I_m\frac{{\rm d}Q_m}{{\rm d}t} = h_i - h_j - R_m Q_m |Q_m| - u(Q_m).
\label{eq:dyn}
\end{equation}
Nodes $i$ and $j$ depend on edge $m$. However, unless we state otherwise, we use $i$ and $j$ here and in the following text to avoid notational abuse.
In Eq.~\eqref{eq:dyn}, $Q_m$ is the time-dependent flow rate through the $m$th pipe, in the direction from the $i$th node to the $j$th node; $t$ represents the time; $h_i$ represents the total energy at the $i$th node, which depends on the time;
$R_m$ is a constant called the pipe coefficient containing the information on the diameter, length and roughness of pipe $m$; $R_m Q_m |Q_m|$ accounts for the nonlinear head loss across the $m$th pipe; $u(Q_m)$ is the input term which comes from a valve or other structure located on the $m$th pipe; $I_m$ is the pipe inertia constant. The dynamical system given by Eq.~\eqref{eq:dyn} is equivalent to that of a nonlinear resistor-inductor electrical circuit, where $Q_m$ is the current, $h_i$ is the voltage, $I_m$ is the inductance, $R_m$ is a nonlinear resistance and $u(Q_m)$ is other nonlinear elements placed on the $m$th edge.

The pipe coefficient is given by
\begin{equation}
R_m = \frac{8 f_m \ell_m}{g \pi^2 d_m^5},
\label{eq:pipe coefficient}
\end{equation}
where $f_m$ is the dimensionless Darcy-Weisbach friction factor for the $m$th pipe and depends on the flow regime (see below for the formula) and other factors,
$\ell_m$ is the length of the $m$th pipe, $g$ is the gravitational acceleration and equal to $9.81$ m/s${}^2$ and $d_m$ is the diameter of the $m$th pipe
\cite{Boulos2006book}. More generally, if the nonlinear head loss term is given by $R_m Q_m
|Q_m|^{\gamma-1}$, one obtains
$R_m = f_m \ell_m / 2gd_m \mathcal{A}_m^{\gamma}$, where
$\mathcal{A}_m$ is the cross-sectional area of the $m$th pipe, i.e., $\mathcal{A}_m = \pi d_m^2/4$ \cite{Todini1988chapter,Boulos2006book}. The pipe inertia constant is given by $I_m = \ell_m/(g \mathcal{A}_m)$.

There are various formula to approximate the friction factor, $f_m$, as a function of the pipe diameter, $d_m$, the roughness coefficient of the pipe (which is determined by the pipe's material, age and other factors), Reynolds number and so forth. Because the Reynolds number is a function of the flow rate, $Q_m$, the friction factor also depends on $Q_m$. Our index developed in section~\ref{sec:theory} involves differentiation of the right-hand side of Eq~\eqref{eq:dyn} with respect to $Q_m$, which requires the derivative of $f_m$ with respect to $Q_m$. In fact, many formula for $f_m$ can be only applicable to particular flow regimes; the flow regime is
specified by the value of the Reynolds number. Because $f_m$ may depend on pipes, different pipes may be in different flow regimes and some pipes may be situated near the boundary of two flow regimes.
Therefore, we select the formula proposed in Ref.~\cite{Bellos2018JHydraulicEng} that is applicable in all flow regimes including laminar and turbulent flows and gives $f_m$ as a continuous function of the Reynolds number (and hence continuous in $Q_m$). The formula is given by
\begin{equation}
f_m = \left(\frac{64}{\mathcal{R}_m}\right)^a \left(0.75 \ln \frac{\mathcal{R}_m}{5.37}\right)^{2(a-1)b}
\left(0.88 \ln \frac{6.82 d_m}{\epsilon}\right)^{2(a-1)(1-b)}.
\label{eq:Bellos}
\end{equation}
In Eq.~\eqref{eq:Bellos}, the Reynolds number for the $m$th pipe is equal to
\begin{equation}
\mathcal{R}_m = \frac{\left|Q_m\right| d_m}{\mathcal{A}_m v},
\end{equation}
where $v = 1.007\times 10^{-6}$ m${}^2/$s is the kinematic viscosity of water, assuming the temperature of 20 degree celsius.
Furthermore, $\epsilon = 2.591 \cdot 10^{-4}$ m is the roughness coefficient of cast iron
\cite{Rossman2000EPANET2} (also see Ref.~\cite{Boulos2006book} for a similar value), and
\begin{align}
a =& \left[1+\left(\frac{\mathcal{R}_m}{2712}\right)^{8.4}\right]^{-1},\\
b =& \left[ 1 + \left( \frac{\epsilon \mathcal{R}_m}{150 d_m}\right)^{1.8}\right]^{-1}.
\end{align}

\section{\label{sec:theory}Local stability index}

In this section, we carry out the linear stability analysis of the steady state of the nonlinear dynamical system given by Eq.~\eqref{eq:dyn}. The system turns out to be always linearly stable, and we propose a local stability index as the magnitude of the eigenvalue corresponding to the most slowly decaying mode. This task is non-trivial because some nodes act as reservoir of water, other nodes act as consumer of water and the friction factor, $f_m$, depends on the dynamical variable, $Q_m$.

\subsection{Removing redundancy from the dynamical equations}

Assume that the total energy at $N_0$ nodes is fixed. The $N_0$ nodes are typically reservoir whose altitude is higher than typical consumer and junction nodes. Without loss of generality, we assume that $h_{N-N_0+1}$, $h_{N-N_0+2}$, $\ldots$, $h_N$ are fixed.
Then, the set of Eq.~\eqref{eq:dyn} has $N-N_0+M$ unknowns, i.e., $h_i$ ($1\le i\le N-N_0$) and $Q_m$ ($1\le m\le M$), whereas Eq.~\eqref{eq:dyn} only provides $M$ differential equations. The remaining constraints come from the Kirchhoff's current law (i.e. conservation of water mass) at each of the $N-N_0$ nodes whose total energy is not fixed. In other words, we obtain
\begin{equation}
\sum_{j=1; j\neq i}^N A_{ij} Q_{ij} = - \tilde{Q}^{\rm ext}_i\quad (i=1, 2, \ldots, N-N_0),
\label{eq:current law}
\end{equation}
where $Q_{ij}$ is the flow rate through pipe $(i, j)$ from the $i$th to $j$th nodes;
$\tilde{Q}^{\rm ext}_i$ is the external demand (i.e., withdrawal) of water at the $i$th node; $A_{ij} (=A_{ji})$ is the element of the adjacency matrix. In other words, $A_{ij}=1$ if there is a pipe between the $i$th and $j$th nodes. Otherwise, $A_{ij} = 0$.

Equations~\eqref{eq:dyn} and \eqref{eq:current law} imply that
there are $N-N_0+M$ unknowns and the same number of equations for solving the steady state and its linear stability.
Now we erase $h_i$ ($1\le i\le N-N_0$) from Eq.~\eqref{eq:dyn} using Eq.~\eqref{eq:current law} to derive a self-contained set of $M$-dimensional dynamical system. It should be noted that we do not have to erase $h_i$ ($N-N_0+1\le i\le N$) from Eq.~\eqref{eq:dyn} because these $h_i$'s are constant. To erase $h_i$ ($1\le i\le N-N_0$), we differentiate Eq.~\eqref{eq:current law} to obtain
\begin{equation}
\sum_{j=1; j\neq i}^N A_{ij} \frac{{\rm d}Q_{ij}}{{\rm d}t} = 0 \quad (i=1, 2, \ldots, N-N_0).
\label{eq:current law/dt}
\end{equation}
By substituting Eq.~\eqref{eq:dyn} in Eq.~\eqref{eq:current law/dt}, one obtains
\begin{equation}
\sum_{j=1; j\neq i}^N A_{ij}
\frac{h_i - h_j - R_{ij} Q_{ij} |Q_{ij}| - u(Q_{ij})}{I_{ij}} = 0\quad (i=1, 2, \ldots, N-N_0),
\label{eq:current law 2}
\end{equation}
where $I_{ij}$ is the inertia constant of pipe $(i, j)$.
Note that Eq.~\eqref{eq:pipe coefficient} implies that $R_{ij} = R_{ji}$.

We define the $M$-dimensional diagonal matrix $D$ by the diagonal elements $D_{mm} = 1/I_m$ ($1\le m\le M$). We also denote the head vector by $\bm h = (h_1, \ldots, h_N)^{\top}$, where $\top$ denotes the transposition. Furthermore, we denote by $\bm q = (q_1, \ldots, q_M)^{\top}$, of which the $m$th element is given by
\begin{equation}
q_m(Q_m) \equiv R_m Q_m |Q_m| + u(Q_m).
\label{eq:def q_m}
\end{equation}

To rewrite Eq.~\eqref{eq:current law 2}, we introduce the incidence matrix, denoted by $C$, which is an $N\times M$ matrix.
By definition, corresponding to each $m$th edge $(i, j)$, the entries of the $m$th column of $C$ are given by $C_{i,m} = 1$, $C_{j,m} = -1$ and $C_{\ell, m} = 0$ ($1\le\ell\le N, \ell\neq i, j$). Although the edges are undirected, we assign $+1$ and $-1$ to $i$ and $j$, respectively, by arbitrarily selecting either node forming the $m$th edge to which $+1$ is assigned, for later convenience. We represent $C$ in the block form given by
\begin{equation}
C = \begin{pmatrix}
C^{(1)}\\ C^{(2)}
\end{pmatrix},
\end{equation}
where the $(N-N_0)\times M$ matrix $C^{(1)}$ is the part of the incidence matrix corresponding to the nodes whose total energy is not fixed, and
the $N_0\times M$ matrix $C^{(2)}$ is the part of the incidence matrix corresponding to the nodes whose total energy is fixed.

We further denote the $N\times N$ combinatorial Laplacian matrix of the network by
\begin{equation}
L_{ij} = \begin{cases}
\sum_{j^{\prime}=1; j^{\prime}\neq i}^N A_{ij^{\prime}} & (i = j),\\
- A_{ij} & (i\neq j).
\end{cases}
\label{eq:Laplacian}
\end{equation}
Using
\begin{equation}
L = C C^{\top},
\end{equation}
we obtain the following block format of the Laplacian matrix:
\begin{equation}
L =
\begin{pmatrix}
L^{(11)} & L^{(12)}\\
L^{(21)} & L^{(22)}
\end{pmatrix}
=
\begin{pmatrix}
C^{(1)} C^{(1)\top} & C^{(1)} C^{(2)\top}\\
C^{(2)} C^{(1)\top} & C^{(2)} C^{(2)\top}
\end{pmatrix},
\end{equation}
where $L^{(11)}$ is the $(N-N_0) \times (N-N_0)$ matrix corresponding to the nodes whose total energy is not fixed, $L^{(12)}$ is an $(N-N_0) \times N_0$ matrix and so forth.

The combinatorial Laplacian matrix for the network with the same set of edges but with the edge weight given by $D_{mm} = 1/I_m$ is called the conductance matrix in electrical circuit theory
\cite{Dorfler2018ProcIeee}. We denote the conductance matrix by $L_{\rm w}$. The conductance matrix is given by
\begin{equation}
L_{\rm w} =
\begin{pmatrix}
L_{\rm w}^{(11)} & L_{\rm w}^{(12)}\\
L_{\rm w}^{(21)} & L_{\rm w}^{(22)}
\end{pmatrix}
= C D C^{\top}.
\end{equation}

Then, we can rewrite Eq.~\eqref{eq:current law 2} as
\begin{equation}
\begin{pmatrix}
L_{\rm w}^{(11)} & L_{\rm w}^{(12)}
\end{pmatrix}
\bm h - C^{(1)} D \bm q = 0.
\label{eq:current law 2 succinct}
\end{equation}
With the notation
\begin{equation}
\bm h = \begin{pmatrix}
\bm{h^{(1)}}\\
\bm{h^{(2)}}
\end{pmatrix},
\end{equation}
where $\bm{h^{(1)}} = (h_1, \ldots, h_{N-N_0})^{\top}$ and $\bm{h^{(2)}} = (h_{N-N_0+1}, \ldots, h_N)^{\top}$,
Eq.~\eqref{eq:current law 2 succinct} is equivalent to
\begin{equation}
C^{(1)}D C^{(1)\top} \bm{h^{(1)}} = - C^{(1)} D C^{(2)\top} \bm{h^{(2)}} + C^{(1)} D \bm q.
\end{equation}
Matrix $C^{(1)}D C^{(1)\top}$ is so-called the loopy Laplacian matrix
\cite{Dorfler2018ProcIeee}. Because its all eigenvalues are positive if the network is a connected network, which we have assumed, $C^{(1)}D C^{(1)\top}$ has the inverse and one obtains
\begin{equation}
\bm{h^{(1)}} =
\left(C^{(1)}D C^{(1)\top}\right)^{-1}\left( - C^{(1)} D C^{(2)\top} \bm{h^{(2)}} + C^{(1)} D \bm q\right).
\label{eq:h^(1) by h^(2)}
\end{equation}
This procedure is essentially the same as the network Kron reduction \cite{Dorfler2018ProcIeee}.

By substituting Eq.~\eqref{eq:h^(1) by h^(2)} in Eq.~\eqref{eq:dyn}, one obtains
\begin{equation}
\frac{{\rm d}\bm Q}{{\rm d}t} = D C^{\top}
\begin{pmatrix}
\left(C^{(1)}D C^{(1)\top}\right)^{-1} \left( - C^{(1)} D C^{(2)\top} \bm{h^{(2)}} + C^{(1)} D \bm q(\bm Q)\right)\\
\bm{h^{(2)}}
\end{pmatrix}
- D \bm q(\bm Q),
\label{eq:dyn 2}
\end{equation}
where $\bm Q \equiv (Q_1, \ldots, Q_M)^{\top}$. Equation~\eqref{eq:dyn 2} is an $M$-dimensional dynamical system. However, it should be noted that the dynamics are confined on a $(M-N+N_0)$-dimensional hyperplane specified by Eq.~\eqref{eq:current law}, respecting the Kirchhoff's current law.

\subsection{Local stability analysis}

In this section, we derive the eigenequation that determines the conventional local linear stability of the steady state of the nonlinear dynamics given by Eq.~\eqref{eq:dyn}.
We use the combination of Eqs.~\eqref{eq:dyn} and \eqref{eq:current law} for solving the steady state and
Eq.~\eqref{eq:dyn 2} to perform the local stability analysis around the obtained steady state.

The steady state is given by setting the left-hand side of Eq.~\eqref{eq:dyn} to zero. By combining the resulting $M$ constraints with the $N-N_0$ constraints provided by Eq.~\eqref{eq:current law}, one can obtain the steady state, which consists of the $N-N_0+M$ unknowns, i.e., $h_i$ ($1\le i\le N-N_0$) and $Q_m$ ($1\le m\le M$).

To set up the Newton-Raphson iteration scheme (e.g. \cite{Boulos2006book}), we rewrite
Eq.~\eqref{eq:current law} as
\begin{equation}
F_i^{\rm node} (\bm Q) \equiv \sum_{j=1; j\neq i}^N A_{ij} Q_{ij} + \tilde{Q}^{\rm ext}_i = 0
\quad (i=1, 2, \ldots, N-N_0).
\label{eq:current law NP}
\end{equation}
We set the left-hand side of Eq.~\eqref{eq:dyn} to 0 to obtain
\begin{equation}
F_m^{\rm pipe}(\bm{h^{(1)}}, \bm Q) \equiv h_{i_m} - h_{j_m} - R_m Q_m |Q_m| - u(Q_m) = 0,
\label{eq:dyn NP}
\end{equation}
where $i_m$ and $j_m$ are the two nodes connected by the $m$th pipe; we use $i_m$ and $j_m$ instead of $i$ and $j$ to avoid possible notational conflict with Eq.~\eqref{eq:current law NP}.
The Jacobian of the right-hand side of Eqs.~\eqref{eq:current law NP} and \eqref{eq:dyn NP} with $N+M$ unknowns, $h_1$, $\ldots$, $h_{N-N_0}$, $Q_1$, $\ldots$, $Q_M$, is given by
\begin{equation}
J = \begin{pmatrix}
O & C^{(1)}\\
\left(C^{(1)}\right)^{\top} & - \text{diag}(q_m^{\prime}(Q_m))
\end{pmatrix},
\end{equation}
where $\text{diag}(q_m^{\prime}(Q_m))$
is the diagonal matrix whose diagonal elements are given by
\begin{equation}
q_m^{\prime}(Q_m) = 2 R_m \text{sgn}(Q_m)Q_m + u^{\prime}(Q_m).
\end{equation}
The update equation for the Newton-Raphson method is given by
\begin{equation}
\begin{pmatrix}
\bm h^{(1)(n+1)}\\
\bm Q^{(n+1)}
\end{pmatrix}
=
\begin{pmatrix}
\bm h^{(1)(n)}\\
\bm Q^{(n)}
\end{pmatrix}
- J^{-1}
\begin{pmatrix}
\bm F_1^{\rm node}(\bm Q^{(n)})\\
\vdots\\
\bm F_{N-N_0}^{\rm node}(\bm Q^{(n)})\\
\bm F_1^{\rm pipe}(\bm h^{(1)(n)}, \bm Q^{(n)})\\
\vdots\\
\bm F_M^{\rm pipe}(\bm h^{(1)(n)}, \bm Q^{(n)})
\end{pmatrix}.
\end{equation}

Denote the obtained steady state by $(h_1^*, \ldots, h_{N-N_0}^*, Q_1^*, \ldots, Q_M^*)$.
We now carry out the local stability analysis around this steady state.
The Jacobian matrix of the right-hand side of Eq.~\eqref{eq:dyn 2}, denoted by
$J^{\rm dyn} = (J^{\rm dyn}_{m m^{\prime}})$, is given by
\begin{equation}
J^{\rm dyn} = \left[ D C^{(1)\top} \left(C^{(1)}DC^{(1)\top}\right)^{-1} C^{(1)} D -  D \right]
\text{diag}\left(q_m^{\prime}(Q_m^*) \right).
\label{eq:J^dyn}
\end{equation}

\bigskip

\noindent
\begin{proposition}
Matrix $B \equiv D C^{(1)\top} \left(C^{(1)}DC^{(1)\top}\right)^{-1} C^{(1)} D -  D$ has $(N-N_0)$-fold zero eigenvalues and $M-N+N_0$ negative eigenvalues.
\end{proposition}

\bigskip

\begin{proof}
A direct substitution verifies
\begin{equation}
B (C^{(1)})^{\top} = 0.
\end{equation}
Therefore, $B$ has $(N-N_0)$-fold zero eigenvalues, and the corresponding $N-N_0$ right eigenvectors are the columns of $(C^{(1)})^{\top}$.

Because $\left(C^{(1)}DC^{(1)\top}\right)^{-1}$ is an $(N-N_0) \times (N-N_0)$ matrix, its rank is at most $N-N_0$. Therefore, the $M\times M$ matrix $B^{\prime} \equiv C^{(1)\top} \left(C^{(1)}DC^{(1)\top}\right)^{-1} C^{(1)} D$ has at most rank $N-N_0$ and therefore has at least $M-N+N_0$ zero eigenvalues.

We denote the eigenvalues of a symmetric matrix $X$ by $\lambda_1(X) \le \cdots \le \lambda_M(X)$.
Using the fact that $D$ is positive definite because all its diagonal elements are positive (i.e., $= 1/I_m$), the
Weyl's theorem on eigenvalues applied to the sum of $B^{\prime}$ and $-D$, which is equal to $B$, implies that
\begin{equation}
\lambda_m(B) < \lambda_m(B^{\prime}).
\label{eq:B<B'}
\end{equation}
Note that $\lambda_{m_0}(B^{\prime}) = \lambda_{m_0+1}(B^{\prime}) = \cdots = \lambda_{m_0+M-N+N_0-1}(B^{\prime}) = 0$ for an $m_0$ ($1\le m_0\le N-N_0+1$). Therefore, Eq.~\eqref{eq:B<B'} yields
$\lambda_{m_0}(B), \ldots, \lambda_{m_0+M-N+N_0-1}(B) < 0$. Because this relationship must be consistent with the fact that $B$ has $(N-N_0)$-fold zero eigenvalues, we obtain $m_0=1$ and
$\lambda_1(B), \ldots, \lambda_{M-N+N_0}(B) < 0, \lambda_{M-N+N_0+1}(B) = \cdots = \lambda_M(B) = 0$.
\end{proof}

Let us assume for any pipe $m$ that $q_m^{\prime}(Q_m^*) > 0$, which means that the head loss owing to friction along the pipe increases as the steady-state flow rate, $Q_m^*$, increases. We verified $q_m^{\prime}(Q_m^*) > 0$ held true in all the numerical simulations carried out in the next section. Then, $J^{\rm dyn}$ has $M-N+N_0$ negative eigenvalues and $N-N_0$ zero eigenvalues as matrix $B$ does (Appendix~A).

Therefore, the dynamics given by Eq.~\eqref{eq:dyn 2} are locally neutrally stable in the subspace spanned by the column vectors of $\left[\text{diag}\left(q_m^{\prime}(Q_m^*) \right)\right]^{-1} (C^{(1)})^{\top}$. However, perturbations along these directions are infeasible because such a perturbation is not on the hyperplane defined by Eq.~\eqref{eq:current law} on which the dynamics are constrained.
Note that the normal vectors of the hyperplane are given by the row vectors of $C^{(1)}$, i.e., column vectors of $(C^{(1)})^{\top}$.
Therefore, in the local stability analysis, we should ignore these $N-N_0$ zero eigenvalues and inspect the other $M-N+N_0$ eigenvalues. With this consideration, we conclude that the steady state constrained by the Kirchhoff's current law (i.e., Eq.~\eqref{eq:current law}) is always locally stable.

To define a local stability index,
we calculate the largest eigenvalue of $J^{\rm dyn}$ except the $(N-N_0)$-fold zero eigenvalue. Note that all eigenvalues of $J^{\rm dyn}$ are real although $J^{\rm dyn}$ is not a symmetric matrix in general
(Appendix~A).
We consider the
the largest non-zero eigenvalue of $J^{\rm dyn}$ (which is negative) and refer to its absolute value as the local stability index, denoted by $\rho$. A large $\rho$ value implies that dynamics perturbed by a small amount from the steady state would rapidly return to the steady state. Mathematically, the perturbations are given to the flow rate for some pipes. However, we expect that $\rho$ serves as an index that can characterise dynamical stability of water distribution networks in wider contexts.

\section{Numerical results}

In this section, we calculate the local stability index for various networks and examine its relevance to stability and resilience in simulated dynamics of water flow.

We use the synthesized networks used in our previous study \cite{Meng2018WaterRes}. As supplementary materials of Ref.~\cite{Meng2018WaterRes}, the full network structure and properties of system components (e.g., pipes, nodal demands, elevations) were made open to the public in the EPANET2 format, except for the one network. For completeness, we have made the complete data for all the networks and their system components used in the present study available in the MATLAB format
at Github, alongside the code for calculating the local stability index.

Of the 85 networks used in our previous study \cite{Meng2018WaterRes}, five networks include pumps, which would pin the flow rate to a constant value. Because our formalism does not directly cover the case of the pinned flow rate, we exclude these five networks. Of the remaining 80 networks, 10 networks have $N=102$ nodes, $N_0=2$ reservoir nodes among the $N$ nodes and $M=110$ pipes; 11 networks have $(N, N_0, M) = (102, 2, 130)$; 10 networks have $(N, N_0, M) = (204, 4, 223)$; 10 networks have $(N, N_0, M) = (204, 4, 263)$; 9 networks have $(N, N_0, M) = (306, 6, 335)$; 9 networks have $(N, N_0, M) = (306, 6, 395)$; 5 networks have $(N, N_0, M) = (404, 4, 443)$; 5 networks have $(N, N_0, M) = (404, 4, 523)$; 5 networks have $(N, N_0, M) = (406, 6, 525)$; 5 networks have $(N, N_0, M) = (408, 8, 447)$; 1 network has $(N, N_0, M) = (506, 6, 554)$.
We calculated $\rho$ for each of the 80 networks.

The diameter of a majority of pipes was equal to 400 mm. There were also a few larger values of the diameter (e.g., 900 mm). Such wider pipes were incident (i.e., directly connected) to a reservoir.
The diameter of the pipes was automatically generated by the software \cite{Decorte2014WaterResourcesMan}.

We set the total energy of all the $N_0$ reservoirs in each network to 65 m \cite{Decorte2014WaterResourcesMan}. In addition to these $N_0$ nodes, some nodes in each network had zero water demand (i.e., $\tilde{Q}_i^{\rm ext} = 0$) because they are pure junctions connecting different pipes.
We set $u(Q_m)$ for each pipe to zero because the data are from the 80 networks that do not have valves or pumps.

In our previous study, we numerically simulated dynamics of water flow in these networks using
software EPANET2 \cite{Mugume2015WaterRes,Diao2016WaterRes} and measured six strain indices to characterise resilience of water distribution systems. These indices are as follows
(see \cite{Meng2018WaterRes} for fuller definitions). First, the time to strain is defined as the time between the application of stress and the start of service failure, which is when the level of service at a node drops below a predefined threshold. Second, the failure duration is defined as the time needed for the system to recover to normal performance since the service failure triggered by the application of stress. Third, the failure magnitude is defined as the most severe drop at a node in terms of the system service performance as a result of the administered stress, averaged over all the nodes. Fourth, the failure rate is defined as the failure magnitude divided by the time between the start of the failure and the occurrence of the worst system performance. Fifth, the recovery rate is defined as the failure magnitude divided by the time between the worst system performance and the return to the performance threshold value used in the definition of the time to strain. Sixth, the severity is defined as the threshold minus the system performance (which is positive during the system failure) integrated with respect to time over the period of system failure.

The relationship between each of the six strain indices and the local stability index, $\rho$, is shown in
Fig.~\ref{fig:vs local stability index}. In the figure, each circle represents one of the 80 networks,
$r$ represents the Pearson correlation coefficient and $p$ represents the $p$-value for the Pearson correlation coefficient. We found a positive correlation between $\rho$ and the recovery rate (Fig.~\ref{fig:vs local stability index}(e)). This result is intuitive
because the recovery rate quantifies how the system returns to the normal performance after being perturbed by external stress, which agrees with the concept of linear stability when the given external stress is infinitesimally small.
The obtained correlation is only moderate (i.e., $r=0.420$, $p\approx 1.04\times 10^{-4}$) probably due to various nonlinearity of the system, which our linear stability analysis does not aim to explain. However, we consider that $\rho$ captures a strain index derived from EPANET2, which is an involved numerical simulation package, reasonably well.

Figure~\ref{fig:vs local stability index} also indicates that
$\rho$ is strongly correlated with the failure magnitude (Fig.~\ref{fig:vs local stability index}(c)) and the failure rate (Fig.~\ref{fig:vs local stability index}(d)).
We do not have an active interpretation of this result.
In fact, it apparently sounds contradictory that a large $\rho$ value, which predicts that the
network is resilient, is associated with a large failure magnitude and failure rate.
However, the failure magnitude and failure rate quantify the magnitude of failure response induced by external stress, which does not have to do with our local stability analysis even in the limit of infinitesimally small external stress.
Therefore, we conclude that the failure magnitude and the failure rate characterise different aspects of the system's resilience from what the recovery rate and $\rho$ do.
This view is consistent with our additional observation that,
within our previous numerical results \cite{Meng2018WaterRes}, there is positive correlation between the failure magnitude and the recovery rate
($r=0.917$, $p<10^{-32}$)
and between the failure rate and the recovery rate
($r=0.677$, $p\approx 5.60\times 10^{-12}$).

In our previous study, we measured eight structural properties of the water distribution network and examined the association between each of them and each of the six strain indices \cite{Meng2018WaterRes}.
Here we measured the Pearson correlation between each of these structural properties, whose definitions are in
Appendix~B, and the recovery rate. We found that none of the eight structural properties of the network was as strongly correlated with the recovery rate as $\rho$ was (link density: $r = 0.145$, $p = 0.204$; algebraic connectivity: $r = 0.114$, $p = 0.316$;
diameter: $r = -0.143$, $p = 0.207$;
average path length: $r = -0.204$, $p = 0.0717$;
central point dominance: $r = -0.0515$, $p = 0.652$;
heterogeneity: $r = 0.335$, $p = 0.0026$;
spectral gap: $r = -0.0509$, $p = 0.656$;
clustering coefficient: $r = 0.198$, $p = 0.0798$).
Therefore, the local stability index is better at capturing the recovery rate than these structural properties. This is probably because these structural properties do not directly quantify the speed at which the water distribution network responds to a perturbation to the steady state.

\section{Discussion}

We carried out a local stability analysis of water flow in pipe networks including reservoir nodes and consumer nodes. The steady state was shown to be always linearly stable. We used the eigenvalue corresponding to the slowest relaxation mode as the local stability index. The proposed index was moderately correlated with the recovery rate in response to external stress, which was numerically obtained from an involved simulator commonly used in water engineering research community.

We used a particular formula for the friction factor
\cite{Bellos2018JHydraulicEng}. However, the present framework accepts other formulae for the friction factor. Because the differentiation of the friction factor with respect to the Reynolds number is used in calculating the local stability index and the value of the friction factor depends on the pipe in general, one should use a formula that is continuous and covers a wide range of the Reynolds number and hence different flow schemes. Some possible choices of the formula different from the one used in the present paper are found in the literature \cite{Swamee1993JTransportEng,Cheng2008JHydraulicEng,Brkic2018ApplSci}. It is straightforward to extend our local stability index to the case of different functional forms of the friction factor.

In 1972, May proposed to examine the eigenvalue spectra of interaction networks to characterise dynamical stability of complex systems, in particular ecological systems \cite{May1972Nature-stability}. This is
a landmark example of local stability analysis.
While the original analysis and many studies that ensued were for random networks \cite{May1972Nature-stability}, the understanding has been extended to the case of interaction networks with general network structure \cite{Okuyama2008EcolLett,Allesina2012Nature}. This family of method is similar to the one proposed in the present study. Main differences are that we have started from a particular set of hydraulic equations modelling water flow in pipe networks and that we have found that our system is always locally stable and hence used the relaxation mode to quantify the extent of the local stability of the system.

The proposed local stability index was significantly correlated with the recovery rate, which was obtained from numerical simulations in which the applied shock was not necessarily small. However,
because the present analysis is a local stability analysis, it does not tell us the behaviour of the system when a shock injected to the steady state is not small. In water engineering applications, the magnitude of perturbation is not necessarily small \cite{Gullick2005JWaterSupplyResTechAqua}. A recent study in network science has developed a formalism that reduces a class of dynamical systems on interaction networks
into a one-dimensional effective dynamical system and assesses the resilience of the dynamical system against perturbation that is not necessarily small \cite{GaoBarzelBarabasi2016Nature}.
Therefore, adapting their resilience formalism to the case of transient and steady-state dynamics of water flow in pipe networks is of a practical relevance. In Eq.~\eqref{eq:dyn}, there are usually only a fraction of nodes whose energy value is fixed (i.e., reservoirs) and different nodes have different water demands. In this sense, the flow dynamics on pipe networks are heterogeneous in addition to being heterogeneous in the network structure and edge weight. Other structures such as valves would add more complexity and heterogeneity to the system. For such a heterogeneous system, it is unclear whether one can transform Eq.~\eqref{eq:dyn} into a form that is compatible with the resilience function formalism \cite{GaoBarzelBarabasi2016Nature,Tu2017PhysRevE}. This issue is left as future work.

\enlargethispage{20pt}

\section*{Data accessibility}

The data on the network structure and system components are available on Github.

\section*{Competing interests}
The authors declare that they have no competing interests.

\section*{Authors' contributions}

NM and FM conceived of and designed the study. FM provided the data on the network structure and system components. NM carried out the theoretical analysis and performed the data analysis.  NM and FM discussed the results and drafted the manuscript. All authors read
and approved the manuscript.

\section*{Funding}
N.M. and F.M. acknowledge the support provided through
BRIM (Building Resilience into Risk Management) project, EPSRC (EP/N010329/1).

\section*{\label{sec:stability of J^dyn}Appendix A: Eigenvalues of $J^{\rm dyn}$}

Here we provide an elementary proof that $J^{\rm dyn}$ has $M-N+N_0$ zero eigenvalues and $N-N_0$ negative eigenvalues when $q_m^{\prime}(Q_m^*) > 0$, $1\le m\le M$.

Let us introduce a short-hand notation $\overline{D} \equiv \text{diag}\left(q_m^{\prime}(Q_m^*) \right)$. Then, one obtains $J^{\rm dyn} = B \overline{D}$. Because $J^{\rm dyn}$ is similar to the symmetric matrix $\overline{D}^{1/2} B \overline{D}^{1/2}$, the min-max theorem
for the $k$th smallest eigenvalue implies that
\begin{align}
\lambda_k (J^{\rm dyn}) = \lambda_k \left(\overline{D}^{1/2} B \overline{D}^{1/2}\right)
=& \min_{\substack{S\subset \mathbb{R}^M\\ \dim (S) = k}}
\left[ \max_{x\in S\backslash \{0\}} \frac{\left(\overline{D}^{1/2} B \overline{D}^{1/2}x, x\right)}{(x, x)}\right]\notag\\
=& \min_{\substack{S\subset \mathbb{R}^M\\ \dim (S) = k}}
\left[ \max_{x\in S\backslash \{0\}} \frac{\left(\overline{D}^{1/2} B \overline{D}^{1/2}x, x\right)}{\left(\overline{D}^{1/2}x, \overline{D}^{1/2}x\right)}
\cdot \frac{\left(\overline{D}x, x\right)}{(x, x)}\right]\notag\\
\le & \min_{\substack{S\subset \mathbb{R}^M\\ \dim (S) = k}}
\left[ \max_{x\in S\backslash \{0\}} \frac{\left(\overline{D}^{1/2} B \overline{D}^{1/2}x, x\right)}{\left(\overline{D}^{1/2}x, \overline{D}^{1/2}x\right)}
\cdot
\max_{x^{\prime}\in \mathbb{R}^M\backslash \{0\}}
\frac{\left(\overline{D}x^{\prime}, x^{\prime}\right)}{(x^{\prime}, x^{\prime})}\right]\notag\\
=& \lambda_k(B) \lambda_M\left(\overline{D}\right).
\end{align}
Similarly, one obtains
\begin{align}
\lambda_k (J^{\rm dyn})
=& \min_{\substack{S\subset \mathbb{R}^M\\ \dim (S) = k}}
\left[ \max_{x\in S\backslash \{0\}} \frac{\left(\overline{D}^{1/2} B \overline{D}^{1/2}x, x\right)}{\left(\overline{D}^{1/2}x, \overline{D}^{1/2}x\right)}
\cdot \frac{\left(\overline{D}x, x\right)}{(x, x)}\right]\notag\\
\ge & \min_{\substack{S\subset \mathbb{R}^M\\ \dim (S) = k}}
\left[ \max_{x\in S\backslash \{0\}} \frac{\left(\overline{D}^{1/2} B \overline{D}^{1/2}x, x\right)}{\left(\overline{D}^{1/2}x, \overline{D}^{1/2}x\right)}
\cdot
\min_{x^{\prime}\in \mathbb{R}^M\backslash \{0\}}
\frac{\left(\overline{D}x^{\prime}, x^{\prime}\right)}{(x^{\prime}, x^{\prime})}\right]\notag\\
=& \lambda_k(B) \lambda_1\left(\overline{D}\right).
\end{align}
By combining $\lambda_1\left(\overline{D}\right) > 0$, $\lambda_M\left(\overline{D}\right) > 0$ and
$\lambda_1(B), \ldots, \lambda_{M-N+N_0} < 0, \lambda_{M-N+N_0+1} = \cdots = \lambda_M = 0$, one obtains
$\lambda_1(J^{\rm dyn}), \ldots, \lambda_{M-N+N_0}(J^{\rm dyn}) < 0, \lambda_{M-N+N_0+1}(J^{\rm dyn}) = \cdots = \lambda_M(J^{\rm dyn}) = 0$.

\section*{Appendix B: Structural properties of networks}

The eight structural properties of the network that we measured in our previous study \cite{Meng2018WaterRes} and used in the present article for comparison purposes are as follows.

The link density is given by $2M/\left[ N(N-1)\right]$. The algebraic connectivity is the smallest positive eigenvalue of the Laplacian matrix, $L$, of the undirected and unweighted network. The version of the clustering coefficient used in
\cite{Meng2018WaterRes} is
given by $3\times (\text{number of triangles in the network}) / (\text{number of connected triples of nodes in the network})$. The average path length is equal to the shortest path length between two nodes, which is averaged over all the possible $N(N-1)/2$ node pairs. The central point dominance is given by
$\sum_{i=1}^N (b_{\max}-b_i)/(N-1)$, where $b_i$ is the betweenness centrality of the $i$th node and $b_{\max} = \max_{j=1, \ldots, N} b_j$. The heterogeneity is that in the degree and given by
$\sum_{i=1}^N \left( k_i - \langle k\rangle\right)^2 / \langle k\rangle$, where $k_i$ is the degree of the $i$th node and $\langle k\rangle = \sum_{i=1}^N k_i/N$ is the average degree of the node. The version of the spectral gap used in \cite{Meng2018WaterRes} is given by the difference between the two largest eigenvalues of $L$.
The modularity is a quantity representing the quality of community structure detected in the network and approximate modularity maximisation was carried out using Newman's algorithm
\cite{Newman2004PhysRevE-fast}.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
\vspace{-1em}
Deep networks, despite their high predictive accuracy, are notoriously vulnerable to adversarial attacks \citep{goodfellow2014explaining,biggio2013evasion,szegedy2013intriguing,papernot2016limitations}. While many defense methods have been proposed to increase a model's \textit{robustness} to adversarial examples, they were typically observed to hamper its \textit{accuracy} on original clean images. \cite{tsipras2019robustness} first pointed out the inherent tension between the goals of adversarial robustness and standard accuracy in deep networks, whose provable existence was shown in a simplified setting. \cite{zhang2019theoretically} theoretically quantified the accuracy-robustness trade-off, in terms of the gap between the risk for adversarial examples versus the risk for non-adversarial examples.

It is intriguing to consider whether and why the model accuracy and robustness have to be at odds. \cite{schmidt2018adversarially} demonstrated that the
number of samples needed to achieve adversarially robust generalization is polynomially larger than that needed for standard generalization, under the adversarial training setting. A similar conclusion was concurred by \cite{sun2019towards} in the standard training setting. \cite{tsipras2019robustness} considered the accuracy-robustness trade-off as an inherent trait of the data distribution itself, indicating that this phenomenon persists even in the limit of
infinite data. \cite{nakkiran2019adversarial} argued from a different perspective, that the complexity (e.g. capacity) of a robust classifier must be higher than that of a standard classifier. Therefore, replacing a larger-capacity classifier might effectively alleviate the trade-off. Overall, those existing works appear to suggest that, while accuracy and robustness are likely to trade off for a fixed classification model and on a given dataset, such trade-off might be effectively alleviated (``win-win''), if supplying more training data and/or replacing a larger-capacity classifier.

On a separate note, deep networks also face the pressing challenge to be deployed on resource-constrained platforms due to the prosperity of smart Internet-of-Things (IoT) devices. Many IoT applications naturally demand security and trustworthiness, e.g.,\ , biometrics and identity verification, but can only afford limited latency, memory and energy budget. Hereby we extend the question: \textit{can we achieve a triple-win, i.e.,\ , an accurate and robust classfier while keeping it efficient?}

This paper makes an attempt in providing a positive answer to the above question. Rather than proposing a specific design of robust light-weight models, we reduce the average computation loads by input-adaptive routing to achieve triple-win. To this end, we introduce the input-adaptive dynamic inference  \citep{teerapittayanon2016branchynet,wang2017skipnet}, an emerging efficient inference scheme in contrast to the (non-adaptive) model compression, to the adversarial defense field for the first time. Given any deep network backbone (e.g.,\ , ResNet, MobileNet), we first follow \citep{teerapittayanon2016branchynet} to augment it with multiple early-branch output layers in addition to the original final output. Each input, regardless of clean or adversarial samples, adaptively chooses which output layer to take for its own prediction. Therefore, a large portion of input inferences can be terminated early when the samples can already be inferred with high
confidence.

Up to our best knowledge, no existing work studied adversarial attacks and defenses for an adaptive multi-output model, as the multiple sources of losses provide much larger flexibility to compose attacks (and therefore defenses), compared to the typical single-loss backbone. We present a systematical exploration on how to (white-box) attack and defense our proposed multi-output network with adaptive inference, demonstrating that the composition of multiple-loss information is critical in making the attack/defense strong. Fig. \ref{fig:branchy} illustrates our proposed \textit{Robust Dynamic Inference Networks} (RDI-Nets). We show experimentally that the input-adaptive inference and multi-loss flexibility can be our friend in achieving the desired ``triple wins''. With our best defended RDI-Nets, we achieve better accuracy and robustness, yet with over 30\% inference computational savings, compared to the defended original models as well as existing solutions co-designing robustness and efficiency \citep{gui2019model,guo2018sparse}. The codes can be referenced from \changeurlcolor{blue}\href{https://github.com/TAMU-VITA/triple-wins}{\textit{https://github.com/TAMU-VITA/triple-wins.}}

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.8\linewidth]{figures/branchynetv2.pdf}
    \vspace{-1em}
    \caption{Our proposed RDI-Net framework, a \textbf{defended} multi-output network enabling dynamic inference. Each image, being it clean or adversarially perturbed, adaptively picks one branch to exit.}
    \label{fig:branchy}
    \vspace{-1em}
\end{figure}

\section{Related Work}
\vspace{-0.5em}
\subsection{Adversarial Defense}
\vspace{-0.5em}
A magnitude of defend approaches have been proposed \citep{kurakin2016adversarial,xu2017feature,song2017pixeldefend,liao2018defense}, although many were quickly evaded by new attacks \citep{carlini2017towards,baluja2017adversarial}. One
strong defense algorithm that has so far not been fully compromised is adversarial training \citep{madry2017towards}. It searches for adversarial images to augment the training procedure, although at the price of higher training costs (but not affecting inference efficiency). However, almost all existing attacks and defenses focus on a single-output classification (or other task) model. We are unaware of prior studies directly addressing attacks/defenses to more complicated networks with multiple possible outputs.

One related row of works are to exploit model ensemble \citep{tramer2017ensemble,strauss2017ensemble} in adversarial training. The gains of the defended ensemble compared to a single model could be viewed as the benefits of either the benefits of diversity (generating stronger and more transferable perturbations), or the increasing model capacity (consider the ensembled multiple models as a compound one). Unfortunately, ensemble methods could amplify the inference complexity and be detrimental for efficiency. Besides, it is also known that injecting randomization at inference time helps mitigate adversarial effects \citep{xie2017mitigating,cohen2019certified}. Yet up to our best knowledge, no work has studied non-random, but rather input-dependent inference for defense.

\vspace{-0.5em}
\subsection{Efficient Inference}
\vspace{-0.5em}

Research in improving deep network efficiency could be categorized into two streams: the \textit{static} way that designs compact models or compresses heavy models, while the compact/compressed models remain fixed for all inputs at inference; and the \textit{dynamic} way, that at inference the inputs can choose different computational paths adaptively, and the simpler inputs usually take less computation to make predictions. We briefly review the literature below.

\noindent \underline{\textit{Static: Compact Network Design and Model Compression.}} Many compact architectures have been specifically designed for resource-constrained applications, by adopting lightweight depthwise convolutions \citep{sandler2018mobilenetv2}, and group-wise convolutions with channel-shuffling \citep{ZhangZLS2018shufflenet}, to just name a few. For model compression, \cite{han2015learning} first proposed to sparsify deep models by removing non-significant synapses and then re-training to restore performance. Structured pruning was later on introduced for more hardware friendliness \citep{wen2016learning}. Layer factorization \citep{tai2015convolutional,yu2017compressing}, quantization \citep{wu2016quantized}, model distillation~\citep{wang2018adversarial} and weight sharing \citep{wu2018deep} have also been respectively found effective.

\noindent \underline{\textit{Dynamic: Input-Adaptive Inference.}} Higher inference efficiency could be also accomplished by enabling input-conditional execution. \cite{teerapittayanon2016branchynet,huang2018multiscale,kaya2019shallow} leveraged intermediate features to augment
multiple side branch classifiers to enable early predictions. Their methodology sets up the foundation for our work. Other efforts \citep{figurnov2017spatially,wang2017skipnet,wang2018energynet,wang2019dual} allow for an input to choose between passing through or skipping each layer. The approach could be integrated with RDI-Nets too, which we leave as future work.

\vspace{-0.5em}
\subsection{Bridging Robustness with Efficiency}
\vspace{-0.5em}

A few studies  recently try to link deep learning robustness and efficiency. \cite{guo2018sparse} observed that in a sparse deep network, appropriately sparsified weights improve robustness, whereas over-sparsification (e.g., less than 5\% nonzero weights) in turn makes the model more fragile. Two latest works \citep{ye2019adversarial,gui2019model} examined the robustness of compressed models, and concluded similar observations that the relationship between mode size and robustness depends on compression methods and are often non-monotonic. \cite{lin2019defensive} found that activation quantization may hurt robustness, but can be turned into effective defense if enforcing continuity constraints.

Different from above methods that tackle robustness from \textit{static} compact/compressed models, the proposed RDI-Nets are the first to address robustness from the \textit{dynamic} input-adaptive inference. Our experiment results demonstrate the consistent superiority of RDI-Nets over those static methods (Section 4.3). Moreover, applying \textit{dynamic} inference top of those static methods may further boost the robustness and efficiency, which we leave as future work.

\section{Approach}
\vspace{-0.5em}
With the goal of achieving inference efficiency, we first look at the setting of multi-output networks and the specific design of RDI-Net in Section \ref{sec:3-1}. Then we define three forms of adversarial attacks for multi-output networks in Section \ref{sec:3-2} and their corresponding defense methods in Section \ref{sec:3-3}.

Note that RDI-Nets achieve ``triple wins''via reducing the average computation loads through input-adaptive routing. It is not to be confused with any specifically-designed robust light-weight model.

\vspace{-0.5em}
\subsection{Designing RDI-Nets for Higher Inference Efficiency}\label{sec:3-1}
\vspace{-0.5em}
Given an input image $x$, an $N$-output network can produce a set of predictions $[\hat{y_{1}}, ..., \hat{y_{N}}]$ by a set of transformations $[f_{\theta_{1}}(\cdot), ..., f_{\theta_{N}}(\cdot)]$. $\theta_{i}$ denote the model parameter of $f_{\theta_{i}}$, $i = 1, ..., N$, and $f_{\theta_{i}}$s will typically share some weights. With an input $x$, one can express $\hat{y_{i}}=f_{\theta_{i}}(x)$. We assume that the final prediction will be \underline{one} \textbf{chosen} (NOT fused) from $[\hat{y_{1}}, ..., \hat{y_{N}}]$ via some \underline{deterministic} strategy.

We now look at RDI-Nets as a specific instance of multi-output networks, specifically designed for the goal of more efficient, input-adaptive inference. As shown in Fig. \ref{fig:branchy}, \textbf{for any deep network} (e.g.,\ , ResNet, MobileNet), we could append $K$ side branches (with negligible overhead) to allow for early-exit predictions. In other words, it becomes a $(K+1)$-output network, and the \textit{subneworks} with the $K+1$ exits, from the lowest to the highest (the original final output), correspond to $[f_{\theta_{1}}(\cdot), ..., f_{\theta_{K+1}}(\cdot)]$. They share their weights in a \textit{nested} fashion: $\theta_{1} \subseteq \theta_{2} ... \subseteq \theta_{K+1}$, with $\theta_{K+1}$ including the entire network's parameters.

Our deterministic strategy in selecting one final output follows \citep{teerapittayanon2016branchynet}. We set a confidence threshold $t_k$ for each $k$-th exit, $k = 1, ..., K+1$, and each input $x$ will terminate inference and output its prediction in the earliest exit (smallest $k$), whose softmax entropy (as a confidence measure) falls below $t_k$. All computations after the $k$-th exit will not be activated for this $x$. Such a progressive and early-halting mechanism effectively saves unnecessary computation for most easier-to-classify samples, and applies in both training and inference. Note that, if efficiency is not the concern, instead of choosing (the earliest one), we could have designed an adaptive or randomized fusion of all $f_{\theta_{i}}$ predictions: but that falls beyond the goal of this work.

The training objective for RDI-Nets could be written as
\begin{equation}
L_{RDI} = \sum_{i=1}^{K+1} w_{i}[(\phi(f_{i}(\theta_{i}|x), y) + \phi(f_{i}(\theta_{i}|x^{adv}), y)], \label{eqn:weighted-loss}
\end{equation}
For each exit loss, we minimize a hybrid loss of accuracy (on clean $x$) and robustness (on $x^{adv}$). The $K+1$ exits are balanced with a group of weights $\{w_{i}\}_{i=1}^{K+1}$. More details about RDI-Net structures, hyperparameters, and inference branch selection can be founded in Appendix \ref{sec:6-1}, \ref{sec:6-2}, and \ref{sec:6-3}.

In what follows, we discuss three ways to generate $x^{adv}$ in RDI-Nets, and then their defenses.

\vspace{-0.5em}
\subsection{Three Attack Forms on Multi-Output Networks}\label{sec:3-2}
\vspace{-0.5em}
We consider white box attacks in this paper. Attackers have access to the model's parameters, and aim to generate an adversarial image $x^{adv}$ to fool the model by perturbing an input $x$ within a given magnitude bound.

We next discuss three \textit{attack forms} for an $N$-output network. Note that they are independent of, and to be distinguished from \textit{attacker algorithms} (e.g.,\ , PGD, C\&W, FGSM): the former depicts the \textit{optimization formulation}, that can be solved any of the attacker algorithms.

\vspace{-0.5em}
\paragraph{Single Attack}
Naively extending from attacking single-output networks, a single attack is defined to maximally fool one $f_{\theta_{i}}(\cdot)$ only, expressed as:
\begin{equation}
\hspace{-0.01in}
x_{i}^{adv} = \argmax_{x'\in |x'-x|_{\infty}\leq\epsilon}|\phi(f_{\theta_{i}}(x'), y)|,\label{eqn:single-attack}
\end{equation}
where $y$ is the ground truth label, and $\phi$ is the loss for $f_{\theta_{i}}$ (we assume softmax for all). $\epsilon$ is the perturbation radius and we adopt $\ell_\infty$ ball for an empirically strong attacker. Naturally, an $N$-output network can have $N$ different single attacks. However, each single attack is derived without being aware of other parallel outputs. The found $x_{i}^{adv}$ is not necessarily transferable to other $f_{\theta_{j}}$s ($j \neq i$), and therefore can be easily bypassed if $x$ is re-routed through other outputs to make its prediction.

\vspace{-0.5em}
\paragraph{Average Attack}
Our second attack maximizes the average of all $f_{\theta_{i}}$ losses, so that the found $x^{adv}$ remains in effect no matter which one $f_{\theta_{i}}$ is chosen to output the prediction for $x$:
\begin{equation}
\hspace{-0.01in} x_{avg}^{adv} = \argmax_{x'\in |x'-x|_{\infty}\leq\epsilon}|\frac{1}{N}\sum_{j=1}^{N}\phi(f_{\theta_{j}}(x'), y)|,\label{eqn:avg-attack}
\end{equation}
The average attack addresses takes into account the attack transferablity and involves all $\theta_{j}$s into optimization. However, while only one output will be selected for each sample at inference, the average strategy might weaken the individual defense strength of each $f_{\theta_{i}}$.

\vspace{-0.5em}
\paragraph{Max-Average Attack} Our third attack aims to emphasize individual output defense strength, more than simply maximizing an all-averaged loss. We first solve the $N$ single attacks $x^{adv}_{i}$ as described in Eqn. \ref{eqn:single-attack}, and denote their collection as $\Omega$. We then solve the max-average attack via the following:
\begin{equation}
\hspace{-0.01in} x_{max}^{adv} \leftarrow x_{i^*}^{adv} \mbox{, where } x^{adv}_{i^{*}}\in\Omega\mbox{ and } i^* = \arg\max_{i}|\frac{1}{N}\sum_{j=1}^{N}\phi(f_{\theta_{j}}(x_{i}^{adv}), y)|.\label{eqn:max-attack}
\end{equation}

Note Eqn. \ref{eqn:max-attack} differs from Eqn. \ref{eqn:avg-attack} by adding an $\Omega$ constraint to balance between ``commodity" and ``specificity". The found $x_{max}^{adv}$ \textbf{both} strongly increases the averaged loss values from all $f_{i}$s (therefore possessing transferablity), \textbf{and} maximally fools one individual $f_{\theta_{i}}$s as it is selected from the collection $\Omega$ of single attacks.

\vspace{-0.5em}
\subsection{Defence on Multi-Output Networks}\label{sec:3-3}
\vspace{-0.5em}

For simplicity and fair comparison, we focus on adversarial training \citep{madry2017towards} as our defense framework, where the three above defined attack forms can be plugged-in to generate adversarial images to augment training, as follows ($\Theta$ is the union of learnable parameters):
\begin{equation}
\hspace{-0.01in} \theta_{i}\in\Theta \mbox{, where } \theta_{i} = \argmin_{\theta'}|\phi(f_{\theta_{i}}(x), y) + \phi(f_{\theta_{i}}(x^{adv}), y) |.\label{eqn:sur-2},
\end{equation}
where $x^{adv} \in \{x_{i}^{adv}, x_{avg}^{adv}, x_{max}^{adv}\}$. As $f_{i}$s partially share their weights $\theta_{i}$ in a multi-output network, the updates from different $f_{i}$s will be averaged on the shared parameters.

\vspace{-0.5em}

\section{Experimental Results}\label{sec:performance}
\vspace{-0.5em}
\subsection{Experimental Setup}
\vspace{-0.5em}
\paragraph{Evaluation Metrics}
We evaluate accuracy, robustness, and efficiency, using the metrics below:
\vspace{-0.1in}
\begin{itemize}[leftmargin=15pt, itemsep=1pt]
\item Testing Accuracy (\textbf{TA}): the classification accuracy on the original clean test set.
\item Adversarial Testing Accuracy (\textbf{ATA}): Given an attacker, ATA stands for the classification accuracy on the attacked test set. It is the same as the ``robust accuracy" in \citep{zhang2019theoretically}.
\item Mega Flops (\textbf{MFlops}): The number of million floating-point multiplication operations consumed on the inference, averaged over the entire testing set.
\end{itemize}
\vspace{-0.1in}

\vspace{-0.5em}

\paragraph{Datasets and Benchmark Models} We evaluate three representative CNN models on two popular datasets: SmallCNN on MNIST \citep{qi2018nipe}; ResNet-38 \citep{he2016resnet} and MobileNet-V2 \citep{sandler2018mobilenetv2} on CIFAR-10. The three networks span from simplest to more complicated, and covers a compact backbone. All three models are defended by adversarial training, constituting \textbf{strong baselines}. Table \ref{tab:benchmark} reports the models, datasets, the attacker algorithm used in attack \& defense, and thee TA/ATA/MFlops performance of three defended models.

\vspace{-0.5em}
\paragraph{Attack and Defense on RDI-Nets}
We build RDI-Nets by appending side branch outputs for each backbone. For SmallCNN, we add two side branches ($K = 2$). For ResNet-38 and MobileNet-V2, we have $K = 6$ and $K = 2$, respectively. The branches are designed to cause negligible overheads: more details of their structure and positions can be referenced in Appendix \ref{sec:6-2}. We call those result models RDI-SmallCNN, RDI-ResNet38 and RDI-MobileNetV2 hereinafter.

    We then generate attacks using our three defined forms. Each attack form could be solved with various attacker algorithms (e.g.,\  PGD, C\&W, FGSM), and by default we solve it with the same attacker used for each backbone in Table \ref{tab:benchmark}. If we fix one attacker algorithm (e.g.,\  PGD), then TA/ATA for a single-output network can be measured without ambiguity. Yet for ($K$+1)-output RDI-Nets, there could be at least $K$+3 different ATA numbers for one defended model, depending on what attack form in Section \ref{sec:3-1} to apply ($K$+1 single attacks, 1 average attack, and 1 max-average attack). For example, we denote by \textbf{ATA (Branch1)} the ATA number when applying the single attack generated from the first side output branch (e.g.,\  $x_{1}^{adv}$); similarly elsewhere.

We also defend RDI-Nets using adversarial training, using the forms of adversarial images to augment training. By default, we adopt three adversarial training defense schemes: \textbf{Main Branch} (single attack using $x_{K+1}^{adv}$)\footnote{We tried adversarial training using other $K$ earlier side branch single attacks, and found their TA/ATA to be much more deteriorated compared to the main branch one. We thus report this only for compactness.},\textbf{Average} (using $x_{avg}^{adv}$),and \textbf{Max-Average} (using $x_{avg}^{max}$), in addition to the undefended RDI-Nets (using standard training) denoted as \textbf{Standard}.

We cross evaluate ATAs of different defenses and attacks, since an ideal defense shall protect against all possible attack forms. To faithfully indicate the actual robustness, we choose the \textit{lowest number} among all $K+3$ ATAs, denoted as \textbf{ATA (Worst-Case)}, as the robustness measure for an RDI-Net.

\begin{table}[ht]
\centering
\vspace{-1.5em}
\caption{Benchmarking results of adverserial training of three networks. PGD-40 denotes running the projected gradient descent attacker \citep{madry2017towards} for 40 iterations. We set the perturbation size as $0.3$ for MNIST and $8/255$ for CIFAR-$10$ in $\ell_\infty$ norm (adopted by all following experiments).}
\vspace{0.1in}
\label{tab:benchmark}
\resizebox{0.7\textwidth}{!}{
\begin{tabular}{lllllll}
\hline
\textbf{Model} & \textbf{Dataset} & \textbf{Defend} & \textbf{Attack} & \textbf{TA} & \textbf{ATA} & \textbf{MFlops} \\ \hline
SmallCNN       & MNIST            & PGD-40          & PGD-40          & 99.49\%     & 96.31\%      & 9.25           \\
ResNet-38      & CIFAR-10         & PGD-10          & PGD-20          & 83.62\%     & 42.29\%      & 79.42          \\
MobileNetV2   & CIFAR-10         & PGD-10          & PGD-20          & 84.42\%     & 46.92\%      & 86.91    \\ \hline
\end{tabular}}
\vspace{-1em}
\end{table}

\vspace{-0.5em}
\subsection{Evaluation and Analysis}
\vspace{-0.5em}
\paragraph{MNIST Experiments}

The MNIST experimental results on RDI-SmallCNN are summarized in table \ref{tab:smallcnn-branch}, with several meaningful observations to be drawn. \underline{First}, the undefended models (Standard) are easily compromised by all attack forms. \underline{Second}, The single attack-defended model (Main Branch) achieves the best ATA against the same type of attack, i.e.,\  ATA (Main Branch), and also seems to boost the closest output branch's robustness, i.e.,\  ATA (Branch 2). However, its defense effect on the further-away Branch 1 is degraded, and also shows to be fragile under two stronger attacks (Average, and Max-Average). \underline{Third}, both Average and Max-Average defenses achieve good TAs, as well as ATAs against all attack forms (and therefore Worst-Case), with Max-Average slightly better at both (the margins are small due to the data/task simplicity; see next two).

Moreover, compared to the strong baseline of SmallCNN defended by PGD (40 iterations)-based adversarial training, RDI-SmallCNN with Max-Average defense wins in terms of both TA and ATA. Impressively, that comes together with 34.30\% computational savings compared to the baseline. Here the different defense forms do not appear to alter the inference efficiency much: they all save around 34\% - 36\% MFlops compared to the backbone.

\vspace{-0.15in}
\begin{table}[ht]
\centering
\caption{The performance of RDI-SmallCNN. The "\textbf{Average MFlops}" is calculated by averaging the total flop costs consumed over the inference of the entire set (different samples take different FLOPs due to input-adaptive inference). The perturbation size and step size are $0.3$ and $0.01$, respectively.}
\vspace{0.1in}
\label{tab:smallcnn-branch}
\resizebox{0.8\textwidth}{!}{
\begin{tabular}{|l|l|l|l|l|}
\hline
Defense Method & Standard  & Main Branch & Average  & Max-Average  \\ \hline
\hline
\textbf{TA} & 99.48\% & 99.50\% & 99.51\% & \textbf{99.52\%} \\ \hline
\hline
ATA (Branch 1) & 6.60\% & 60.50\% & 98.69\% & 98.52\% \\ \hline
ATA (Branch 2) & 3.16\% & 98.14\% & 97.64\% & 97.62\% \\ \hline
ATA (Main Branch) & 1.32\% & 96.70\% & 96.30\% & 96.43\% \\ \hline
ATA (Average) & 2.61\% & 61.35\% & 97.37\% & 97.42\% \\ \hline
ATA (Max-Average) & 2.10\% & 61.83\% & 96.82\% & 96.89\% \\ \hline
\textbf{ATA (Worst-Case)} & 1.32\% & 60.50\% & 96.30\% & \textbf{96.43\%} \\ \hline
\hline
Average MFlops & 5.89 & 5.89 & 5.95 & 6.08 \\ \hline
\textbf{Computation Saving} & 36.40\% & 36.40\% & 35.70\% & 34.30\% \\ \hline
\end{tabular}
}
\vspace{-0.5em}
\end{table}

\vspace{-0.5em}
\paragraph{CIFAR-10 Experiments}\label{sec:4-2}
The results on RDI-ResNet38 and RDI-MobileNetV2 are presented in Tables \ref{tab:bnet38-branch} and \ref{tab:mnv2-branch}, respectively. Most findings seem to concur with MNIST experiments. \underline{Specifically}, on the more complicated CIFAR-10 classification task, Max-Average defense achieves much more obvious margins over Average defense, in terms of ATA (Worst-Case): 2.79\% for RDI-ResNet38, and 1.06\% for RDI-MobileNetV2. Interestingly, the Average defense is not even the strongest in defending average attacks, as Max-Average defense can achieve higher ATA (Average) in both cases. We conjecture that averaging all branch losses might ``over-smooth'' and diminish useful gradients.

Compared to the defended ResNet-38 and MobileNet-V2 backbones, RDI-Nets with Max-Average defense achieve higher TAs and ATAs for both. Especially, the ATA (Worst-Case) of RDI-ResNet-38 surpasses the ATA of ResNet-38 defended by PGD-adversarial training by \textbf{1.03\%}, while saving around \textbf{30\%} inference budget. We find that different defenses on CIFAR-10 have more notable impacts on computational saving. Seemingly, a stronger defense (Max-Average) requires inputs to go through the scrutiny of more layers on average, before outputting  confident enough predictions: a sensible observation as we expect.

\vspace{-0.5em}
\paragraph{Visualization of Adaptive Inference Behaviors} We visualize the exiting behaviors of RDI-ResNet38 in Fig \ref{fig:exit_behavior}. We plot each branch exiting percentage on clean set and adversarial sets (worst-case) of examples. A few interesting observations can be found. \underline{First}, we observe that the single-attack defended model can be easily fooled as adversarial examples can be routed through other less-defended outputs (due to the limited transferability of attacks between different outputs). \underline{Second}, the two stronger defenses (Average and Max-Average) show much more uniform usage of multiple outputs. Their routing behaviors for clean examples are almost identical. For adversarial examples, Max-Average tends to call upon the full inference more often (i.e., more ``conservative'').

\begin{table}[ht]
\centering
\vspace{-1em}
\caption{The performance evaluation on RDI-ResNet38. The perturbation size and step size are $8/255$ and $2/255$, respectively.}
\vspace{0.1in}
\label{tab:bnet38-branch}
\resizebox{0.8\textwidth}{!}{
\begin{tabular}{|l|l|l|l|l|}
\hline
Defence Method & Standard & Main Branch & Average & Max-Average \\ \hline
\hline
\textbf{TA} & 92.43\% & 83.74\% & 82.42\% & \textbf{83.79\%} \\ \hline
\hline
ATA (Branch1) & 0.12\% & 12.02\% & 71.56\% & 69.71\% \\ \hline
ATA (Branch2) & 0.01\% & 5.58\% & 66.67\% & 63.11\% \\ \hline
ATA (Branch3) & 0.04\% & 42.73\% & 60.65\% & 60.72\% \\ \hline
ATA (Branch4) & 0.06\% & 34.95\% & 50.17\% & 47.82\% \\ \hline
ATA (Branch5) & 0.06\% & 41.77\% & 44.83\% & 45.53\% \\ \hline
ATA (Branch6) & 0.11\% & 41.68\% & 45.83\% & 44.12\% \\ \hline
ATA (Main Branch) & 0.13\% & 42.74\% & 47.52\% & 49.82\% \\ \hline
ATA (Average) & 0.01\% & 9.14\% & 42.09\% & 43.32\% \\ \hline
ATA (Max-Average) & 0.01\% & 7.15\% & 40.53\% & 43.43\% \\ \hline
\textbf{ATA (Worst-Case)} & 0.01\% & 5.58\% & 40.53\% & \textbf{43.32\%} \\ \hline
\hline
Average MFlops & 29.41 & 48.27 & 56.90 & 57.81 \\ \hline
\textbf{Computation Saving} & 62.96\% & 39.20\% & 28.35\% & 27.20\% \\ \hline
\end{tabular}}
\vspace{-1em}
\end{table}

\begin{table}[ht]
\centering
\caption{The performance evaluation on RDI-MobilenetV2. The perturbation size and step size are $8/255$ and $2/255$, respectively.}
\vspace{0.1in}
\label{tab:mnv2-branch}
\resizebox{0.8\textwidth}{!}{
\begin{tabular}{|l|l|l|l|l|}
\hline
Defence Method & Standard & Main Branch & Average & Max-Average \\ \hline
\hline
\textbf{TA} & 93.22\% & \textbf{85.28\%} & 82.14\% & 84.91\% \\ \hline
\hline
ATA (Branch1) & 0.35\% & 37.40\% & 67.65\% & 71.78\% \\ \hline
ATA (Branch2) & 0\% & 47.35\% & 50.38\% & 50.15\% \\ \hline
ATA (Main Branch) & 0\% & 46.69\% & 49.33\% & 46.99\% \\ \hline
ATA (Average) & 0\% & 35.20\% & 45.93\% & 47.00\% \\ \hline
ATA (Max-Average) & 0\% & 36.66\% & 49.33\% & 50.18\% \\ \hline
\textbf{ATA (Worst-Case)} & 0\% & 35.20\% & 45.93\% & \textbf{46.99\%} \\ \hline
\hline
Average MFlops & 49.78 & 52.81 & 58.23 & 60.84 \\ \hline
\textbf{Computation Saving} & 42.72\% & 39.23\% & 33.00\% & 29.99\% \\ \hline
\end{tabular}}
\vspace{-0.5em}
\end{table}

\begin{figure}[ht]
\begin{tabular}{ccc}
\subfigure[]{
\centering
   \includegraphics[height=0.14\textheight]{figures/exit_main.pdf}
   \label{fig:exit_main}
}
\quad
\centering
\subfigure[]{
   \includegraphics[height=0.14\textheight]{figures/exit_avg.pdf}
   \label{fig:exit_avg}
}
\quad
\subfigure[]{
\centering
   \includegraphics[height=0.14\textheight]{figures/exit_max.pdf}
   \label{fig:exit_max_avg}
}
\end{tabular}
\label{fig:exit_behavior}
\vspace{-0.2in}
\caption{ The exiting behaviours of RDI-ResNet38 defended by (a) Single attack defense (Main Branch); (b) Average defense; and (c) Max-Average defense.}
\vspace{-1em}
\end{figure}

\vspace{-0.5em}
\subsection{Comparison with Defended Sparse Networks}
\vspace{-0.5em}

An alternative to achieve accuracy-robust-efficiency trade-off is by defending a sparse or compressed model.
Inspired by \citep{guo2018sparse,gui2019model}, we compare RDI-Net with Max-Average defense to the following baseline: first compressing the network with a state-of-the-art model compression method~\citep{huang2018spase}, and then defend the compressed network using the PGD-10 adversarial training.
We sample different sparsity ratios in~\citep{huang2018spase} to obtain models of different complexities.
Fig. \ref{fig:comparision} in Appendix visualizes the comparison on ResNet-38: for either method, we sample a few models of different MFLOPs. At similar inference costs (e.g., 49.38M for pruning + defense, and 48.35M for RDI-Nets), our proposed approach consistently achieves higher ATAs ($>$ 2\%) than the strong pruning + defense baseline, with higher TAs.

\begin{wraptable}{r}{7cm}
\vspace{-2em}
\caption{Performance comparison between RDI-ResNet38 and ATMC.}\label{tab:atmc}
\vspace{-0.5em}
\begin{center}
\resizebox{0.5\textwidth}{!}{
\begin{tabular}{lccc}
\multicolumn{1}{c}{\bf Methods} & \multicolumn{1}{c}{\bf TA} & \multicolumn{1}{c}{\bf ATA} & \multicolumn{1}{c}{\bf MFlops}
\\ \hline
ATMC (\cite{gui2019model}) & \textbf{83.81} & 43.02 & 56.82 \\
RDI-ResNet-38 (\textbf{Worst-Case}) & 83.79 & \textbf{43.32} & 57.81 \\
\end{tabular}}
\end{center}
\vspace{-1em}
\end{wraptable}

We also compare with the latest ATMC algorithm \citep{gui2019model} that jointly optimizes robustness and efficiency, applied the same ResNet-38 backbone. As shown in Table~\ref{tab:atmc}, at comparable MFlops, RDI-ResNet-38 surpasses ATMC by 0.3\% in terms of ATA, with a similar TA.

\vspace{-0.5em}
\subsection{Generalized Robustness Against Other Attackers}
\vspace{-0.5em}
In the aforementioned experiments, we have only evaluated on RDI-Nets against ``deterministic" PGD-based adversarial images. We show that RDI-Nets also achieve better generalized robustness against other ``randomized" or unseen attackers. We create the new ``random attack": that attack will randomly combine the multi-exit losses, and summarize the results in Table \ref{tab:random}. We also follow the similar setting in \cite{gui2019model} and report the results against FGSM \citep{goodfellow2014explaining} and WRM \citep{simha2018wrm} attacker, in Tables \ref{tab:fgsm} and \ref{tab:wrm} respectively (more complete results can be found in Appendix \ref{sec:6-4}).

\begin{table}[ht]
\centering
\vspace{-0.5em}
\caption{Performance on RDI-ResNet38 against random attack. The perturbation size and step size are $8/255$ and $2/255$, respectively. More details of random attack can be referenced in Appendix \ref{sec:6-4}.}
\vspace{0.1in}
\label{tab:random}
\resizebox{0.8\textwidth}{!}{
\begin{tabular}{|l|l|l|l|l|}
\hline
Defence Method & Standard & Main Branch & Average & Max-Average \\ \hline
\hline
\textbf{TA} & 92.43\% & 83.74\% & 82.42\% & \textbf{83.79\%} \\ \hline
\hline
ATA (Random) & 0.01\% & 10.33\% & 43.11\% & \textbf{44.86}\% \\ \hline
\hline
Average MFlops & 27.33 & 52.36 & 55.21 & 56.54 \\ \hline
\textbf{Computation Saving} & 65.58\% & 34.07\% & 30.48\% & 28.80\% \\ \hline
\end{tabular}}
\vspace{-0.8em}
\end{table}

\begin{table}[ht]
\centering
\vspace{-0.5em}
\caption{Performance on RDI-ResNet38 (defended with PGD) against FGSM attack (perturbation size is $8/255$). The original defended ResNet38 by PGD under the same attack has ATA $51.11\%$.}
\vspace{0.1in}
\label{tab:fgsm}
\resizebox{0.8\textwidth}{!}{
\begin{tabular}{|l|l|l|l|l|}
\hline
Defence Method & Standard & Main Branch & Average & Max-Average \\ \hline
\hline
\textbf{TA} & 92.43\% & 83.74\% & 82.42\% & \textbf{83.79\%} \\ \hline
\hline
ATA (Main Branch) & 11.51\% & 51.45\% & 53.64\% & 54.72\% \\ \hline
ATA (Average) & 11.41\% & 50.21\% & 51.81\% & 53.20\% \\ \hline
ATA (Max-Average) & 2.09\% & 47.53\% & 50.63\% & 52.40\% \\ \hline
\textbf{ATA (Worst-Case)} & 2.09\% & 47.53\% & 50.63\% & \textbf{51.05\%} \\ \hline
\hline
Average MFlops & 65.74 & 55.27 & 58.27 & 59.67 \\ \hline
\textbf{Computation Saving} & 17.21\% & 30.40\% & 26.40\% & 24.86\% \\ \hline
\end{tabular}}
\vspace{-0.8em}
\end{table}

\begin{table}[!ht]
\centering
\vspace{-0.5em}
\caption{Performance on RDI-ResNet38 (defended with PGD) against WRM attack (perturbation size is $0.3$). The original defended ResNet38 by PGD under the same attack has ATA $83.35\%$.}
\vspace{0.1in}
\label{tab:wrm}
\resizebox{0.8\textwidth}{!}{
\begin{tabular}{|l|l|l|l|l|}
\hline
Defence Method & Standard & Main Branch & Average & Max-Average \\ \hline
\hline
\textbf{TA} & 92.43\% & 83.74\% & 82.42\% & \textbf{83.79\%} \\ \hline
\hline
ATA (Main Branch) & 34.42\% & 83.74\% & 82.42\% & 83.78\% \\ \hline
ATA (Average) & 26.48\% & 83.69\% & 82.36\% & 83.77\% \\ \hline
ATA (Max-Average) & 23.51\% & 83.73\% & 82.40\% & 83.78\% \\ \hline
\textbf{ATA (Worst-Case)} & 23.51\% & 83.69\% & 82.36\% & \textbf{83.77\%} \\ \hline
\hline
Average MFlops & 50.05 & 50.46 & 52.89 & 52.38 \\ \hline
\textbf{Computation Saving} & 36.98\% & 36.46\% & 33.40\% & 34.04\% \\ \hline
\end{tabular}}
\vspace{-0.8em}
\end{table}

\section{Discussion and Analysis}
\vspace{-0.5em}
\paragraph{Intuition: Multi-Output Networks as Special Ensembles} Our intuition on defending multi-output networks arises from the success of ensemble defense in improving both accuracy and robustness \citep{tramer2017ensemble,strauss2017ensemble}, which also aligns with the model capacity hypothesis \citep{nakkiran2019adversarial}. A general multi-output network \citep{xu2019survey} could be decomposed by an ensemble of single-output models, with weight re-using enforced among them. It is thus more compact than an ensemble of independent models, and the extent of sharing weight calibrates ensemble diversity versus efficiency. Therefore, we expect a defended multi-output network to (mostly) inherit the strong accuracy/robustness of ensemble defense, while keeping the inference cost lower.
\vspace{-0.5em}
\paragraph{Do "Triple Wins" Go Against the Model Capacity Needs?} We point out that our seemingly ``free'' efficiency gains (e.g.,\  not sacrificing TA/ATA) do not go against the current belief that a more accurate and robust classifier relies on a larger model capacity \citep{nakkiran2019adversarial}. From the visualization, there remains to be a portion of clean/adversarial examples that have to utilize the full inference to predict well. In other words, the full model capacity is still necessary to achieve our current TAs/ATAs. Meanwhile, just like in standard classification \citep{wang2017skipnet}, not all adversarial examples are born equally. Many of them can be predicted using fewer inference costs (taking earlier exits). Therefore, RDI-Nets reduces the ``effective model capacity'' averaged on all testing samples for overall higher inference efficiency, while not altering the full model capacity.

\section{Conclusion}
\vspace{-0.5em}
This paper targets to simultaneously achieve high accuracy and robustness and meanwhile keeping inference costs lower. We introduce the multi-output network and input-adaptive dynamic inference, as a strong tool to the adversarial defense field for the first time. Our RDI-Nets achieve the ``triple wins'' of better accuracy, stronger robustness, and around 30\% inference computational savings. Our future work will extend RDI-Nets to more dynamic inference mechanisms.

\section{Acknowledgement}
\vspace{-0.5em}
We would like to thank Dr. Yang Yang from Walmart Technology for highly helpful discussions throughout this project.
\section{Learning Details of RDI-Nets}\label{sec:6-1}
\paragraph{MNIST}
We adopt the network architecture from \citep{qi2018nipe} with four convolutions and three full-connected layers. We train for $13100$ iterations with a batch size of $256$. The learning rate is initialized as $0.033$ and is lowered by $10$ at $12000$th and $12900$th iteration. For hybrid loss, the weights $\{w_{i}\}_{i=1}^{N+1}$ are set as $\{1, 1, 1\}$ for simplicity. For adversarial defense/attack, we perform 40-steps PGD for both defense and evaluation. The perturbation size and step size are set as $0.3$ and $0.01$.
\paragraph{CIFAR-10}
We take ResNet-38 and MobileNetV2 as the backbone architectures. For RDI-ResNet38, we initialize learning rate as $0.1$ and decay it by a factor of 10 at $32000$th and $48000$th iteration. The learning procedure stops at $55000$ iteration. For RDI-MobileNetV2, the learning rate is set to $0.05$ and is lowered by $10$ times at $62000$th and $70000$th iteration. We stop the learning procedure at $76000$ iteration. For hybrid loss, we follow the discussion in \citep{hu2019anytime} and set $\{w_{i}\}_{i=1}^{N+1}$ of RDI-ResNet38 and RDI-MobileNetV2 as $\{0.5, 0.5, 0.7, 0.7, 0.9, 0.9, 2\}$ and $\{0.5, 0.5, 1\}$, respectively. For adversarial defense/attack, the perturbation size and step size are set as $8/255$ and $2/255$. 10-steps PGD is performed for defense and 20-steps PGD is utilized for evaluation.

\section{Network Structure of RDI-Nets}\label{sec:6-2}
To build RDI-Nets, we follow the similar setting in \cite{teerapittayanon2016branchynet} by appending additional branch classifiers at equidistant points throughout a given network, as illustrated in Fig \ref{fig:smallCNN}, Fig \ref{fig:res38} and Fig \ref{fig:mob}. A few pooling operations, light-weight convolutions and fully-connected layers are appended to each branch classifiers. Note that the extra flops introduced by side branch classifiers are less than 2\% than the original ResNet-38 or MobileNetV2.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.7\linewidth]{figures/smallCNN.pdf}

    \caption{Network architecture of RDI-SmallCNN. Two branch classifiers are inserted after $1st$ convolutional layer and $3$rd convolutional layer in the original SmallCNN.}
    \label{fig:smallCNN}

\end{figure}

\begin{figure}[ht]
    \centering
    \includegraphics[width=1\linewidth]{figures/resnet38.pdf}
    \vspace{-1em}
    \caption{Network architecture of RDI-ResNet38. In each residual block group, two branch classifiers are inserted after $1st$ residual block and $4$th residual block.}
    \label{fig:res38}
    \vspace{-1em}
\end{figure}
\begin{figure}[ht]
    \centering
    \includegraphics[width=1\linewidth]{figures/mob.pdf}
    \vspace{-1em}
    \caption{Network architecture of RDI-MobilenetV2. Two branch classfiers are inserted after $3rd$ inverted residual block and $11$th inverted residual block in the orignal MobilenetV2.}
    \label{fig:mob}
    \vspace{-1em}
\end{figure}
\section{Input-Adaptive Inference for RDI-Nets}\label{sec:6-3}
Similar to the deterministic strategy in \cite{teerapittayanon2016branchynet}, we adopt the entropy as the measure of the prediction confidence. Given a prediction vector $y\in{\mathbb{R}}^{C}$, where $C$ is the number of class, the entropy of $y$ is defined as follow,
\begin{equation}
-\sum_{c=1}^{C}(y_{c} + \epsilon)log(y_{c} + \epsilon),
\end{equation}
where $\epsilon$ is a small positive constant used for robust entropy computation. To perform fast inference on a ($K$+1)-output RDI-Net, we need to determine $K$ threshold numbers, i.e.,\  $\{t_{i}\}_{i=1}^{K}$, so that the input $x$ will exit at $i$th branch if the entropy of $y_{i}$ is larger than $t_{i}$. To choose $\{t_{i}\}_{i=1}^{K}$, \cite{huang2018multiscale} provides a good starting point  by fixing exiting probability of each branch classifiers equally on validation set so that each sample can equally contribute to inference. We follow this strategy but adjust the thresholds to make the contribution of middle branches slightly larger than the early branches. The threshold numbers for RDI-SmallCNN, RDI-ResNet38, and RDI-MobilenetV2 are set to be $\{0.023, 0.014\}$, $\{0.32, 0.36, 0.39, 0.83, 1.12, 1.35\}$, and $\{0.267, 0.765\}$, respectively.

\begin{figure}[ht]
    \centering

    \includegraphics[width=0.8\linewidth]{figures/comparision.pdf}

    \caption{Performance comparison between RDI-Net and the pruning + defense baseline. Each marker represents a model, whose size is proportional to its MFlops. $\gamma$ is the sparsity trade-off parameter: the larger the sparser (smaller model).
}
    \label{fig:comparision}
\end{figure}

\section{Generalized Robustness}\label{sec:6-4}
Here, we introduce the attack form of random attack and report the complete results against FGSM \citep{goodfellow2014explaining} and WRM \citep{simha2018wrm} attacker under various attack forms, in Tables \ref{tab:complete_fgsm} and \ref{tab:complete_wrm}, respectively.

\paragraph{Random Attack} the attack exploits multi-loss flexibility by randomly fusing all $f_{\theta_{i}}$ losses. Given a $N$-output network, we have a fusion vector $C \in {\mathbb{R}}^{N} \sim {\mathbb{D}}$, where ${\mathbb{D}}$ is some distribution (uniform by default). We denote $c_{j}$ as the $j$th element of $C$ and $x^{adv}_{rdm}$ can be found by:
\begin{equation}
\hspace{-0.01in} x_{rdm}^{adv} = \argmax_{x'\in |x'-x|_{\infty}\leq\epsilon}|\frac{1}{N}\sum_{j=1}^{N}c_{j}\phi(f_{\theta_{j}}(x'), y)|.
\label{eqn:rdm-attack}
\end{equation}
It is expected to challenge our defense, due to the infinitely many ways of randomly fusing outputs.

\begin{table}[ht]
\centering
\vspace{-0.5em}
\caption{The performance evaluation on RDI-ResNet38 (defended with PGD) against FGSM attack. The perturbation size is $8/255$. The ATA of the original defended ResNet38 by PGD under the same attacker is $51.11\%$.}
\vspace{0.1in}
\label{tab:complete_fgsm}
\begin{tabular}{|l|l|l|l|l|}
\hline
Defence Method & Standard & Main Branch & Average & Max-Average \\ \hline
\hline
\textbf{TA} & 92.43\% & 83.74\% & 82.42\% & \textbf{83.79\%} \\ \hline
\hline
ATA (Branch1) & 20.69\% & 66.06\% & 72.77\% & 72.76\% \\ \hline
ATA (Branch2) & 16.15\% & 53.87\% & 70.40\% & 69.71\% \\ \hline
ATA (Branch3) & 8.13\% & 63.70\% & 64.19\% & 65.14\% \\ \hline
ATA (Branch4) & 10.09\% & 56.67\% & 58.45\% & 58.20\% \\ \hline
ATA (Branch5) & 9.45\% & 50.81\% & 52.76\% & 52.96\% \\ \hline
ATA (Branch6) & 10.22\% & 50.34\% & 53.17\% & 51.05\% \\ \hline
ATA (Main Branch) & 11.51\% & 51.45\% & 53.64\% & 54.72\% \\ \hline
ATA (Average) & 11.41\% & 50.21\% & 51.81\% & 53.20\% \\ \hline
ATA (Max-Average) & 2.09\% & 47.53\% & 50.63\% & 52.40\% \\ \hline
\textbf{ATA (Worst-Case)} & 2.09\% & 47.53\% & 50.63\% & \textbf{51.05\%} \\ \hline
\hline
Average MFlops & 65.74 & 55.27 & 58.27 & 59.67 \\ \hline
\textbf{Computation Saving} & 17.21\% & 30.40\% & 26.40\% & 24.86\% \\ \hline
\end{tabular}
\vspace{-0.5em}
\end{table}

\begin{table}[!ht]
\centering
\vspace{-0.5em}
\caption{The performance evaluation on RDI-ResNet38 (defended with PGD) against WRM attack. The perturbation size is $0.3$. The ATA of the original defended ResNet38 by PGD under the same attacker is $83.35\%$.}
\vspace{0.1in}
\label{tab:complete_wrm}
\begin{tabular}{|l|l|l|l|l|}
\hline
Defence Method & Standard & Main Branch & Average & Max-Average \\ \hline
\hline
\textbf{TA} & 92.43\% & 83.74\% & 82.42\% & \textbf{83.79\%} \\ \hline
\hline
ATA (Branch1) & 46.60\% & 83.73\% & 82.42\% & 83.78\% \\ \hline
ATA (Branch2) & 71.33\% & 83.73\% & 82.42\% & 83.79\% \\ \hline
ATA (Branch3) & 23.51\% & 83.73\% & 82.41\% & 83.78\% \\ \hline
ATA (Branch4) & 33.41\% & 83.73\% & 82.42\% & 83.78\% \\ \hline
ATA (Branch5) & 42.35\% & 83.73\% & 82.41\% & 83.78\% \\ \hline
ATA (Branch6) & 47.77\% & 83.74\% & 82.40\% & 83.78\% \\ \hline
ATA (Main Branch) & 34.42\% & 83.74\% & 82.42\% & 83.78\% \\ \hline
ATA (Average) & 26.48\% & 83.69\% & 82.36\% & 83.77\% \\ \hline
ATA (Max-Average) & 23.51\% & 83.73\% & 82.40\% & 83.78\% \\ \hline
\textbf{ATA (Worst-Case)} & 23.51\% & 83.69\% & 82.36\% & \textbf{83.77\%} \\ \hline
\hline
Average MFlops & 50.05 & 50.46 & 52.89 & 52.38 \\ \hline
\textbf{Computation Saving} & 36.98\% & 36.46\% & 33.40\% & 34.04\% \\ \hline
\end{tabular}
\vspace{-0.5em}
\end{table}