[table]capposition=top

# Risk-averse controller design against data injection attacks on actuators
for uncertain control systems

Sribalaji C. Anand1 and André M. H. Teixeira2 *This work is supported by the
Swedish Research Council under the grant 2018-04396 and by the Swedish
Foundation for Strategic Research.1 Sribalaji C. Anand is with the Department
of Electrical Engineering, Uppsala University, PO Box 65, SE-75103, Uppsala,
Sweden<EMAIL_ADDRESS>André M. H. Teixeira is with the
Department of Information Technology, Uppsala University, PO Box 337, SE
-75105, Uppsala, Sweden<EMAIL_ADDRESS>
###### Abstract

In this paper, we consider the optimal controller design problem against data
injection attacks on actuators for an uncertain control system. We consider
attacks that aim at maximizing the attack impact while remaining stealthy in
the finite horizon. To this end, we use the Conditional Value-at-Risk to
characterize the risk associated with the impact of attacks. The worst-case
attack impact is characterized using the recently proposed output-to-output
$\ell_{2}$-gain (OOG). We formulate the design problem and observe that it is
non-convex and hard to solve. Using the framework of scenario-based
optimization and a convex proxy for the OOG, we propose a convex optimization
problem that approximately solves the design problem with probabilistic
certificates. Finally, we illustrate the results through a numerical example.

##  I INTRODUCTION

Cyber-physical systems (CPSs) represent a large class of networked control
systems where the physical world and the digital infrastructure are tightly
coupled, such as smart cities, autonomous systems, transportation networks,
and Internet Of Things. However, the trend towards increased usage of open-
standard communication protocols among control systems has made these systems
vulnerable to online cyber-attacks such as Stuxnet [1], Industroyer [2], etc.
Such cyber-attacks can negatively affect the operation of CPS [3].

Significant work is done in detecting and mitigating cyber-attacks (see [4, 5]
and references therein). For instance, [6] designs an optimal controller in
the presence of covert attacks. The limitation of [6] is, it approximates the
risk metric Conditional Value-at-Risk (CVaR) empirically using samples, and it
parameterizes the controller as a finite family of Finite Impulse Response
(FIR) filters. Although these approximations simplify the problem, the
validity of these approximations is not discussed except in the asymptotic
case i.e., as the number of samples for empirical approximation and the number
of FIR filters tends to infinity.

The article [7] proposes and solves two controller design problems. Firstly,
it proposes a convex design problem such that the volume of the reachable set
of states by the adversary is minimized. Secondly, it proposes a convex design
problem that maximizes the Euclidean distance between the set of states
reachable by the adversary and the set of critical states. A similar approach
was also adopted in [8]. However, both of these works do not consider an
uncertain system. The works [9] and [10] addresses the issue of jointly
designing the controller and detector against false data injection (FDI)
attacks. However, they also assume a deterministic system.

This paper addresses some of the existing limitations in the literature, by
investigating the optimal controller design problem against FDI attacks on
actuators for an uncertain control system. To this end, we adopt the following
setup. We consider a discrete-time (DT) linear time-invariant (LTI) process
with parametric uncertainty, a static output feedback controller, and an
anomaly detector. An adversary with perfect system knowledge injects false
data into the actuators. In reality, it is hard for the adversary to have
perfect system knowledge, but this assumption helps to study the worst-case.
The system operator (or the defender) knows only about the bounds of the
uncertainty. Under this setup, we present the following contributions.

  1. 1.

Firstly, we formulate the risk-averse design problem. Here, for a given
realization of the uncertainty, we use the output-to-output $\ell_{2}$-gain
(OOG) [11] to characterize the worst-case impact. We then use the CVaR to
characterize the risk associated with the attack impact. The advantages of
using the OOG over the classical $H_{\infty}/H_{\\_}$ metrics were
demonstrated in [9]. We also observe that the design problem corresponds to an
untractable infinite non-convex optimization problem.

  2. 2.

Secondly, extending the results of [12], we derive an upper bound for the OOG.
Using this upper bound, we relax the infinite non-convex design problem into
an infinite convex design problem.

  3. 3.

Finally, by adopting the scenario-based approach [13], we modify the infinite
convex optimization problem into its sampled counterpart. We also provide
probabilistic guarantees on the infinite design problem based on the number of
samples used to formulate the sampled optimization problem and the dimension
of the controller. The advantage of using scenario-based approach over other
approaches is discussed in [14].

To the best of the author’s knowledge, the problem of risk-sensitive
controller design for an uncertain control in the finite horizon against FDI
attacks has not been addressed in the literature.

The remainder of this paper is organized as follows. Section II describes the
problem background. The design problem is formulated in Section III. The
problem is relaxed and convexified in Section IV. Section V approximates the
problem empirically using the scenario-based approach. We illustrate the
results using a numerical example in Section VI. Finally, we provide
concluding remarks in Section VII.

##  II Problem background

In this section, we describe the control system structure and the goal of the
adversary. Consider the general description of a finite horizon closed-loop DT
LTI system with a process ($\mathcal{P}$) with parametric uncertainty, a
static output feedback controller ($\mathcal{C}$) and an anomaly detector
($\mathcal{D}$) as shown in Fig. 1. The closed-loop system is represented by

$\displaystyle\mathcal{P}:$
$\displaystyle\left\\{\begin{array}[]{ll}x_{p}[k+1]&=A^{\Delta}x_{p}[k]+B\tilde{u}[k]\\\
y[k]&=Cx_{p}[k]\\\ y_{p}[k]&=C_{J}x_{p}[k]\end{array}\right.$ (4)
$\displaystyle\mathcal{C}:$
$\displaystyle\left\\{\begin{array}[]{ll}{u}[k]&=Ky[k]\end{array}\right.$ (6)
$\displaystyle\mathcal{D}:$
$\displaystyle\left\\{\begin{array}[]{ll}\hat{x}_{p}[k+1]&=A\hat{x}_{p}[k]+Bu[k]+Ly_{r}[k]\\\
y_{r}[k]&=y[k]-C\hat{x}_{p}[k],\;k=0,\dots,N_{h}-1.\end{array}\right.$ (9)

Here $A^{\Delta}\triangleq A+\Delta A(\delta)$ with $A$ representing the
nominal system matrix and $\delta\in\Omega$ denoting the probabilistic
parameter uncertainty with probability space
$(\Omega,\mathcal{D}_{a},\mathbf{P})$. We assume the uncertainty set
$\Omega\subset\mathbb{R}^{v}$ to be closed, bounded, and to include the zero
uncertainty yielding $\Delta A(0)=0$. The state of the process is represented
by $x_{p}[k]$, the output of the process is $y[k]$, $\tilde{u}[k]$ is the
control signal received by the process, $u[k]$ is the control signal generated
by the controller, $y_{p}[k]$ is the virtual performance output, and
$y_{r}[k]$ is the residue generated by the detector. In this paper, we assume
all signals have the same dimension $n_{x}$. That is, we consider a fully
actuated square system. The system is said to have a good performance over the
horizon $N_{h}$, when the energy of the performance output
($||y_{p}||_{\ell_{2},[0,N_{h}]}^{2}$) is small. In the closed-loop system
described above, we consider that an adversary is injecting false data into
the actuators. An attack is said to be detected when the energy of the
detection output ($||y_{r}||_{\ell_{2},[0,N_{h}]}^{2}$) is higher than a
predefined threshold (say $\epsilon_{r}$). We assume that the detection
threshold and the detector $L$ is designed to be robust against all
uncertainties.

Given this setup, we now discuss the resources the adversary has access to.
For clarity, we establish the following:

###### Assumption II.1

$(A^{\Delta},B)$ is controllable $\forall\delta\in\Omega$.
$\hfill\triangleleft$

###### Assumption II.2

Matrices $B$ and $C$ are invertible. $\hfill\triangleleft$

Figure 1: Control system under data injection attack on actuators

###  II-A Disruption and disclosure resources

The adversary can access the control channels and inject data. This is
represented by $\tilde{u}[k]={u}[k]+a[k],$ where $a[k]\in\mathbb{R}^{n_{x}}$
is the data injected by the adversary. The adversary cannot access the sensor
channels. The adversary does not have access to any disclosure (eavesdropping)
resources.

###  II-B System knowledge

We assume that, at design time, the defender knows the bounds of the set
$\Omega$ and that the system matrix $A^{\Delta}$ is known only up to the
nominal system matrix $A$. Next, at operation time, we assume that the
adversary has full system knowledge. That is, the adversary knows the system
matrix $A^{\Delta}$ without any uncertainties. In reality, it is hard for the
adversary to know the system matrices, but this assumption helps to study the
worst case.

The system knowledge is used by the adversary to calculate the optimal data
injection attacks. Defining $e[k]\triangleq x_{p}[k]-\hat{x}_{p}[k]$ and
${x}[k]\triangleq[x_{p}[k]^{T}\;e[k]^{T}]^{T}$, the closed-loop system under
attack with the performance output and detection output as system outputs
becomes

$\mathcal{P}_{cl}:\left\\{\begin{array}[]{ll}{x}[k+1]&={A}_{cl}^{\Delta}{x}[k]+{B}_{cl}a[k]\\\
y_{p}[k]&={C}_{p}{x}[k]\\\ y_{r}[k]&={C}_{r}{x}[k],\\\ \end{array}\right.$
(10) where
$\displaystyle{A}_{cl}^{\Delta}\triangleq\begin{bmatrix}A^{\Delta}+BKC&0\\\
\Delta A&A-LC\end{bmatrix},$
$\displaystyle{B}_{cl}\triangleq\begin{bmatrix}B\\\ B\end{bmatrix},$
$\displaystyle{C}_{p}\triangleq\begin{bmatrix}C_{J}&0\end{bmatrix},$
$\displaystyle{C}_{r}\triangleq\begin{bmatrix}0&C\end{bmatrix}.$

###  II-C Attack goals and constraints

Given the resources the adversary has access to, the adversary aims at
disrupting the system’s behavior whilst remaining stealthy. The system
disruption is evaluated by the increase in energy of the performance output,
and the attack signal is deemed to be stealthy when the energy of the
detection output is less than $\epsilon_{r}$. Next, we discuss the optimal
attack policy of the adversary and the design problem of the defender when the
system is deterministic.

###  II-D Design for a deterministic system

From the previous discussions, it can be understood that the goal of the
adversary is to maximize the energy of the performance output whilst remaining
stealthy. When the system is deterministic ($\Omega=\\{0\\}$), the attack
policy of the adversary can be formulated as

$\displaystyle q(K,0)\triangleq\sup_{a\in\ell_{2e}}$
$\displaystyle\|y_{p}(K,0)\|_{\ell_{2}}^{2}$ (11) s.t.
$\displaystyle\|y_{r}(K,0)\|_{\ell_{2}}^{2}\leq\epsilon_{r},\;x(K,0)[0]=0,$

where the subscript $[0,N_{h}]$ is dropped for clarity. In (11), $q(K,0)$ is
the disruption caused by the attack signal on the nominal system, $y_{p}(K,0)$
and $y_{r}(K,0)$ are the performance output and the detection output under the
given controller $K\in\mathbb{R}^{n_{x}\times n_{x}}$, and $N_{h}$ is the
horizon length. In (11), the constraint $x(K,0)[0]=0$ is introduced since the
system is at equilibrium before the attack commences.

###### Assumption II.3

The system (10) is at equilibrium before the attack commences.
$\hfill\triangleleft$

The aim of the defender then is to design a controller $K$ such that the
disruption caused by the adversary
($q(K,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}0})$)
is minimized. To this end, the design problem can be formulated as

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}K^{*}=\arg\inf_{K}q(K,0)}$
(12)

The design problem (12) is optimal only when (10) is deterministic. By
extending (12), we formulate the design problem when the system is uncertain
in the next section.

##  III Problem Formulation

Consider the data injection attack scenario where the parametric uncertainty
$\delta\in\Omega$ of the system is known to the adversary but not to the
defender. The defender knows only about the probabilistic description of the
set $\Omega$. In reality, it is hard for the adversary to know the system
matrices, but this assumption helps to study the worst case. Under this setup,
the adversary can cause high disruption by remaining stealthy as it will be
able to inject attacks by solving

$\displaystyle q(K,\delta)\triangleq\sup_{a\in\ell_{2e}}$
$\displaystyle\|y_{p}(K,\delta)\|_{\ell_{2}}^{2}$ (13) s.t.
$\displaystyle\|y_{r}(K,\delta)\|_{\ell_{2}}^{2}\leq\epsilon_{r},\;x(K,\delta)[0]=0,$

where $y_{p}({K,\delta})$ and $y_{r}({K,\delta})$ are the performance and
detection output corresponding to the controller $K$ and uncertainty $\delta$.
Since the defender does not know the system completely, $q(K,\delta)$ becomes
a random variable. Thus, from the defenders point of view, the best option is
to choose a feedback policy $K$, such that the risk corresponding to the
impact random variable $q(K,\delta)$ is minimized. This design problem can be
formulated as Problem 1.

###### Problem 1

Find an optimal feedback controller $K^{*}$ s.t:

$K^{*}\triangleq\arg\inf_{K}\mathcal{R}_{\Omega}(q(K,\delta)),$ (14)

where $\mathcal{R}_{\Omega}$ is a risk metric chosen by the defender. The
subscript $\Omega$ denotes that the risk acts on the uncertainty whose
probabilistic description is known to the defender. $\hfill\triangleleft$

Problem 1 searches for a controller $K$ such that the risk is minimized. Let
us consider the setup where the defender evaluates the risk based on the risk
metric CVaR. CVaR is used in the research community due to its numerous
advantages [15] and is defined in Definition III.1.

######  Definition III.1 (CVaR [13])

Given a random variable $X$ and $\alpha\in(0,1)$, the CVaR is defined as
111This Definition assumes the distribution of $X$ has no point masses. For
general definitions of CVaR see [16].

$\text{CVaR}_{\alpha}(X)=\mathbb{E}\\{X|X>\text{VaR}_{\alpha}(X)\\},$
$\text{where}\quad\quad\text{VaR}_{\alpha}(X)=\inf\\{x|\mathbb{P}[X\geq
x]\leq\alpha\\}.$

$\text{CVaR}_{\alpha}(X)=\beta$ implies that $X\leq\beta$ at least
$\alpha\times 100\%$ of the time on average. $\hfill\triangleleft$

In our setting, the defender is interested in determining the controller such
that the $\text{CVaR}_{\alpha}$ (given $\alpha$) of the impact random variable
($q(K,\delta)$) is minimized. To this end, Problem 1 can be reformulated as

$K^{*}=\arg\inf_{K}\mathbb{E}_{\Omega}\\{q(K,\delta)|q(K,\delta)>\text{VaR}_{\alpha}(q(K,\delta))\\}.$
(15)

There are two difficulties in solving (15). Firstly, (13) is non-convex for
any given $\delta$. Secondly, since the operator $\mathbb{E}$ operates over
the continuous space $\Omega$, the optimization problem (15) is
computationally intensive and in general NP-hard. To this end, in Section IV,
we determine a convex approximation for (13). We then use this approximation,
to recast (15) as a convex optimization problem. In Section V, we provide a
method to approximate the expectation operator.

##  IV Design problem formulation using a convex Impact proxy

In this section, we consider the function $q(K,\delta_{j})$ for a given
uncertainty $\delta_{j}\in\Omega$ and prove that it has an upper bound. We
then show that the term of the upper bound that is dependent on the controller
(say $\bar{q}(\cdot)$) is convex in $K$. The main objective of performing this
step is that, once we determine the term $\bar{q}(\cdot)$, it can be used in
(15) instead of $q(K,\delta_{j})$ to formulate a relaxed convex design
problem. To this end, we will refer to $\bar{q}(\cdot)$ as Impact proxy in the
reminder of the paper.

To derive the upper bound, we begin by defining the vectors
$\textbf{a}_{j}\triangleq\begin{bmatrix}a_{j}[0]^{T},\dots,a_{j}[N_{h}-1]^{T}\end{bmatrix}^{T}$,
$\textbf{x}_{p,j}\triangleq\begin{bmatrix}x_{p,j}[1]^{T},\dots,x_{p,j}[N_{h}]^{T}\end{bmatrix}^{T}$,
$\textbf{e}_{j}\triangleq\begin{bmatrix}e_{j}[1]^{T},\dots,e_{j}[N_{h}]^{T}\end{bmatrix}^{T}$,
$\textbf{y}_{p,j}\triangleq\begin{bmatrix}y_{p,j}[1]^{T},\dots,y_{p,j}[N_{h}]^{T}\end{bmatrix}^{T}$,
and
$\textbf{y}_{r,j}\triangleq\begin{bmatrix}y_{r,j}[1]^{T},\dots,y_{r,j}[N_{h}]^{T}\end{bmatrix}^{T}$.
Here $\textbf{a}_{j},\textbf{x}_{p,j},\textbf{y}_{p,j}$ and $\textbf{y}_{r,j}$
are the stacked attack vector, system state, performance output vector and the
detection output vectors corresponding to the uncertainty $\delta_{j}$
respectively. Let us define the matrices
$F_{xa}(K,\delta_{j}),F_{ea}(\delta_{j}),F_{ex}(\delta_{j})\in\mathbb{R}^{n_{x}N_{h}\times
n_{x}N_{h}}$, such that

$\displaystyle\textbf{x}_{p,j}$
$\displaystyle=F_{xa}(K,\delta_{j})\textbf{a}_{j},\;\textbf{e}_{j}=F_{ea}(\delta_{j})\textbf{a}_{j}+F_{ex}(\delta_{j})\textbf{x}_{p,j},$
(16) $\displaystyle\textbf{y}_{p,j}$
$\displaystyle=F_{p}(K,\delta_{j})\textbf{a}_{j},\;\;\textbf{y}_{r,j}=F_{r}(K,\delta_{j})\textbf{a}_{j},$
(17) $\displaystyle F_{p}$
$\displaystyle(K,\delta_{j})\triangleq(I_{N_{h}}\otimes
C_{J})F_{xa}(K,\delta_{j}),$ (18) $\displaystyle F_{r}$
$\displaystyle(K,\delta_{j})\triangleq(I_{N_{h}}\otimes
C)(F_{ea}(\delta_{j})+F_{ex}(\delta_{j})F_{xa}(K,\delta_{j})).$ (19)

Under the uncertainty $\delta_{j}$, let us represent the system matrix of (4)
by $A_{j}$. Then $F_{\upsilon a}(K,\delta_{j}),\upsilon=\\{x,e\\}$ is given by

$\begin{bmatrix}B&0&\dots&0\\\ {A}_{\upsilon,j}B&B&\dots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
{A}_{\upsilon,j}^{N_{h}-1}B&{A}_{\upsilon,j}^{N_{h}-2}B&\dots&B\end{bmatrix},$

where ${A}_{x,j}\triangleq A_{j}+BKC$ and $A_{e,j}\triangleq A_{j}-LC$.
Similarly $F_{ex}(\delta_{j})$ is given by

$\displaystyle\begin{bmatrix}0&0&\dots&0\\\ \Delta A&0&\dots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ {A}_{e,j}^{N_{h}-2}\Delta
A&{A}_{e,j}^{N_{h}-3}B\Delta A&\dots&0\end{bmatrix}$ (20)

Under these definitions, $q(K,\delta_{j})$ can be obtained by the non-convex
optimization problem

$\displaystyle q(K,\delta_{j})\triangleq\sup_{\textbf{a}_{j}}$
$\displaystyle\|F_{p}(K,\delta_{j})\textbf{a}_{j}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\|_{2}^{2}}$
(21) s.t.
$\displaystyle\|F_{r}(K,\delta_{j})\textbf{a}_{j}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\|_{2}^{2}}\leq\epsilon_{r},x(K,\delta_{j})[0]=0.$

Next, we derive the upper bound of $q(K,\delta_{j})$ in Lemma IV.1.

###### Lemma IV.1

Let $\delta_{j}\in\Omega$ and $\kappa\triangleq F_{p}(\cdot)F_{r}^{-1}(\cdot)$
(Here the arguments of $F_{p}$ and $F_{r}$ are dropped for clarity). Let the
matrix $B$ be invertible. Then, it holds that

$q(K,\delta_{j})\leq\mu\bar{q}(K,\delta_{j})^{f},\;\bar{q}(K,\delta_{j})\triangleq\eta||K||_{F}^{2}+\sum_{i=2}^{n_{x}N_{h}}\sigma_{i}(\kappa^{-1}),$
(22)

where $f\triangleq n_{x}N_{h}-1,\mu$ is a term independent of $K$ and $\eta$
is a positive scalar weight on the regularization term.

###### Proof:

See Appendix. ∎

In Lemma IV.1, we formulated an upper bound for ${q}(K,\delta_{j})$. However,
only the term $\bar{q}(K,\delta_{j})$ of the bound is dependent on the
variable $K$. Moreover, since
$\left(\bar{q}(K,\delta_{j})\right)^{n_{x}N_{h}-1}$ is a monotonically
increasing function on $\bar{q}(K,\delta_{j})>0$, $\bar{q}(K,\delta_{j})$ can
be replaced as the term to be optimized. Next, we show that
$\bar{q}(K,\delta_{j})$ is strongly convex in the design variable $K$.

###### Theorem IV.2

For any given $\delta_{j}\in\Omega$, the function $\bar{q}(K,\delta_{j})$ is
strongly convex in the design variable $K$.

###### Proof:

See Appendix. ∎

We have shown in this section that the term $\bar{q}(K,\delta_{j})$ can be
used as the convex proxy objective function for the attack impact
${q}(K,\delta_{j})$. That is, we can recast (15) as

$K^{*}=\arg\inf_{K}\mathbb{E}\\{\bar{q}(K,\delta)|\bar{q}(K,\delta)>\text{VaR}_{\alpha}(\bar{q}(K,\delta))\\}.$
(23)

Although (23) is convex, it is computationally intensive due to the
expectation operator. In Section V, we discuss a method to approximate the
expectation operator.

###### Remark 1

We use Definition III.1 to formulate (23). Thus (23) implicitly assumes that
the the distribution of $\bar{q}$ has no point masses. However, verifying this
conditions is beyond the scope of this paper and is left for future work.

##  V Empirical risk using scenario based approach

The optimization problem (23) is computationally intensive since it involves
an expectation operator which acts on a continuum of uncertainties $\Omega$.
In this section, we provide a method to approximate the expectation operator
using the scenario-based approach [13]. To begin with, let us establish the
following:

###### Assumption V.1

For any $\delta$, $\bar{q}(\cdot,\delta)$ is a convex function in the design
variable $K$. $\hfill\triangleleft$

We have shown in Theorem IV.2 that Assumption V.1 is satisfied. Next, we
approximate the expectation operator in (23) empirically. To do this, let us
begin by sampling the uncertainty set $\Omega$ with $N$ samples. Let us
consider that $(\delta_{1},\dots,\delta_{N})$ is a collection of $N$
independent realizations from $\Omega$ and let
$\Omega_{N}\triangleq\\{1,\dots,N\\}$. Then for any given $K$, and
$i\in\Omega_{N}$, we denote by $\bar{q}_{i}(K)$, the value attained by
$\bar{q}(K,\delta_{i})$, and we denote by $\bar{q}_{(i)}(K)$, the $N-i+1^{th}$
order statistic. That is
$\bar{q}_{(1)}(K)\geq\bar{q}_{(2)}(K)\geq\dots\geq\bar{q}_{(N)}(K)$. Now we
present the first result of the section.

###### Lemma V.1

Let the dimension of the design variable $K$ be $d=n_{x}^{2}$. Let $N\geq d$
and $m\triangleq\lceil N(1-\alpha)\rceil$. Then, under Assumption V.1, the
solution to (23) can be obtained empirically by solving the convex
optimization problem

$K^{*}=\arg\inf_{K}\frac{1}{m}\sum_{i=1}^{m}\bar{q}_{(i)}(K).$ (24)

###### Proof:

It was shown that (24) is the empirical version of (23) in [13] (See equations
$(3)$ and $(5)$ in [13]). ∎

The design problem (24) is formulated using the empirical formulation of the
risk metric CVaR. If we could order the functions $\bar{q}_{(i)}(K)$, then the
problem of obtaining the optimal controller would be simple. However,
$\bar{q}_{(i)}(K)$ is a function of the optimization variable $K$. Hence it is
not possible to explicitly define the order without knowing $K$ beforehand.
Alternatively, the design problem (24) can be modified such that the
controller is optimal to any possible ordering of the functions. To this end,
(24) can be recast as

$\displaystyle\arg\inf_{K,y}$ $\displaystyle y$ (25) s.t.
$\displaystyle\frac{1}{m}\sum_{j=1}^{m}\bar{q}_{i_{j}}(K)\leq y,$
$\displaystyle\forall\text{choice of $m$
indices}\\{i_{1},\dots,i_{m}\\}\subseteq\Omega_{N}.$

Since $\\{i_{1},\dots,i_{m}\\}$ is any subset of $\Omega_{N}$ with cardinality
$m$, the optimization problem (25) has ${N\choose m}$ constraints. To
summarize, (25) is the empirical reformulation of (23).

Using the solution obtained from (25), we provide probabilistic out-of-sample
certificate on the random variable $\bar{q}(K,\delta),\delta\in\Omega$. To
this end, let us represent the optimal controller of (25) as $K_{N}^{*}$, and
let $N\geq m+d$. Let us define the out-of-sample Probability of Shortfall (PS)
as follows:

######  Definition V.1 (Probability of Shortfall)

PS is defined as

$PS(K_{N}^{*})\triangleq\mathbf{P}\\{\delta\in\Omega\;|\;\bar{q}(K_{N}^{*},\delta)\geq\bar{q}_{(m+d)}(K_{N}^{*})\\}.$
$\triangleleft$

By ensuring that the $PS(K_{N}^{*})$ is small (say $\epsilon$), one can ensure
that the impact proxy under the optimal controller for any new uncertainty
drawn from the set $\Omega$, $\bar{q}(K_{N}^{*},\delta)$, exceeds a predefined
valued $\bar{q}_{(m+d)}(K_{N}^{*})$ (also called as the shortfall threshold)
with a small probability $\epsilon$. Now we are ready to present the main
result in Theorem V.2.

###### Theorem V.2

It holds that $\mathbf{P}^{N}\\{PS(K_{N}^{*})\leq\epsilon\\}\triangleq$

$\int_{0}^{\epsilon}\frac{\Gamma(N+1)}{\Gamma(m+d)\Gamma(N+1-m-d)}p^{m+d-1}(1-p)^{N-m-d}dp,$
(26)

where $\Gamma$ is Euler’s Gamma function and $\epsilon\in(0,\;1).$

###### Proof:

See Appendix. ∎

Theorem V.2 provides posteriori results on the confidence with which
$PS(K_{N}^{*})$ is below a small threshold $\epsilon$. In other words, Theorem
V.2 states that, the confidence of the $PS(K_{N}^{*})$ can be evaluated by
knowing the dimension of the decision variable ($d$), the number of samples
($N$), and $m$.

To recall, in this section we proposed an empirical version of (23) in (25).
We also provided probabilistic guarantees on the out-of-sample PS. We conclude
this section by providing Algorithm 1 which depicts the outline for solving
Problem 1 approximately. In the next section, we depict the efficacy of the
proposed algorithm using a numerical example.

Result: $K_{N}^{*},PS$

Initialization: $N_{h},k,d,\Omega,\epsilon$;

  1. 1.

Choose $N$ such that $N\geq m+d$.

  2. 2.

Extract $N$ independent realizations from $\Omega$.

  3. 3.

Build the matrix $\kappa^{-1},\forall i\in\Omega_{N}$.

  4. 4.

Solve the convex optimization problem (25).

  5. 5.

Given $\epsilon$, evaluate the confidence using Theorem V.2.

  6.

Algorithm 1 Risk averse design algorithm

##  VI Numerical example

In this section, the efficacy of the design algorithm is depicted through a
numerical example. Consider the system of the form (10) where $C_{J}=C=I_{3}$,
$A,B$ and $L$ are

$\begin{bmatrix}2&0&1\\\ 1&a&0\\\ 0&1&b\end{bmatrix},\begin{bmatrix}1&1&0\\\
0&0.3&1\\\ 0&0&1\end{bmatrix},\text{and}\;\begin{bmatrix}1.95&0&1\\\
1&0.36&1\\\ 0&1&-0.87\end{bmatrix},$

respectively. Here $a\in[0.5\;1.5]$ and $b\in[-0.5\;0.5]$ are uncertain
parameters. The observer gain $L$ is designed using pole placement method.

|  CVaR${}_{\alpha}(\bar{q}(K,\delta)+0.1||K||_{F}^{2})$ | $\bar{q}_{(m+d)}(K_{N}^{*})$
---|---|---
$m=1$ | 20.7160 | 16.9069
$m=2$ | 20.6436 | 16.8910
TABLE I: Risk and Shortfall threshold  Figure 2: Evaluation of controller
performance: The plot depicts the distribution of the attack impact $q(\cdot)$
and the impact proxy $\bar{q}(\cdot)$ for $100$ different uncertainties under
the nominal and the optimal controller. The optimal controller is obtained
from Algorithm 1 by optimizing CVaR${}_{0.8}(\bar{q}(\cdot))$ with parameters
$N_{h}=5,N=11$, and $m=2$. The nominal controller is obtained by optimizing
$\bar{q}(K,0)$ for the nominal system. On each box, the central mark indicates
the median, and the bottom and top edges of the box indicate the $25^{th}$ and
$75^{th}$ percentiles, respectively. The whiskers extend to the most extreme
data points.

Let $N=11,\epsilon_{r}=1$, and $N_{h}=5$. The risk and the corresponding
shortfall threshold $\bar{q}_{(m+d)}(K_{N}^{*})$ obtained from Algorithm 1 for
different values of $m$ is shown in TABLE I. The corresponding confidence that
the the out-of-sample PS is less than or equal to $\epsilon$ can be evaluated
from the integral in Theorem V.2. In reality, the confidence can be made
higher by increasing $N$ or decreasing $m$.

Next, we depict the efficiency of the proposed design framework using Fig. 2.
It can be seen that: (a) compared to the nominal controller, the controller
obtained from Algorithm 1 lowers the impact proxy $\bar{q}(\cdot)$); (b) The
value of attack impact ($q(\cdot)$) is shifted towards a lower median.

##  VII Conclusion

The problem of risk-optimal controller design against FDI attacks on actuators
of an uncertain control system was studied. We considered an attacker with
perfect system knowledge that maximizes the system disruption whilst remaining
stealthy. We quantified the system disruption using the OOG. We formulate a
design problem where the defender aims at designing a controller such that its
CVaR is minimized. We proposed a convex optimization problem that
approximately solves the design problem with probabilistic certificates.
Finally, we illustrated the results through a numerical example.

## Proof of Lemma IV.1

Before presenting the proof, we provide an intermediate result which helps to
construct the proof.

######  Lemma A.1 (Inequality of arithmetic and geometric means)

Given $N$ real numbers $x_{1},\dots,x_{N}$, it holds that
$\prod_{i=1}^{N}x_{i}\leq\big{(}\frac{\sum_{i=1}^{N}x_{i}}{N}\big{)}^{N}$

###### Proof:

The convex dual problem of (21) can be formulated as (27). Furthermore, it was
shown in [17, Theorem 3.1] that the duality gap is zero.

$\inf\\{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\epsilon_{r}}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{}\gamma\;|\;F_{p}^{T}(\cdot)F_{p}(\cdot)-\gamma
F_{r}^{T}(\cdot)F_{r}(\cdot)\preceq 0.\\}$ (27)

Pre-multiplying the constraint of (27) by $F_{r}^{-T}(\delta_{j})$ and post-
multiplying by $F_{r}^{-1}(\delta_{j})$, (27) can be rewritten as

$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\epsilon_{r}}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{}\inf\\{\gamma\;|\;\kappa^{T}\kappa\preceq\gamma
I.\\}$ (28)

Since $\epsilon_{r}$ is a pre-defined constant, it can be moved outside the
optimization problem (28). Using the definition of singular values, (28) can
be re-written as
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\epsilon_{r}\bar{\sigma}(\kappa)}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{}$.
Thus we have shown that
$q(K,\delta_{j})=\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\epsilon_{r}\bar{\sigma}(\kappa)}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{}$.
Next, we prove that
$\bar{\sigma}(\kappa)\leq\mu\bar{q}(K,\delta_{j})^{n_{x}N_{h}-1}$. To this
end, we can show that the matrix $F_{r}(\cdot)$ is a block lower triangular
with the element $CB$ in the leading diagonal. Then $F_{r}(\delta_{j})^{-1}$
will be $\begin{bmatrix}B^{-1}C^{-1}&\dots&0\\\ \vdots&\ddots&\vdots\\\
*&\dots&B^{-1}C^{-1}\end{bmatrix},$ where $*$ represents that its value is
unimportant for now. Then, matrix $\kappa$ will be of the form
$\begin{bmatrix}C_{J}C^{-1}&\dots&0\\\ \vdots&\ddots&\vdots\\\
*&\dots&C_{J}C^{-1}\end{bmatrix}.$ Since $\kappa$ is block lower triangular,
its determinant is the product of determinant of diagonal blocks [12, Proof of
Theorem A.1]. Thus $\det(\kappa)=|\det(C_{J}C^{-1})^{N_{h}}|\triangleq c$.
Since the product of singular values of a matrix is equal to its determinant,
we get

$\bar{\sigma}(\kappa)=\frac{c}{\prod_{i=1}^{n_{x}N_{h}-1}\sigma_{i}(\kappa)}\stackrel{{\scriptstyle
1}}{{=}}c\prod_{i=2}^{n_{x}N_{h}}\sigma_{i}(\kappa^{-1}).$ (29)

The equality $1$ in (29) follows since the singular values of $\kappa$ and
$\kappa^{-1}$ are reciprocals of each other. By using the result of Lemma A.1,
the term $\epsilon_{r}\bar{\sigma}(\kappa)$ can be bounded as

$\displaystyle\epsilon_{r}\bar{\sigma}(\kappa)$
$\displaystyle\leq\epsilon_{r}|\det(C_{J}C)^{N_{h}}|\Big{(}\frac{\sum_{i=2}^{n_{x}N_{h}}\sigma_{i}(\kappa^{-1})}{n_{x}N_{h}-1}\Big{)}^{n_{x}N_{h}-1}$
(30)
$\displaystyle=\underbrace{\frac{\epsilon_{r}|\det(C_{J}C)^{N_{h}}|}{(n_{x}N_{h}-1)^{n_{x}N_{h}-1}}}_{\mu}\Big{(}\sum_{i=2}^{n_{x}N_{h}}\sigma_{i}(\kappa^{-1})\Big{)}^{n_{x}N_{h}-1}$
(31)
$\displaystyle\leq\mu\Big{(}\underbrace{\sum_{i=2}^{n_{x}N_{h}}\sigma_{i}(\kappa^{-1})+\eta||K||_{F}^{2}}_{\bar{q}(K,\delta_{j})}\Big{)}^{n_{x}N_{h}-1}$
(32)

where the last inequality follows since $\mu$ and $\eta||K||_{F}^{2}$ are non-
negative terms. In (32), only the term $\bar{q}(K,\delta_{j})$ is dependent on
the design variable $K$. This concludes the proof.∎

## Proof of Theorem IV.2

###### Proof:

Since $\kappa=F_{p}(\cdot)F_{r}(\cdot)^{-1}$, it follows that $\kappa^{-1}=$

$F_{r}(K,\delta_{j})F_{xa}(K,\delta_{j})^{-1}(I_{N_{h}}\otimes C_{J})^{-1}$
(33) $=(I_{N_{h}}\otimes
C)(F_{ea}(\delta_{j})F_{xa}(K,\delta_{j})^{-1}+F_{ex}(\delta_{j}))(I_{N_{h}}\otimes
C_{J})^{-1}.$ (34)

Here all the matrices are independent of the design variable $K$ except
$F_{xa}(K,\delta_{j})^{-1}$. It can be verified by matrix multiplication that
$F_{xa}(K,\delta_{j})^{-1}$ is of the form

$\begin{bmatrix}B^{-1}&0&0&\dots&0\\\ -B^{-1}{A}_{x,j}&B^{-1}&0&\dots&0\\\
0&-B^{-1}{A}_{x,j}&B^{-1}&\dots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&0&\dots&B^{-1}\end{bmatrix}$ (35)

The matrix ${A}_{x,j}$ is affine in $K$ and thus the same holds for
$\kappa^{-1}$. It can also be shown that the sum of singular values of a
matrix, $X\to\sum\sigma(X)$, is convex [18]. THus we have proven that the term
$\sum\sigma(\kappa^{-1})$ is convex in $K$. Since the regularization term
$||K||_{F}^{2}$ is strictly convex in $K$, we then conclude that
$\bar{q}_{i_{j}}(K,\delta)+\eta||K||_{F}^{2}$ is a strictly convex function in
$K$ which concludes the proof. ∎

## Proof of Theorem V.2

Before presenting the proof, we provide an intermediate result Theorem A.2:
which follows from applying [13, Theorem 3.1] to our problem setup.

###### Theorem A.2

Let us suppose that $(a)$ a solution to (25) exists and it is unique almost
surely, and $(b)$ for a sample $(\delta_{1},\dots,\delta_{N})$ of independent
realizations from $\Omega$, the event that {
$\exists\;K|\bar{q}_{i}(K),\forall\;i\in\Omega_{N}$ has the same value} has
zero probability. Then $\mathbf{P}^{N}\\{PS(K_{N}^{*})\leq\epsilon\\}$ has a
Beta($m+d,N+1-m-d$) distribution.

That is, Theorem A.2 states that the result of Theorem V.2 follows if $(a)$
and $(b)$ hold: which we prove next.

###### Proof:

Proof that $(a)$ holds: We intend to show that the (25) has a unique solution.
However, as discussed in the paper, (25) is simply a reformulation of (24).
Thus, we can equivalently show that the solution to (24) is unique.

In Theorem IV.2, we have shown that the term $\bar{q}_{i_{j}}(K)$ is strongly
convex function in $K$. The uniqueness of the optimal solution, $K_{N}^{*}$,
then follows from the strict convexity of the objective function.
Proof that $(b)$ holds: We prove by contradiction. Let us assume that $\exists
K,l$ and independent samples $(\delta_{1},\dots,\delta_{N})$ such that
$\textbf{P}^{N}\\{\bar{q}_{i}(K)=l,i=1,2,\dots,N\\}=\theta\geq 0$. Then by
definition, $\exists\;\delta_{j}$ such that
$\textbf{P}\\{\delta_{j}|\bar{q}_{j}(K)=l\\}\neq 0$. However this contradicts
the assumption that the distribution of $\bar{q}(\cdot)$ has no point masses.
This concludes the proof. ∎

## References

  * [1] R. Langner, “Stuxnet: Dissecting a cyberwarfare weapon,” IEEE Security & Privacy, vol. 9, no. 3, pp. 49–51, 2011.
  * [2] A. Cherepanov and R. Lipovsky, “Industroyer: Biggest threat to industrial control systems since stuxnet,” WeLiveSecurity, ESET, vol. 12, 2017.
  * [3] D. U. Case, “Analysis of the cyber attack on the ukrainian power grid,” Electricity Inf. Sharing and Analy. Center (E-ISAC), vol. 388, 2016.
  * [4] S. M. Dibaji, M. Pirani, D. B. Flamholz, A. M. Annaswamy, K. H. Johansson, and A. Chakrabortty, “A systems and control perspective of cps security,” Annu. Reviews in Contr., vol. 47, pp. 394–411, 2019.
  * [5] P. Griffioen, S. Weerakkody, B. Sinopoli, O. Ozel, and Y. Mo, “A tutorial on detecting security attacks on cyber-physical systems,” in 2019 18th Eur. Contr. Conf. (ECC), pp. 979–984, IEEE, 2019.
  * [6] M. I. Müller, J. Milošević, H. Sandberg, and C. R. Rojas, “A risk-theoretical approach to $\mathcal{H}_{2}$-optimal control under covert attacks,” in 2018 IEEE Conf. on Decision and Contr. (CDC), pp. 4553–4558, IEEE, 2018.
  * [7] C. Murguia, I. Shames, J. Ruths, and D. Nešić, “Security metrics and synthesis of secure control systems,” Automatica, vol. 115, p. 108757, 2020\.
  * [8] N. Hashemi and J. Ruths, “Gain design via LMIs to minimize the impact of stealthy attacks,” in 2020 Am. Contr. Conf. (ACC), pp. 1274–1279, IEEE, 2020.
  * [9] S. C. Anand and A. M. Teixeira, “Joint controller and detector design against data injection attacks on actuators,” IFAC-PapersOnLine, vol. 53, no. 2, pp. 7439–7445, 2020.
  * [10] S. D. Bopardikar, A. Speranzon, and J. P. Hespanha, “An ${H}_{\infty}$ approach to stealth-resilient control design,” in 2016 Resilience Week (RWS), pp. 56–61, IEEE, 2016.
  * [11] A. Teixeira, H. Sandberg, and K. H. Johansson, “Strategic stealthy attacks: the output-to-output $\ell_{2}$-gain,” in 2015 54th IEEE Conf. on Decision and Contr. (CDC), pp. 2582–2587, IEEE, 2015.
  * [12] K. Dvijotham, E. Todorov, and M. Fazel, “Convex structured controller design in finite horizon,” IEEE Trans. on Contr. of Network Systems, vol. 2, no. 1, pp. 1–10, 2014.
  * [13] F. A. Ramponi and M. C. Campi, “Expected shortfall: Heuristics and certificates,” Eur. Journal of Operational Research, vol. 267, no. 3, pp. 1003–1013, 2018.
  * [14] M. C. Campi and S. Garatti, Introduction to the scenario approach.  SIAM, 2018.
  * [15] R. T. Rockafellar, S. Uryasev, et al., “Optimization of conditional value-at-risk,” Journal of risk, vol. 2, pp. 21–42, 2000.
  * [16] R. T. Rockafellar and S. Uryasev, “Conditional value-at-risk for general loss distributions,” Journal of banking & finance, vol. 26, no. 7, pp. 1443–1471, 2002.
  * [17] S. C. Anand and A. M. H. Teixeira, “Stealthy cyber-attack design using dynamic programming,” in 2021 60th IEEE Conference on Decision and Control (CDC), pp. 3474–3479, 2021.
  * [18] S. Subramani, On sums of singular values.  University of New South Wales, 1993.


# LINNA: Likelihood Inference Neural Network Accelerator

Chun-Hao To 11footnotetext: Corresponding author. Eduardo Rozo  Elisabeth
Krause  Hao-Yi Wu  Risa H. Wechsler  Andrés N. Salcedo

###### Abstract

Bayesian posterior inference of modern multi-probe cosmological analyses
incurs massive computational costs. For instance, depending on the
combinations of probes, a single posterior inference for the Dark Energy
Survey (DES) data had a wall-clock time that ranged from 1 to 21 days using a
state-of-the-art computing cluster with 100 cores. These computational costs
have severe environmental impacts and the long wall-clock time slows
scientific productivity. To address these difficulties, we introduce LINNA:
the Likelihood Inference Neural Network Accelerator. Relative to the baseline
DES analyses, LINNA reduces the computational cost associated with posterior
inference by a factor of 8–50. If applied to the first-year cosmological
analysis of Rubin Observatory’s Legacy Survey of Space and Time (LSST Y1), we
conservatively estimate that LINNA will save more than US $\$300,000$ on
energy costs, while simultaneously reducing $\rm{CO}_{2}$ emission by $2,400$
tons. To accomplish these reductions, LINNA automatically builds training data
sets, creates neural network surrogate models, and produces a Markov chain
that samples the posterior. We explicitly verify that LINNA accurately
reproduces the first-year DES (DES Y1) cosmological constraints derived from a
variety of different data vectors with our default code settings, without
needing to retune the algorithm every time. Further, we find that LINNA is
sufficient for enabling accurate and efficient sampling for LSST Y10 multi-
probe analyses. We make LINNA publicly available at
https://github.com/chto/linna, to enable others to perform fast and accurate
posterior inference in contemporary cosmological analyses.

##  1 Introduction

Modern cosmological analyses typically employ a Bayesian framework, where
posteriors are sampled using Markov Chain Monte Carlo (MCMC) techniques.
Despite significant progress in the development of MCMC sampling algorithms
[e.g. 1, 2], the computational requirements associated with performing these
calculations can be significant. This is especially true for multi-probe
analyses [3, 4, 5, 6, 7], which typically require the introduction of many
tens of nuisance parameters to account for systematic effects presented in the
data. As a specific example, the joint analysis of clusters, weak lensing, and
galaxy clustering of [7] required 26 nuisance parameters. To achieve a well-
sampled posterior, the chains in that analysis required at least five million
likelihood evaluations, totaling $\sim 50k$ CPU hours. That is, a single run
had a wall-clock time of three weeks using a state-of-the-art computing
cluster with 100 cores.

The primary difficulty in this type of analysis is that the theory model used
to evaluate the likelihood is computationally expensive. One way to bypass
this limitation is to rely on fast surrogates of said model [e.g. 8, 9, 10,
11, 12]. However, the surrogates are typically built in pre-defined parameter
spaces that are narrower than those used in most cosmological analyses. That
is, the use of surrogates that are constructed a priori either limits the
choice of priors that can be employed, or requires extrapolations of the
surrogate outside the parameter regions over which the surrogate was
calibrated.

To address these shortcomings, as well as to avoid having to build surrogates
of theory models on a case-by-case basis, we present a new posterior inference
tool: the Likelihood Inference Neural Network Accelerator (LINNA). LINNA is
designed to enable fast and accurate posterior inference with arbitrary priors
through the use of automatically built surrogates of the theory model. Given
the priors of the parameters, a theory model, the corresponding posterior
function, and an efficient MCMC sampling algorithm, LINNA automatically
generates training data, trains a neural network emulator of the model
predictions, and samples the posterior using the provided MCMC sampling
algorithm. The training data for the neural network emulator is generated in
an iterative fashion designed to achieve accurate posteriors.

In this paper, we demonstrate that LINNA succeeds in accurately recovering the
parameter posteriors derived from three different data vectors derived from
the Dark Energy Survey [13] year 1 analyses (DES Y1):

  1. 1.

3$\times$2pt, a joint analysis of galaxy clustering, galaxy–galaxy lensing,
and cosmic shear [3].

  2. 2.

4$\times$2pt+N, a joint analysis of cluster–galaxy cross correlations, cluster
lensing, cluster clustering, and cluster abundances [7, 14].

  3. 3.

6$\times$2pt+N, a joint analysis of data vectors in 3$\times$2pt and
4$\times$2pt+N. [7]

Further, we explicitly demonstrate that LINNA accurately reproduces the
forecasted constraints for Rubin Observatory’s Legacy Survey of Space and Time
(LSST) Year 10 data set (see appendix C for details).

This paper is organized as follows: In section 2, we describe the design of
LINNA, the architecture of the neural network, the generation of training
data, and the optimization of neural network parameters. In section 3, we
demonstrate that LINNA can properly sample the posteriors of 6$\times$2pt+N
(section 3.1), 4$\times$2pt+N (section 3.2), and 3$\times$2pt (section 3.2).
We further compare the performance of LINNA to a promising fast sampler
Bayesfast. The amount of reduction on $\rm{CO}_{2}$ emission by using LINNA is
estimated in section 3.3. A concluding remark is presented in section 4.

##  2 Methods

###  2.1 LINNA design

Figure 1: Outline of LINNA structure. Green blocks indicate inputs provided by
users and orange blocks represent outputs of LINNA. Blocks within the blue
rectangular are run iteratively. The symbol $\bigoplus$ represents the
procedure of expanding the training data by randomly sampling the MCMC chain
from the previous iteration.

LINNA uses neural networks to create fast surrogates of computationally
expensive theory models. This emulator can then be used in combination with
efficient posterior sampling algorithms to quickly recover accurate
approximations of the parameter posteriors. There are three main reasons why
this approach is attractive:

  1. 1.

Within the context of posterior estimation, the parameter space of interest
has a well-defined boundary: the prior. This boundary makes it possible to
generate a training set that spans the entire parameter space. This ensures
that the LINNA emulator never extrapolates outside the domain over which it
has been trained.

  2. 2.

The neural network surrogate is differentiable, allowing the use of more
efficient sampling algorithms, such as the Hamiltonian Monte Carlo (HMC). In
practice, existing HMC codes are not parallelized, which incurs a significant
penalty in the wall-clock time required to generate an MCMC chain relative to
less efficient but fully parallelized sampling schemes.

  3. 3.

Because the posteriors are sampled using standard techniques, we may readily
use well-tested convergence tests of the resulting chains.

Despite these features, a significant problem remains. In high-dimensional
parameter spaces, the volume of reasonably defined priors can greatly exceed
the volume of the posteriors. Consequently, uniform sampling of the prior is
extremely inefficient: it is likely that only a few training points fall
within one- to three- sigma region of the posteriors, where high accuracy of
the neural network surrogates is most desirable. We solve this problem using
an iterative approach. In particular, a crude neural network is trained on a
training set that uniformly samples the priors. This neural network is then
fed into a posterior sampling tool to produce a MCMC chain, from which
additional training data is generated. In this way, we construct a training
set that covers the entire prior volume and is also capable of accurately
reconstructing the parameter posteriors with comparatively little training
data.

In the first few iterations, the neural network surrogates are rough at best,
so their high-confidence regions can be offset from those of the true
posterior. To overcome this problem, for the first few iterations of the
algorithm, we enlarge the parameter posteriors using a temperature parameter
$T$ such that

$\chi^{2}_{T}=\chi^{2}_{\rm{NN}}/T^{2},$ (2.1)

where $\chi^{2}_{\rm{NN}}$ is the $\chi^{2}$ value of the data calculated
using neural network surrogates. We can think of this scaling as simply
increasing all observational errors by a factor of T. The resulting broadening
allows for the posteriors of the neural network surrogate model to be offset
from the true posterior while still completely covering the same. In this way,
when we generate new training data from the surrogate posterior we are able to
improve the performance of the network over the true parameter posterior, even
if the latter was offset from that of the approximate model. Initially, the
temperature $T$ should be set to a large value to ensure that the training
data spans a wide enough parameter space. In the subsequent iterations, $T$
can be gradually decreased due to the increasing accuracy of the neural
network model over the region of interest. Empirically, we find that setting
$T=4$ and $2$ for the first two iterations and $T=1$ for all subsequent
iterations leads to good posteriors.

We summarize the workflow for LINNA in figure 1. First, we use Latin Hypercube
sampling of the prior volume space to generate an initial set of training
points. At each training point, we use the provided theory model to evaluate
the expectation value of the observable. The resulting training points and
model expectations are used to train a neural network to arrive at our first
surrogate model. This surrogate model in combination with the provided
likelihood function, which is scaled by the appropriate temperature parameter,
is fed to a standard sampling algorithm to arrive at MCMC samples of the
posterior. The resulting chain is then randomly sampled to expand the number
of points in the parameter space at which the neural network is trained. The
procedure is iterated, decreasing the temperature from iteration to iteration
as described above.

We run LINNA using the two publicly available posterior sampling code emcee
[15] and Zeus [2]. The convergence criteria are set based on the integrated
autocorrelation time $\tau$ [16]. Specifically, let $\tau_{i}$ be the
integrated autocorrelation time estimated using a chain of length $i$ for each
walker. A chain is considered converged if the following three criteria are
met:

  1. 1.

The autocorrelation time estimation is stable, i.e.,
$\tau_{i}/\tau_{i-100}-1<\eta$, where $\eta$ is a user-specified parameter.

  2. 2.

The number of steps per walker in units of the autocorrelation time is large,
i.e., $i/\tau_{i}>\alpha$, where $\alpha$ is user-specified criterion.

  3. 3.

The mean and error of each parameter estimated from the MCMC chain are stable.
We first select the last $\beta\times\tau_{i}$ samples for each walker to
remove the burn-in steps. We compute the estimated mean of each parameter
($p$) using the first and second half of the selected MCMC samples
($\mu_{1,p}$ and $\mu_{2,p}$ respectively). The error of each parameter is
also estimated using the standard deviation of the first and second half of
the MCMC samples ($\sigma_{1,p}$ and $\sigma_{2,p}$ respectively). The chain
is converged if
$\rm{max}_{p}(\|\mu_{1,p}-\mu_{2,p}\|/\sigma_{2,p})<\delta_{\mu}$ and
$\rm{max}_{p}(\|\sigma_{1,p}/\sigma_{2,p}-1\|)<\delta_{\sigma}$.

The above criteria are checked for every $100$ evaluation per walker. The
parameters $\eta$, $\alpha$, $\delta_{\mu}$, and $\delta_{\sigma}$ set the
convergence criteria for the posterior sampling. We remove burn-in steps by
keeping only the last $\beta\times\tau_{i}$ steps of each chain. Empirically,
we have found the following settings provide accurate posteriors using 128
walkers:

  * •

For emcee, the values of the parameters $\eta$, $\alpha$, and $\beta$ are set
to $\eta=[0.03,0.03,0.02,0.01]$, $\alpha=[5,5,10,15]$, and $\beta=[2,2,5,4]$,
where $[...]$ denotes the value for each iteration. $\delta_{\mu}$ and
$\delta_{\sigma}$ are set to $0.2$ and $0.15$ respectively for all iterations.

  * •

For Zeus, all parameters are kept the same, except for the value of $\alpha$
in the last iteration. The speed of Zeus allows us to run a much longer chain
with a modest increase of the overall run time. We therefore adopt $\alpha=50$
in the last iteration.

###  2.2 Neural network

####  2.2.1 Neural network architecture

Figure 2: The architecture of the neural network in LINNA. Green blocks
represent input of the network and orange blocks represent the output. Numbers
in paratheses show the size of input and output vectors of each block. Hidden
size is set to the minimum value of 32 times output vector size and $1000$.
Relu represents the Rectified Linear Unit, a widely used nonlinear function.
The symbol $\bigoplus$ denotes element-wise addition.

An artificial neural network (or simply neural network) is an interconnection
of nodes that mimic the neural structures of the brain. The basic components
of a neural network are nodes and connections between nodes. A node receives
multiple inputs, calculates their weighted sums, adds an optional bias term,
and passes the weighted sum through a non-linear function to produce the
output. Connections between nodes specify the weights used in the weighted sum
calculations. A neural network architecture specifies how nodes are connected,
and it can drastically affect the performance of a neural network on a
specific problem [e.g. 17]. One commonly used neural network architecture
simply arranges nodes into multiple groups (often called layers) and connects
all the nodes from one group to another. This network structure is often
referred to as a multilayer perceptron and has been widely used in the
cosmology community [e.g. 18, 19, 20, 12].

Despite their broad applicability, neural networks with multiple layers are
difficult to train and their accuracy is not necessarily better than neural
networks with fewer layers. This is known as the Degradation Problem: with
increasing numbers of layers, the accuracy of a neural network gets saturated
and can degrade [21]. This problem can be mitigated via introductions of
“shortcut” connections [21], which provide connections that skip one or more
layers. Inspired by [21], we modify the widely used multilayer perceptron
architectures, which consists of multiple fully connected layers, by adding
three shortcut connections for LINNA. The resulting neural network
architecture is shown in figure 2.

####  2.2.2 Loss function

The loss function is judiciously chosen based on the task at hand, namely
recovering accurate posteriors from the data. To that end, the loss function
is motivated by two observations.

First, the accuracy requirement of the surrogate model is set by the noise of
the data, i.e., noisy data do not require accurate models. Consequently,
instead of using a uniformly weighted loss function, we weight the differences
between the training data and the surrogate model by the inverse of the
covariance matrix of the data.

Second, for the purposes of sampling posteriors, the surrogate model needs
only to achieve high accuracy over high-likelihood regions. That is, whether a
point in the parameter space has a posterior value that is $10^{-10}$ or
$10^{-20}$ is irrelevant; we only need to know that this point is ruled out at
high confidence. For this reason, we upweight the importance of those points
in the parameter space that are a good fit to the data by dividing their
contribution to the loss function by the $\chi^{2}$ of the model. Our loss
function is therefore defined as

$Loss=\left\langle\frac{(NN-M)^{T}C^{-1}(NN-M)}{(M-d)^{T}C^{-1}(M-d)}\right\rangle,$
(2.2)

where $NN$ denotes the output of the neural network, $M$ represents the
(actual) model prediction, $C$ is the covariance matrix, $d$ is the data
vector, and $\langle\dots\rangle$ refers to the mean over training points.

####  2.2.3 Neural network optimization

The accuracy of the neural network increases as the training set expands, at
the cost of increasing the computation expense of generating training data and
training the neural network. We find that in each iteration of our algorithm,
expanding the set of points used for training and validation by $10,000$ and
$500$ respectively leads to a good compromise between accuracy and speed (see
appendix B for details).

We also find that the performance improves if the training data are
standardized before they are used to train the neural network. That is, the
input of the neural network (i.e., the cosmological and nuisance parameters)
is rescaled such that the entire training set has zero means and unit standard
deviations. The model predictions are also rescaled judiciously. Specifically,
let $\vec{d}_{\alpha}$ be the model prediction at the training point $\alpha$,
and $d_{\alpha}^{i}$ be the $i^{{\rm th}}$ element of the data vector
$\vec{d}_{\alpha}$. We rescale each training point via

$d_{\alpha}^{i}\rightarrow\frac{d_{\alpha}^{i}}{\sigma^{i}}$ (2.3)

where $(\sigma^{i})^{2}$ is the $i^{{\rm th}}$ diagonal entry of the
covariance matrix used to compute the likelihood function. We further define

$\displaystyle\mu^{i}$ $\displaystyle=$ $\displaystyle{\rm
median_{\alpha}}(d_{\alpha}^{i})$ (2.4) $\displaystyle s^{i}$ $\displaystyle=$
$\displaystyle{\rm MAD_{\alpha}}(d_{\alpha}^{i}),$ (2.5)

where MAD refers to the median absolute deviation. Notice that the median
values that define $\mu^{i}$ and $s^{i}$ are over the set of training data
(index $\alpha$), and are computed for each individual data point (index $i$).
With these quantities in hand, we further rescale the training data via

$d_{\alpha}^{i}\rightarrow\frac{d_{\alpha}^{i}-\mu^{i}}{s^{i}}.$ (2.6)

We find that rescaling by the median/MAD performs significantly better than
rescaling by the mean and standard deviation, as the latter two are much more
sensitive to a few outliers in the training data. The original data vector
$d^{i}$ can be recovered from the emulated data vector $\tilde{d}^{i}$ by
applying the inverse transformation, i.e.

$d^{i}=(\tilde{d}^{i}s^{i}+\mu^{i})\sigma^{i}.$ (2.7)

During training, the training set is randomly shuffled and split into batches,
each containing $500$ points. The loss function is estimated using each batch
of the training data and weights of the neural network are updated using the
estimated loss function. We update the weights of the neural network using
adamw [22], with an initial learning rate determined following methods
described in [23] and an initial weight decay parameter of $10^{-4}$.

The stop condition for the training of the neural network is set using the
validation data. Specifically, after the network is updated using all batches
of the training data, we evaluate the loss function on the validation data.
The learning rate is decreased by half when the minimum loss in the validation
data is not updated after $450$ steps. During the training, we use 200 steps
to estimate the local derivative of the loss at a step. If the derivative is
negative on the training set but positive on the validation set, the weight
decay parameter is doubled to prevent over-fitting. The training process is
automatically terminated when the minimum loss in the validation data is not
updated after $500$ steps.

The posterior is sampled using the trained neural network with the use of
Intel Math Kernel Library for Deep Neural Networks
222https://github.com/rsdubtso/mkl-dnn, which increases the evaluation speed
by more than a thousand times relative to not using it.

##  3 Results and Discussion

We now test the performance of LINNA in a real-world scenario. Specifically,
we compare the posterior estimated using LINNA to that estimated using emcee
to sample the posterior function implemented in CosmoLike [24] (hereafter, the
“brute force method”). Consequently, for the purpose of this comparison, the
results presented here are obtained using emcee as the LINNA sampler and
CosmoLike as the LINNA theory model. In practice, further improvements can be
achieved using faster posterior sampling algorithms such as Zeus (see figure
10 for details). The comparisons described in the main body of the paper are
mostly visual (see figures 4, 5, and 6). A more rigorous comparison of the
LINNA and “brute force method” chains can be found in Appendix A.

We perform all analyses using $128$ CPU cores and $1$ NVIDIA GeForce RTX 2080
Ti GPU on the Sherlock supercomputer333https://www.sherlock.stanford.edu/. We
consider three survey cosmology data vectors measured from first-year
observations of DES (see a summary of these data vectors in figure 1 of [7]).
Ranked by their computational cost, the three data sets are:

  1. 1.

6$\times$2pt+N [7]: a joint analysis of all data vectors in 3$\times$2pt and
4$\times$2pt+N.

  2. 2.

4$\times$2pt+N [7, 14]: a joint analysis of galaxy clustering, cluster–galaxy
cross correlations, cluster clustering, cluster lensing, and cluster
abundances.

  3. 3.

3$\times$2pt [3]: a joint analysis of galaxy clustering, galaxy–galaxy
lensing, and cosmic shear.

In these analyses, we sample 32 (28, 26) cosmological and nuisance parameters
for 6$\times$2pt+N (4$\times$2pt+N, 3$\times$2pt) with the same priors
presented in [7].

Figure 3: An illustration of the LINNA algorithm in the
$A_{\rm{s}}$–$\Omega_{\rm{m}}$ parameter subspace, when applied to the
6$\times$2pt+N data set. Top row: distributions of training data. Training
data from previous iterations are shown as grey dots while newly added
training data are shown as blue dots. Black contours show the $68\%$ and
$95\%$ contours of the targeted posteriors as evaluated using emcee to sample
the posterior function implemented in CosmoLike (the brute force method).
Bottom rows: comparison of LINNA posteriors (orange) to the target posterior
(black) at each iteration. LINNA shows strong consistency with the brute force
method after three iterations. We note that while we only show
$\Omega_{\rm{m}}$ and $A_{\rm{s}}$ in this plot, LINNA shows strong
consistency with the brute force method after three iterations in all other
cosmological and nuisance parameters as well.

###  3.1 6$\times$2pt+N analysis

figure 3 shows the $68\%$ and $95\%$ constraints on $A_{\rm{s}}$ and
$\Omega_{{\rm m}}$ using LINNA at each iteration. The blue points in the top
row correspond to the $10,000$ newly added training data points at each
iteration, while the grey points represent the cumulative training data from
previous iterations. In the first iteration, the training set is generated by
uniformly sampling the prior volume using a Latin Hypercube. For the purposes
of Latin Hypercube sampling, the boundaries of the parameters with Gaussian
priors are set to $\pm n\sigma$ around their mean values, where $n$ is a user-
specified parameter. Throughout the analysis, we adopt $n=3$ because the
Gaussian priors are only applied to nuisance parameters, whose priors are
conservatively chosen. We see that LINNA converges in just $3$ iterations,
though we advocate running an additional high-precision iteration to verify
the convergence. For this reason, we adopt four iterations as our standard
LINNA setting.

Figure 4: Constraints on important parameters when fitting to 6$\times$2pt+N
after marginalizing over additional $24$ nuisance parameters. Black contours
are obtained with the brute force method (emcee+CosmoLike), orange contours
are obtained with LINNA, and green contours are obtained using Bayesfast.
Contours show $68\%$ and $95\%$ constraints.

figure 4 compares the final constraints obtained with LINNA (orange contours)
to those derived using the brute force method (black contours). The comparison
is restricted to the most important parameters in the analyses: $\Omega_{{\rm
m}}$, $A_{\rm{s}}$, two nonlinear intrinsic alignment parameters
($A_{\rm{IA}},\eta_{\rm{IA}}$), and four richness–mass relation parameters
($\rm{ln}\lambda_{0},A_{\rm{ln}\lambda},B_{\rm{ln}\lambda},\sigma_{\rm{ln}\lambda}$).
We find excellent agreement between the two posteriors, with a shift in the
mean parameter constraints of less than $0.2\sigma$. The one-sigma errors on
individual parameters all agree to better than $12\%$. A more detailed
comparison of the means and errors between the two posteriors is shown in
figure 7. Critically, LINNA ran $\sim 10$ hours, compared to $\sim 3$ weeks
for the brute force method. A breakdown of total run time in each component of
LINNA is further shown in table 1. As noted earlier, the run time can be
further reduced by decreasing the number of training data and by replacing
emcee with Zeus in LINNA. Specifically, LINNA can reach similar performances
on all constrained cosmological and nuisance parameters with the run time down
to $5.5$ hours.

Figure 5: Similar to figure 4 but on 4$\times$2pt+N data vector.

We also compare our results to those obtained using Bayesfast [1]
444https://github.com/HerculesJack/bayesfast, another posterior inference
acceleration code currently in development (green contours in figure 4).
Similar to LINNA, Bayesfast builds surrogate models by fitting quadratic
polynomials to training data set in an iterative fashion. At each iteration,
Bayesfast samples the approximate posterior using the No U-Turn Hamiltonian
Monte Carlo algorithm [25]. After the last iteration, it importance samples
the chain by comparing the likelihood evaluated using the surrogate model and
the likelihood evaluated using the full theory model. Using the same MCMC
convergence criteria and hardware as LINNA, we find that Bayesfast converges
in $5.8$ hours, a run time comparable to that of LINNA using Zeus with $2000$
training points. However, figure 4 demonstrates that the accuracy in the
recovered posteriors of Bayesfast is slightly worse than LINNA.

###  3.2 4$\times$2pt+N analysis and 3$\times$2pt analysis

figure 5 compares the posteriors derived for the 4$\times$2pt+N data vector
using LINNA (orange contours), the brute force method (emcee+CosmoLike, black
contours) and Bayesfast (green contours). As before, we restrict the figure to
the most important parameters for the 4$\times$2pt+N data vector. All samplers
are run with exactly the same setting as for the 6$\times$2pt+N analysis. The
agreement between LINNA and the brute force method is remarkable and
demonstrates that LINNA’s hyperparameters do not need to be retuned when
running in subsets of the original data used to tune LINNA. This is not a
trivial statement: Bayesfast posteriors, which were only modestly biased for
6$\times$2pt+N data set, are now much more strongly biased than before. figure
6 compares the posteriors derived for the 3$\times$2pt data vector from LINNA
(orange contours), the brute force method (emcee+CosmoLike, black contours),
and Bayesfast (green contours). Here, LINNA and Bayesfast are both in
excellent agreement with the brute force method.

As was the case for 6$\times$2pt+N, using LINNA for these analyses leads to
significant reductions in run time relative to the brute force method
(emcee+CosmoLike). With the baseline setting, LINNA takes 11.5 (5.4) hours to
sample the 4$\times$2pt+N (3$\times$2pt) posteriors, while the brute force
method takes $\sim 4$ (1) days respectively. Again, LINNA can be further sped
up while maintaining the same accuracy on all constrained parameters by
decreasing the number of training data and by using faster samplers. With the
same fast setting as used for 6$\times$2pt+N, the run time of LINNA is $4.2$
and $3.5$ hours for 4$\times$2pt+N and 3$\times$2pt respectively.

Figure 6: Similar to figure 4 but on 3$\times$2pt data vector.

###  3.3 Environmental Impact

Supercomputers have been recognized as a significant contributor to the
$\rm{CO}_{2}$ footprint of astronomical research [26, 27, 28]. To estimate the
potential reductions of environmental impact by using LINNA, we first
calculate the $\rm{CO}_{2}$ emission savings associated with running one
posterior inference analysis using LINNA versus the current standard approach
(e.g. emcee+CosmoLike). [26] estimates that the energy consumption per CPU
core is $\sim 53W$. The maximum electricity consumption of a NVIDIA GeForce
RTX 2080 Ti GPU on the Sherlock supercomputer is $280$W. Assuming the GPU is
run with the maximum power consumption during the entire analysis, we find
that LINNA saves 12 billion Joule of electricity consumption per posterior
inference analysis for 6$\times$2pt+N data vector. Assuming an average energy
cost of US $\$0.10$ per kW-hour, as appropriate for the United States, this
corresponds to US $\$330$ per posterior inference analysis. We translate this
energy consumption to the amount of $\rm{CO}_{2}$ emission using the EPA
Carbon Offset Calculator 555https://www.epa.gov/energy/greenhouse-gas-
equivalencies-calculator and find a reduction of 2.4 tons $\rm{CO}_{2}$
emission per analysis.

Modern cosmology analyses typically require hundreds of posterior inference
analyses to validate the analysis choices and modeling pipelines. For example,
the DES collaboration performed over 900 posterior inference analyses to
analyze the 3$\times$2pt data set measured from the first three years
observations of DES. Therefore, we estimate that the publications of the
cosmological constraints derived from first-year observations of Rubin
observatory’s Legacy Survey of Space and Time (LSST Y1) will require $\sim
1,000$ posterior inference analyses. Assuming each posterior inference
analysis of LSST Y1 has similar computational costs to those of the
6$\times$2pt+N data set of DES Y1, we estimate that adopting LINNA for the
LSST Y1 analysis will save upwards of US $\$300,000$ in energy costs, while
simultaneously reducing the $\rm{CO}_{2}$ footprint of the analysis by $2,400$
tons of $\rm{CO}_{2}$. This is roughly equivalent to $2,400$ transatlantic
flights, or the yearly carbon footprint of $65$ astronomers [29]. In practice,
we expect the number of chains is likely to exceed our estimate, as we
anticipate the data will be analyzed using many different theoretical
approaches. In other words, relative to the DESY3 analyses, we have not scaled
our estimates with collaboration size or model complexity in any way.

##  4 Conclusions

Posterior inference for modern survey cosmological analyses is computationally
expensive. This is mostly due to the increasing complexity of the theory
models needed to yield unbiased cosmological constraints from unprecedentedly
precise measurements. This situation becomes more challenging when one
considers multi-probe analyses, in which one analyzes different cosmological
probes simultaneously to yield more robust and precise cosmological
constraints. The complex data vector of multi-probe analyses requires one to
sample more nuisance parameters that describe systematics of each individual
probe, thereby severely increasing the dimensionality of the posteriors. These
difficulties are especially pronounced in the case of cluster cosmology.
Modeling of cluster-related two-point correlation functions requires
marginalization over observable–mass relation and cluster-related systematics,
such as selection biases [e.g. 14], cluster mis-centering [e.g. 30], and
correlation coefficients between different cluster mass–observable connections
[e.g. 31]. These additional modeling complexities significantly increase the
computational cost of posterior inference. For instance, the run time for the
3$\times$2pt analysis of the first-year data from the Dark Energy Survey [3]
was $\sim$ 1 day. With the same hardware, the corresponding 6$\times$2pt+N
analysis [7] took $\sim 3$ weeks.

In this work, we expedite the posterior inference by constructing a surrogate
of expensive cosmological models using a deep neural network (see figure 2).
We design an iterative process to build a training data set that both spans
the entire prior volume and reasonably samples the posterior with a limited
amount of total training data (summarized in figure 1). The trained surrogate
model can be combined with standard posterior sampling tools to perform
accurate Bayesian posterior inference analyses. We test the developed
posterior inference accelerator LINNA on data vectors measured from the Dark
Energy Survey first-year observations. In particular, we consider three data
vectors: 3$\times$2pt (a joint analysis of galaxy clustering, galaxy–galaxy
lensing, and cosmic shear), 4$\times$2pt+N (a joint analysis of cluster–galaxy
cross correlations, cluster lensing, cluster clustering, and cluster
abundances), and 6$\times$2pt+N (a joint analysis of data vectors in
3$\times$2pt and 4$\times$2pt+N). The results of our tests can be summarized
as follows:

  1. 1.

In all data vectors, the posteriors obtained from LINNA agree with those from
the brute force method that uses emcee to sample the posterior implemented in
CosmoLike. The hyperparameters of LINNA are fine-tuned for 6$\times$2pt+N, but
stay the same for 4$\times$2pt+N and 3$\times$2pt.

  2. 2.

In the 6$\times$2pt+N (4$\times$2pt+N, 3$\times$2pt) analysis, the brute force
method took $\sim 3$ weeks ($\sim 4$ days, $\sim 1$ day) to sample the
posterior using 128 CPUs while LINNA took 10.1 (11.5, 5.4) hours using 128
CPUs and $1$ NVIDIA GeForce RTX 2080 Ti GPU. With the same performance on all
constrained parameters, the run time of LINNA can be further reduced to $5.5$
($4.2,3.5$) hours by reducing the number of training data and adopting Zeus as
the posterior sampler.

  3. 3.

The reduced computational demands lead to a reduction of 12 billion Joules of
electricity consumption per posterior inference analysis. We estimate that
using LINNA for the LSST Y1 cosmological analysis will save $\$300,000$ in
energy costs, while simultaneously reducing $\rm{CO}_{2}$ emission by $2,400$
tons, equivalent to the annual carbon footprint of $\sim 65$ astronomers.

While we have demonstrated the performance of LINNA using DES Y1 data set, we
have also explicitly verified that LINNA accurately reproduces the forecasted
constraints for LSST Y10 for the 3$\times$2pt and 6$\times$2pt+N data sets
(see appendix C for details). Further, although the run times measured in this
paper are based on CosmoLike and emcee, LINNA allows users to replace these
components with other modeling codes and different samplers such as Zeus [2]
or MultiNest [32] to further speed up the analysis. We make the LINNA package
publicly available at https://github.com/chto/linna, in the hope that it will
prove useful for ongoing and future astronomical and cosmological analyses.

For discussions on an independently developed likelihood acceleration tool for
LSST Year 1 3$\times$2pt analyses, we refer the reader to [33].

## Acknowledgements

CHT thanks Wei-Lin Chiang for helpful discussions on GPUs and neural networks
and David Weinberg for insightful discussions. This work received support from
the United States Department of Energy, Office of High Energy Physics under
Award Number DE-SC-0011726 and under contract number DE-AC02-76SF00515 to SLAC
National Accelerator Laboratory. ER is supported by DOE grant DE-SC0009913,
NSF grant 2009401, and a Cottrell Scholar award. HYW is supported by DOE Grant
DE-SC0021916 and NASA grant 15-WFIRST15-0008.

We acknowledge the use of GetDist [34], Pytorch [35], Matplotlib [36], and
NumPy [37] for the analyses. Some of the computing for this project was
performed on the Sherlock cluster at Stanford. We would like to thank Stanford
University, the Kavli Institute for Particle Astrophysics and Cosmology, and
the Stanford Research Computing Center for providing computational resources
and support that contributed to these research results.

## References

  * [1] H. Jia and U. Seljak, BayesFast: A Fast and Scalable Method for Cosmological Bayesian Inference, in prep. (2022).
  * [2] M. Karamanis, F. Beutler and J. A. Peacock, zeus: A python implementation of ensemble slice sampling for efficient bayesian parameter inference, arXiv preprint arXiv:2105.03468 (2021).
  * [3] T. M. C. Abbott, F. B. Abdalla, A. Alarcon, J. Aleksić, S. Allam, S. Allen et al., Dark Energy Survey year 1 results: Cosmological constraints from galaxy clustering and weak lensing, Phys. Rev. D 98 (Aug., 2018) 043526, [arXiv:1708.01530].
  * [4] E. van Uitert, B. Joachimi, S. Joudaki, A. Amon, C. Heymans, F. Köhlinger et al., KiDS+GAMA: cosmology constraints from a joint analysis of cosmic shear, galaxy-galaxy lensing, and angular clustering, MNRAS 476 (June, 2018) 4662–4689, [arXiv:1706.05004].
  * [5] S. Joudaki, C. Blake, A. Johnson, A. Amon, M. Asgari, A. Choi et al., KiDS-450 + 2dFLenS: Cosmological parameter constraints from weak gravitational lensing tomography and overlapping redshift-space galaxy clustering, MNRAS 474 (Mar., 2018) 4894–4924, [arXiv:1707.06627].
  * [6] DES Collaboration, T. M. C. Abbott, M. Aguena, A. Alarcon, S. Allam, O. Alves et al., Dark Energy Survey Year 3 Results: Cosmological Constraints from Galaxy Clustering and Weak Lensing, arXiv e-prints (May, 2021) arXiv:2105.13549, [arXiv:2105.13549].
  * [7] C. To, E. Krause, E. Rozo, H. Wu, D. Gruen, R. H. Wechsler et al., Dark Energy Survey Year 1 Results: Cosmological Constraints from Cluster Abundances, Weak Lensing, and Galaxy Correlations, Phys. Rev. Lett. 126 (Apr., 2021) 141301, [arXiv:2010.01138].
  * [8] T. Auld, M. Bridges and M. P. Hobson, COSMONET: fast cosmological parameter estimation in non-flat models using neural networks, MNRAS 387 (July, 2008) 1575–1582, [astro-ph/0703445].
  * [9] S. Agarwal, F. B. Abdalla, H. A. Feldman, O. Lahav and S. A. Thomas, PkANN - II. A non-linear matter power spectrum interpolator developed using artificial neural networks, MNRAS 439 (Apr., 2014) 2102–2121, [arXiv:1312.2101].
  * [10] M. Pellejero-Ibañez, R. E. Angulo, G. Aricó, M. Zennaro, S. Contreras and J. Stücker, Cosmological parameter estimation via iterative emulation of likelihoods, MNRAS 499 (Dec., 2020) 5257–5268, [arXiv:1912.08806].
  * [11] A. Spurio Mancini, D. Piras, J. Alsing, B. Joachimi and M. P. Hobson, COSMOPOWER: Emulating cosmological power spectra for accelerated Bayesian inference from next-generation surveys, MNRAS (Jan., 2022).
  * [12] J. DeRose, S.-F. Chen, M. White and N. Kokron, Neural Network Acceleration of Large-scale Structure Theory Calculations, arXiv e-prints (Dec., 2021) arXiv:2112.05889, [arXiv:2112.05889].
  * [13] A. Drlica-Wagner, I. Sevilla-Noarbe, E. S. Rykoff, R. A. Gruendl, B. Yanny, D. L. Tucker et al., Dark Energy Survey Year 1 Results: The Photometric Data Set for Cosmology, ApJS 235 (Apr., 2018) 33, [arXiv:1708.01531].
  * [14] C.-H. To, E. Krause, E. Rozo, H.-Y. Wu, D. Gruen, J. DeRose et al., Combination of cluster number counts and two-point correlations: validation on mock Dark Energy Survey, MNRAS 502 (Apr., 2021) 4093–4111, [arXiv:2008.10757].
  * [15] D. Foreman-Mackey, D. W. Hogg, D. Lang and J. Goodman, emcee: The MCMC Hammer, PASP 125 (Mar., 2013) 306, [arXiv:1202.3665].
  * [16] J. Goodman and J. Weare, Ensemble samplers with affine invariance, Communications in Applied Mathematics and Computational Science 5 (Jan., 2010) 65–80.
  * [17] A. Canziani, A. Paszke and E. Culurciello, An analysis of deep neural network models for practical applications, CoRR abs/1605.07678 (2016) [arXiv:1605.07678].
  * [18] T. Auld, M. Bridges, M. P. Hobson and S. F. Gull, Fast cosmological parameter estimation using neural networks, MNRAS 376 (Mar., 2007) L11–L15, [astro-ph/0608174].
  * [19] S. Agarwal, F. B. Abdalla, H. A. Feldman, O. Lahav and S. A. Thomas, PkANN - I. Non-linear matter power spectrum interpolation through artificial neural networks, MNRAS 424 (Aug., 2012) 1409–1418, [arXiv:1203.1695].
  * [20] C. Modi, Y. Feng and U. Seljak, Cosmological reconstruction from galaxy light: neural network based light-matter connection, J. Cosmology Astropart. Phys 2018 (Oct., 2018) 028, [arXiv:1805.02247].
  * [21] K. He, X. Zhang, S. Ren and J. Sun, Deep Residual Learning for Image Recognition, arXiv e-prints (Dec., 2015) arXiv:1512.03385, [arXiv:1512.03385].
  * [22] Z. Xie, I. Sato and M. Sugiyama, Understanding and Scheduling Weight Decay, arXiv e-prints (Nov., 2020) arXiv:2011.11152, [arXiv:2011.11152].
  * [23] L. N. Smith, Cyclical learning rates for training neural networks, in 2017 IEEE winter conference on applications of computer vision (WACV), pp. 464–472, IEEE, 2017.
  * [24] E. Krause and T. Eifler, cosmolike - cosmological likelihood analyses for photometric galaxy surveys, MNRAS 470 (Sept., 2017) 2100–2112, [arXiv:1601.05779].
  * [25] M. D. Homan and A. Gelman, The no-u-turn sampler: Adaptively setting path lengths in hamiltonian monte carlo, J. Mach. Learn. Res. 15 (jan, 2014) 1593–1623.
  * [26] A. R. H. Stevens, S. Bellstedt, P. J. Elahi and M. T. Murphy, The imperative to reduce carbon emissions in astronomy, Nature Astronomy 4 (Sept., 2020) 843–851, [arXiv:1912.05834].
  * [27] L. Burtscher, H. Dalgleish, D. Barret, T. Beuchert, A. Borkar, F. Cantalloube et al., Forging a sustainable future for astronomy, Nature Astronomy 5 (Sept., 2021) 857–860.
  * [28] F. van der Tak, L. Burtscher, S. Portegies Zwart, B. Tabone, G. Nelemans, S. Bloemen et al., The carbon footprint of astronomy research in the Netherlands, Nature Astronomy 5 (Dec., 2021) 1195–1198.
  * [29] J. Knödlseder, S. Brau-Nogué, M. Coriat, P. Garnier, A. Hughes, P. Martin et al., Estimate of the carbon footprint of astronomical research infrastructures, arXiv e-prints (Jan., 2022) arXiv:2201.08748, [arXiv:2201.08748].
  * [30] T. McClintock, T. N. Varga, D. Gruen, E. Rozo, E. S. Rykoff, T. Shin et al., Dark energy survey year 1 results: weak lensing mass calibration of redmapper galaxy clusters, MNRAS 482 (Jan., 2019) 1352–1378, [arXiv:1805.00039].
  * [31] S. Bocquet, J. P. Dietrich, T. Schrabback, L. E. Bleem, M. Klein, S. W. Allen et al., Cluster Cosmology Constraints from the 2500 deg2 SPT-SZ Survey: Inclusion of Weak Gravitational Lensing Data from Magellan and the Hubble Space Telescope, ApJ 878 (June, 2019) 55, [arXiv:1812.01679].
  * [32] F. Feroz, M. P. Hobson and M. Bridges, MULTINEST: an efficient and robust Bayesian inference tool for cosmology and particle physics, MNRAS 398 (Oct., 2009) 1601–1614, [arXiv:0809.3437].
  * [33] S. S. Boruah, T. Eifler and V. Miranda, Accelerating cosmological inference with Gaussian processes and neural networks, in prep. (2022).
  * [34] A. Lewis, GetDist: a Python package for analysing Monte Carlo samples, arXiv:1910.13970.
  * [35] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan et al., Pytorch: An imperative style, high-performance deep learning library, Advances in neural information processing systems 32 (2019) 8026–8037.
  * [36] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing in Science & Engineering 9 (2007), no. 3 90–95.
  * [37] S. van der Walt, S. C. Colbert and G. Varoquaux, The numpy array: A structure for efficient numerical computation, Computing in Science Engineering 13 (March, 2011) 22–30.
  * [38] The LSST Dark Energy Science Collaboration, R. Mandelbaum, T. Eifler, R. Hložek, T. Collett, E. Gawiser et al., The LSST Dark Energy Science Collaboration (DESC) Science Requirements Document, arXiv e-prints (Sept., 2018) arXiv:1809.01669, [arXiv:1809.01669].

##  Appendix A Multidimensional comparison of two MCMC chains

We wish to compare LINNA chains to those obtained using the brute force
method. To do so, we compare the means and errors for each of the parameters
in the model, as shown in figure 7. In the top row of figure 7, we show the
difference in the mean between LINNA and the brute force method for each model
parameter, in units of the standard deviation of the same (as evaluated using
the brute force method). We find that the difference between the two means is
always less than $0.2\sigma$, with the median difference being $0.2\sigma$.
The figure also compares the emcee chains from the brute force method to those
generated using MultiNest. Notice the LINNA posteriors are in better agreement
with the brute force method than those estimated using MultiNest.

We further compare the uncertainty estimations from LINNA and those from the
brute force method in the bottom row of figure 7. To make this comparison, we
compute the parameter covariance matrix using chains from LINNA and chains
from the brute force method. We rotate the covariance matrix estimated using
the LINNA chain to the eigen-space of the covariance matrix estimated using
the chain from the brute force method (reference covariance matrix hereafter).
If the two covariance matrices are consistent, the rotated covariance matrix
will be mostly diagonal with the diagonal value being identical to the eigen-
value of the reference covariance matrix. In the bottom two panels of figure
7, we show that the rotated covariance matrix is mostly diagonal and the
diagonal value is nearly identical to the eigen-value of the reference
covariance matrix. Similar comparisons for 4$\times$2pt+N and 3$\times$2pt
data vectors are shown in figure 8 and 9.

Figure 7: A comparison of multidimensional posteriors obtained with LINNA and
with emcee from the 6$\times$2pt+N data vector. Top: a comparison of parameter
means obtained with different samplers to $1\sigma$ uncertainties using emcee
chains. Bottom left: a comparison of eigen-value of parameter covariance
matrix estimated using emcee chains (black) and diagonal value of parameter
covariance matrix estimated using LINNA ( Bayesfast) chains after they are
rotated to the eigen-space of former (orange (green)). Bottom right: the
correlation matrix of parameters estimated using LINNA chains after they are
rotated to the eigen-space of parameter covariance matrix estimated using
emcee chains. Figure 8: Similar to figure 7 but on the 4$\times$2pt+N data
vector. Figure 9: Similar to figure 7 but on the 3$\times$2pt data vector.

##  Appendix B Optimizations of the number of training data and run-time
breakdowns

In figure 10, we compare the performance of LINNA with different numbers of
training data points added in each iteration. We find that as long as the
number of training data points added in each iteration is more than 2000,
LINNA leads to less than $0.2\sigma$ shifts in all parameters that are
meaningfully constrained by the data. With 10000 training data points added in
each iteration, all parameters are shifted by less than $0.2\sigma$ when
compared to the brute force method. We thus adopt adding $10000$ training data
points in each iteration as the default option for LINNA, but note that $2000$
training data points could lead to a reasonable performance in all constrained
parameters. In table 1, we show the breakdown of the total run time of LINNA
with $10000$ and $2000$ training data points added in each iteration
respectively. The total run time reduces from 10.1 hours to $6.89$ hours when
reducing the number of training data points from $10000$ to $2000$. We further
note that in the case of $2000$ training data points, the LINNA run time is
dominated by the MCMC sampling. One could further speed up LINNA using faster
samplers such as Zeus. We find that using Zeus reduces the run time for the
$2000$ training data point case to $5.52$ hours while maintaining the same
performance.

Figure 10: A comparison of LINNA performance with different numbers of training data points. The $y$ axis shows the shifts in parameter mean estimations using LINNA and using the brute force method in the unit of $1\sigma$ uncertainties. Different color lines correspond to LINNA having $1000-10000$ additional training data points per iteration. Table 1: Breakdowns of LINNA run time in each component using the default option. Model shows the amount of time spent in generating the training data using CosmoLike. Train shows the time spent in training the neural network. MCMC shows the time spent in sampling the posterior using emcee. Numbers in each panel are in units of hours. Numbers in the parentheses show run time for LINNA with 2000 training data points in each iteration.  Analysis: 6$\times$2pt+N |  |  |
---|---|---|---
| iteration 1 | iteration 2 | iteration 3 | iteration 4 | Total
Model | 1.12 (0.46) | 0.67 (0.18) | 0.74 (0.19) | 0.78 (0.19) | 3.31 (1.02)
Train | 0.51 (0.56) | 0.61 (0.61) | 1.23 (0.56) | 1.25 (0.69) | 3.60 (2.42)
MCMC | 0.29 (0.30) | 0.28 (0.25) | 0.51 (0.58) | 2.08 (2.33) | 3.16 (3.46)
Total | 1.93 (1.32) | 1.55 (1.05) | 2.48 (1.33) | 4.11 (3.21) | 10.07 (6.89)
Analysis: 4$\times$2pt+N |  |  |
| iteration 1 | iteration 2 | iteration 3 | iteration 4 | Total
Model | 1.05 (0.74) | 0.69 (0.25) | 0.71 (0.26) | 0.72 (0.26) | 3.18 (1.51)
Train | 0.73 (0.22) | 0.89 (0.38) | 1.09 (0.52) | 1.57 (0.52) | 4.28 (1.64)
MCMC | 0.26 (0.44) | 0.30 (0.48) | 0.61 (1.57) | 2.82 (6.76) | 3.99 (9.25)
Total | 2.05 (1.40) | 1.88 (1.11) | 2.42 (2.35) | 5.11 (7.54) | 11.45 (12.40)
Analysis: 3$\times$2pt |  |  |
| iteration 1 | iteration 2 | iteration 3 | iteration 4 | Total
Model | 0.55 (0.35) | 0.14 (0.04) | 0.14 (0.08) | 0.14 (0.04) | 0.99 (0.51)
Train | 0.48 (0.26) | 0.63 (0.62) | 0.80 (0.33) | 1.25 (0.46) | 3.16 (1.67)
MCMC | 0.20 (0.21) | 0.21 (0.20) | 0.29 (0.28) | 0.56 (0.57) | 1.25 (1.26)
Total | 1.23 (0.82) | 0.98 (0.86) | 1.24 (0.70) | 1.95 (1.06) | 5.39 (3.44)

##  Appendix C LSST Year 10 forecast

We evaluate the expected performance of LINNA on ten-year LSST (LSST Y10) data
using simulated data vectors. We construct the data vector using CosmoLike,
using parameters described in [38]. For details, we refer the reader to To,
Krause et al. in prep. Here, we provide a short description. The lens galaxy
samples are split into 10 equally spaced photometric redshift bins from
$z=0.2$ to 1.2. The number density is $48$ $\rm{arcmin}^{-2}$. For the source
galaxy samples, we assume the number density to be $27$ $\rm{arcmin}^{-2}$ and
split them into five redshift bins from $z=0.2$ to $4$. The cluster samples
are generated using the richness–mass relation obtained in [7], and are split
into four equally spaced redshift bins from $z=0.2$ to $1$. In the analysis,
we assume a wCDM cosmological model with $\rm{w}_{0}$ and $\rm{w}_{a}$
describing the dark energy equation of state and its redshift dependence.

In figures 11 and 12, we compare the result of LINNA with the result obtained
from the brute force method. We find that the level of agreement is consistent
with what we found in analyzing DES Y1 data. These results indicate that LINNA
is expected to be applicable on the analysis of data with LSST Y10 precision.

Figure 11: Performance of LINNA on the 3$\times$2pt analysis using the
simulated LSST Y10 data vector. Details of this plot are similar to figure 7.
Figure 12: Performance of LINNA on the 6$\times$2pt+N analysis using the
simulated LSST Y10 data vector. Details of this plot are similar to figure 7.