# Higher-Order Moment-Based Anomaly Detection

Venkatraman Renganathan, Navid Hashemi, Justin Ruths, and Tyler H. Summers
This work was partially supported by the Army Research Office and was
accomplished under Grant Number: W911NF-17-1-0058 and the Air Force Office of
Scientific Research under award number FA2386-19-1-4073.The authors are with
the Department of Mechanical Engineering, The University of Texas at Dallas,
800 W. Campbell Rd, Richardson, TX, USA. Email: (vrengana, navid.hashemi,
jruths<EMAIL_ADDRESS>
###### Abstract

The identification of anomalies is a critical component of operating complex,
large-scale and geographically distributed cyber-physical systems. While
designing anomaly detectors, it is common to assume Gaussian noise models to
maintain tractability; however, this assumption can lead to the actual false
alarm rate being significantly higher than expected. Here we design a
distributionally robust threshold of detection using finite and fixed higher-
order moments of the detection measure data such that it guarantees the actual
false alarm rate to be upper bounded by the desired one. Further, we bound the
states reachable through the action of a stealthy attack and identify the
trade-off between this impact of attacks that cannot be detected and the
worst-case false alarm rate. Through numerical experiments, we illustrate how
knowledge of higher-order moments results in a tightened threshold, thereby
restricting an attacker’s potential impact.

##  I INTRODUCTION

From critical infrastructures and industrial process control to autonomous
driving and various biomedical applications, dynamical control systems are
increasingly able to be instrumented with new sensing and actuation
capabilities. These cyber-physical systems (CPS) comprise growing webs of
interconnected feedback loops and must operate efficiently and resiliently in
dynamic and uncertain environments. As these systems become large, devising
both model-based [1, 2] and data-driven [3, 4, 5] methods for detecting
anomalies (such as component failures or malicious attacks) are critical for
their robust and efficient operation. Such critically important cyber-physical
networks have become an attractive target to attackers. These systems are
large and complex and are often not monitored well enough, enabling attackers
to manipulate the system without being detected and cause damage [6, 7].

To simplify the analysis and design, often such complex cyber-networks are
modeled as a discrete-time linear time invariant system with Gaussian noises.
However, this can lead to a significant miscalculation of probabilities and
risk if the underlying processes behave differently, for example due to
various nonlinearities or malicious attacks. In the context of attacks, it is
possible for an attacker to modify the sensor outputs and effectively generate
aggressive and strategic noise profiles to sabotage the operation of the
system. With stochastic optimization techniques, particularly using the
emerging area of distributionally robust optimization (DRO) approaches [8],
these limitations can be recognized and addressed. DRO enables modelers to
explicitly incorporate inherent ambiguity in probability distributions into
optimization problems. DRO approaches can be categorized based on the form of
the ambiguity set. There are several different parameterizations, including
those based on moments, support, directional derivatives [9], and Wasserstein
balls [10]. In practice, we have access to only a finite amount of historical
data. However, it is possible to use finite historical data and guarantee
resiliency in such critical cyber-physical networks. Here, we propose to use
moment-based DRO methods to improve modeling and reduce false alarm rates in
cyber-physical networks.

In the context of attacks, the detector tuning has a direct implication on the
effect an attacker can have while still remaining stealthy [11, 2]. A model-
based approach to attack detection uses a detector that raises alarms when
there is a large enough discrepancy between the actual and predicted
measurements, a statistic termed the residual. The detector’s sensitivity can
be increased by decreasing the threshold of detection, but there is an
inherent trade-off between sensitivity and the rate at which false alarms are
generated. Keeping false alarms to a manageable level requires adjusting
sensitivity and the tuning of the detector threshold is typically informed by
the distribution of the residual. Authors in [12] used Gaussian Mixture Model
to approximate the arbitrary noise distributions and obtained a detector
threshold corresponding to a desired false alarm rate.

However, when noise distributions are only known to an ambiguity set,
traditional tools and approximations no longer suffice to select the threshold
and so we turn to a distributionally robust approach. Interest in a DRO-
informed perspective on detector tuning is supported by recent work on using a
Wasserstein metric [5]. Our moment-based ambiguity set formulation includes
all distributions with fixed moments up to some order. The problem of
designing anomaly detector thresholds subject to moment constraints of system
uncertainties can be addressed using Generalized Moment Problems described in
[13, 14]. Further, the conditions for a truncated (finite) moment sequence to
represent a probability measure were studied in [15, 16]. Authors in [17]
proposed semidefinite programs to compute a probability bound for a random
variable lying in a set with known moments up to some order. These techniques
can be utilized to design detector thresholds for residual distributions
consistent with finite fixed moments up to some order.

Contributions: This paper is a significant extension of our previous work [3]
where we used a moment-based ambiguity set formulation with fixed first two
moments (Proposition 2) to obtain a detector threshold via generalized
Chebyshev inequality. The main contributions of the present paper are:

  * •

We propose an approach to construct moment-based ambiguity sets with fixed
moments up to $k\textsuperscript{th}$ order for the anomaly detection measure
and design an anomaly detector threshold for CPSs that exhibit non-Gaussian
uncertainties. The approach can utilize either residual moments obtained
through a dynamic model or estimated directly from residual data.

  * •

We use a semidefinite program (SDP) defined using higher-order moments of the
detection measure to find a sharper probability bound for classifying the
residual with an improved detector threshold (Theorem 4).

  * •

We illustrate empirically that inclusion of higher-order moments results in
tightened threshold (Lemma 3), thereby restricting an attacker’s potential
impact, and also prove that the volume of the attack-reachable set shrinks
with a tightened threshold (Corollary 5).

While anomaly detection is widely studied in CPS literature, our
distributionally robust approach using higher-order moments of the detection
measure data marks the novel contribution of this paper. The rest of the paper
is organized as follows. Section II formulates the problem using moment-based
ambiguity sets. In Section III, the design of anomaly detector threshold is
discussed. Section IV describes the procedure to find the boundary of the
reachable sets obtained using distributionally robust tuned detector. Section
V presents the numerical results with inferences. Finally, Section VI
concludes and summarizes the future research directions.

##  II Problem Formulation Using Higher-Order Moment-Based Ambiguity Sets

In this section, we propose a framework for designing an anomaly detector
threshold for cyberphysical systems. The approach can utilize either model-
based propagation of residual moments or data-driven estimation of detection
measure moments directly from data.

###  II-A Model-based Problem Formulation

Here, we model an uncertain cyber-physical system as a stochastic discrete-
time linear system

$\displaystyle x_{t+1}$ $\displaystyle=Ax_{t}+Bu_{t}+w_{t},\quad
t\in\mathbb{N}$ (1) $\displaystyle y_{t}$ $\displaystyle=Cx_{t}+v_{t},$ (2)

where $x_{t}\in\mathbb{R}^{n},u_{t}\in\mathbb{R}^{m}$ are the system state and
input respectively at time $t$. The matrices $A$ and $B$ denote the system
matrix and control input matrix, respectively. The output
$y_{t}\in\mathbb{R}^{p}$ aggregates a linear combination of the states with
the observation matrix $C\in\mathbb{R}^{p\times n}$. We assume that the pair
$(A,C)$ is detectable and $(A,B)$ is stabilizable. The process noise
$w_{t}\in\mathbb{R}^{n}$ and the sensor noise $v_{t}\in\mathbb{R}^{p}$ are
modeled as zero-mean random vectors independent and identically distributed
across time with covariance matrix $\Sigma_{w},\Sigma_{v}$ respectively. Let
$\kappa=(k_{1},\dots,k_{n})^{\top}$ with $k_{j}\in\mathbb{Z}_{+}$ non-negative
integers and $J_{k}=\\{\kappa|\sum^{n}_{i=1}k_{i}\leq k\\}$. Further, assume
that the feasible first $k$ moment sequence of the distributions of the
$w_{t},v_{t}$ namely $M^{\kappa}_{w},M^{\kappa}_{v},\kappa\in J_{k}$
respectively are known. The distributions $P_{w}$ of $w_{t}$ and $P_{v}$ of
$v_{t}$ are unknown (and not necessarily Gaussian and possibly heavy-
tailed111For the purposes of this paper, we consider heavy-tailed
distributions as those whose moments above a certain order may be infinite, in
which case their tails are heavier than a Gaussian. We assume the moments up
to order $k$ are finite.) and will be assumed to belong to the $k$-moments-
based ambiguity sets of distributions $\mathcal{P}^{w}_{k}$ and
$\mathcal{P}^{v}_{k}$ respectively defined as follows

$\displaystyle\mathcal{P}^{w}_{k}=\\{P_{w}\mid\mathbb{E}[w^{k_{1}}_{t_{1}}w^{k_{2}}_{t_{2}}\dots
w^{k_{n}}_{t_{n}}]=M^{\kappa}_{w},\ \kappa\in J_{k}\\},$ (3)
$\displaystyle\mathcal{P}^{v}_{k}=\\{P_{v}\mid\mathbb{E}[v^{k_{1}}_{t_{1}}v^{k_{2}}_{t_{2}}\dots
v^{k_{p}}_{t_{p}}]=M^{\kappa}_{v},\ \kappa\in J_{k}\\}.$ (4)

When the actual measurement $y_{t}$ is corrupted by an additive attack,
$\delta_{t}\in\mathbb{R}^{p}$, the true output of the system fed to the
controller becomes

$\displaystyle\bar{y}_{t}=y_{t}+\delta_{t}=Cx_{t}+v_{t}+\delta_{t}.$ (5)

We utilize a steady-state Kalman filter to construct a state estimate
$\hat{x}_{t}$ to minimize the squared norm of the estimation error
$e_{t}=x_{t}-\hat{x}_{t}$ where,

$\displaystyle\hat{x}_{t+1}$
$\displaystyle=A\hat{x}_{t}+Bu_{t}+L(\bar{y}_{t}-C\hat{x}_{t}),$ (6)

and the estimation error evolves as

$\displaystyle e_{t+1}$ $\displaystyle=(A-LC)e_{t}+w_{t}-Lv_{t}-L\delta_{t}.$
(7)

In the absence of attacks (that is, $\delta_{t}=0$) and when the covariance
matrices of the noises are fixed and known, the Kalman gain
$L=PC^{\top}(CPC^{\top}+\Sigma_{v})^{-1}$ minimizes the steady state
covariance matrix

$\displaystyle P$
$\displaystyle:=\lim_{t\rightarrow\infty}P_{t}:=\mathbf{E}[e_{t}e^{\top}_{t}].$
(8)

Since $(A,C)$ is assumed to be detectable, the existence of $P$ is guaranteed
and it can be found through the solution of an algebraic Ricatti equation. We
define a residual sequence $r_{t}$ as the difference between the actual
received output $\bar{y}_{t}$ and the predicted output $C\hat{x}_{t}$ as,

$\displaystyle r_{t}$
$\displaystyle=\bar{y}_{t}-C\hat{x}_{t}=Ce_{t}+v_{t}+\delta_{t},$ (9)

and in the attack free setting, $r_{t}$ falls according to a zero mean
distribution with covariance

$\displaystyle\Sigma_{r}=\mathbf{E}[r_{t}r^{\top}_{t}]=CPC^{\top}+\Sigma_{v}.$
(10)

Since (9) is linear, it is possible to obtain the fixed and first $k$ moments
of the random variable $r_{t}$ by propagating the corresponding moments of the
primitive random variables $w_{t},v_{t}$. Thus, the distribution $P_{r}$ of
$r_{t}$ (not necessarily Gaussian) belongs to an ambiguity set
$\mathcal{P}^{r}_{k}$ given by

$\displaystyle\mathcal{P}^{r}_{k}=\\{P_{r}\mid\mathbb{E}[r^{k_{1}}_{t_{1}}r^{k_{2}}_{t_{2}}\dots
r^{k_{p}}_{t_{p}}]=M^{\kappa}_{r},\ \kappa\in J_{k}\\}.$ (11)

We define detection measure $q_{t}$ as a quadratic function of $r_{t}$

$\displaystyle q_{t}=r^{\top}_{t}\Sigma^{-1}_{r}r_{t},$ (12)

which will be compared to a threshold for anomaly detection.

###  II-B Data-Driven Moment Estimation from Residual Data

An alternative to obtaining residual moments by propagating the moments of
primitive random variables through the system model is to instead collect
residual data $r_{t}$ (from attack-free operation) and estimate residual
moments or the moments of $q_{t}$ directly from the data. Such a data-driven
approach allows our proposed tuning approaches to be used in much broader
settings where it is difficult to propagate moments through a model (or even
to obtain a model), but where residual data is easily generated from sensors
and a state estimator. Higher-order moments require increasingly more data to
obtain accurate estimates. Determining the required amount of data is possible
using finite-sample measure concentration results given as in [18], but we
leave such an analysis for future work. The proposed methods for setting
thresholds for anomaly detection can be used together with ambiguity sets
built upon data-driven residual estimates, although the false alarm rates will
also be affected by sampling errors. Moment estimation uncertainty could be
accommodated using the same generalized moment problem computations we propose
here, but with assumed _bounds_ on moment estimates rather than having them
fixed to exact known values. Such problems can also be reformulated in a
computationally tractable manner, as described in, e.g., Section 3 of [19],
which incorporate uncertainty sets for estimated moments. These formulations
would simply impose further constraints on our primal problem defined in
subsection III.B of our revised manuscript. It is possible and would be
interesting to explore how to use statistical confidence intervals for data-
driven moment estimates to inform bounds used in the SDP and obtain end-to-end
guarantees on false alarm rates associated with certain thresholds. This will
be pursued in future work.

In either setting, the feasible first $k$ moment sequence of the distribution
of $q_{t}$ denoted by $M^{k}_{q}$ is assumed to be known either from its
primitive variables through the model-based approach or estimated from data.
Then, the $k$-moments-based ambiguity set of the scalar random variable
$q_{t}$ is defined as

$\displaystyle\mathcal{P}^{q}_{k}=\\{P_{q}\mid\mathbb{E}[q^{k}_{t}]=M^{k}_{q}\\}.$
(13)

######  Proposition 1.

[17] Consider a univariate random variable $X$ defined on
$\Omega=\mathbb{R}_{+}$, endowed with its Borel sigma algebra of events. A
sequence $\bar{\sigma}=(M_{1},M_{2},\dots,M_{k})^{\top}$ is a feasible
$(1,k,\Omega)-$moment vector of the random variable $X$, if and only if the
matrices $R_{k}\succeq 0$ and $R_{k-1}\succeq 0$, where for any integer $l\geq
0$, the matrices are defined as

$R_{2l}=\begin{bmatrix}1&M_{1}&\dots&M_{l}\\\ M_{1}&M_{2}&\dots&M_{l+1}\\\
\vdots&\vdots&\ddots&\vdots\\\
M_{l}&M_{l+1}&\dots&M_{2l}\end{bmatrix},R_{2l+1}=\begin{bmatrix}M_{1}&\dots&M_{l+1}\\\
M_{2}&\dots&M_{l+2}\\\ \vdots&\ddots&\vdots\\\
M_{l+1}&\dots&M_{2l+1}\end{bmatrix}.$ (14)

##  III Design of Anomaly Detector Thresholds

Given a threshold 222The first and second subscripts in the threshold
separated by a comma denote the random variable and number of moments
respectively.$\alpha_{q,k}>0$ and the distance measure $q_{t}$, alarm time(s)
$t^{\star}$ are produced according to the following rules,

$\displaystyle\begin{cases}q_{t}\leq\alpha_{q,k},&\text{no alarm}\\\
q_{t}>\alpha_{q,k},&\text{alarm: }t^{\star}=t.\end{cases}$ (15)

Even in the absence of attacks, the detector is expected to generate false
alarms due to the infinite support of $v_{t}$, because some values drawn from
$P_{q}$ will exceed the threshold $\alpha_{q,k}$. If $P_{q}$ is known, then it
is possible to extract an optimum threshold value $\alpha_{q}^{*}$ from the
corresponding cumulative distribution function $F_{q}$ for a desired false
alarm rate, $\mathcal{A}$. For example, if $r_{t}$ is Gaussian, $q_{t}$ would
be a chi-squared random variable, and the optimum threshold
$\alpha^{\star}_{q}$ corresponding to the desired false alarm rate
$\mathcal{A}=\mathcal{A}^{\star}$ is then

$\alpha^{\star}_{q}=\alpha_{\chi^{2}}:=2\mathbb{P}^{-1}\left(1-\mathcal{A}^{\star},\frac{p}{2}\right),$
(16)

where $\mathbb{P}^{-1}(\cdot,\cdot)$ denotes the inverse regularized lower
incomplete gamma function [1]. When the complete distribution is not
available, tuning methods using (16) may design thresholds that generate
actual false alarm rates significantly higher than what is desired. With the
distributionally robust approach, we aim to achieve a false alarm rate less
than $\mathcal{A}$, and the detector threshold $\alpha^{\star}_{r,k}$ is
selected such that

$\displaystyle\sup_{P_{r}\in\mathcal{P}^{r}_{k}}P_{r}(r^{\top}_{t}\Sigma^{-1}_{r}r_{t}\leq\alpha^{\star}_{r,k})$
$\displaystyle=1-\mathcal{A}.$ (17)

When the first two moments $(k=2)$ of $r_{t}$ are known, the following
proposition using the generalized Chebyshev inequality explained in [20, 21]
can be used to obtain the worst case detector threshold $\alpha^{\star}_{r,2}$
satisfying (17).

######  Proposition 2.

[3] Given a desired false alarm rate $\mathcal{A}$ and $r_{t}\sim
P_{r}\in\mathcal{P}^{r}_{k}$ with $k=2$, the optimal distributionally robust
threshold $\alpha^{\star}_{r,2}$ satisfying (17) is

$\displaystyle\alpha^{\star}_{r,2}=\frac{p}{\mathcal{A}}.$ (18)

###  III-A Improved Detector Threshold With Higher-Order Moments

It is possible to obtain a sharpened detector threshold than
$\alpha^{\star}_{r,2}$, if higher-order moments are taken into account. For
example, the skewness and kurtosis parameters convey asymmetry and heaviness
of tails of the distribution, respectively. We can leverage such information
about the true but unknown distribution revealed by the higher-order moments
to tighten the required probability bound and thereby obtain an improved
detector threshold. It is possible to use $r_{t}$ with higher order moments
(first $k$ moments) to obtain a sharpened detector threshold
$\alpha^{\star}_{r,k}$. However, it was shown in [17] that for
$r_{t}\in\mathbb{R}^{p}$ with support $\Omega=\mathbb{R}^{p}$ and $k\geq 4$,
it is NP-hard to find tight bounds for the corresponding moment bound problem
with rational problem data. On the other hand, they provide a semidefinite
optimization problem in $k+1$ dimension for the case of a univariate random
variable with $k$ moments. Hence, rather looking for higher-order moments of
$r_{t}$, instead we look for higher-order moments of scalar random variable
$q_{t}$. Subsequently, we use a bisection algorithm to obtain a sharpened
detection threshold $\alpha^{\star}_{q,k}$ for a given $\mathcal{A}$.

###  III-B Estimating Probability Using $(1,k,\Omega)-$Moment Bound

Figure 1: A moment based polynomial bounding the indicator function
$\mathbf{1}_{\mathcal{S}_{k}}$ representing the set
$\mathcal{S}_{k}=\mathbb{R}_{>\alpha_{q,k}}$ is shown here.

Given the first $k$ moments of random variable $q_{t}$ with support
$\Omega=\mathbb{R}_{\geq 0}$ and the set
$\mathcal{S}_{k}=\mathbb{R}_{>\alpha_{q,k}}$ representing an alarm event as
shown in Figure 1, the infinite dimensional $(1,k,\Omega)-$moment bound primal
problem is given by

$\displaystyle\sup_{P_{q}\in\mathcal{P}^{q}_{k}}$ $\displaystyle
P_{q}(q_{t}\in\mathcal{S}_{k})$ (19) s.t.
$\displaystyle\int_{\Omega}q^{r}d\mu=M^{r}_{q},\quad r=0,1,\dots,k.$

Since (19) is infinite dimensional, it is difficult to solve efficiently.
However, we can use the linear programming duality theory to associate a dual
variable $y_{r},r=0,1,\dots,k$ with each equality constraint of (19) to get
the corresponding dual problem

min $\displaystyle\sum^{k}_{r=0}y_{r}M^{r}_{q}$ (20) s.t. $\displaystyle
p(q)=\sum^{k}_{r=0}y_{r}q^{r}\geq 1,\quad\forall q\in\mathcal{S}_{k},$
$\displaystyle p(q)=\sum^{k}_{r=0}y_{r}q^{r}\geq 0,\quad\forall q\in\Omega.$

The probability obtained as the solution to (20) is an upper bound to the
probability associated with (19). In general, if the moment vector
$\bar{\sigma}_{k}=(M^{0}_{q},M^{1}_{q},\dots M^{k}_{q})$ is an interior point
of the set $\mathcal{M}_{k}$ of all feasible moment vectors, then strong
duality exists between (19) and (20) enabling us to obtain a _tight_ bound on
$P_{q}(q_{t}\in\mathcal{S}_{k})$. To achieve a desired false alarm rate
$\mathcal{A}$, we can tune the threshold $\alpha_{q,k}$ defining the set
$\mathcal{S}_{k}$ such that solution to (20) is $\mathcal{A}$. The following
lemma establishes the general trend observed between the values of the tuned
threshold $\alpha_{q,k}$ for increasing values of $k$.

######  Lemma 3.

Let $\mathcal{A}$ be the desired false alarm rate and $q_{t}\sim
P_{q}\in\mathcal{P}^{q}_{k}$ as in (13). Then, $\forall j,k\in\mathbb{Z}_{\geq
1}$ with $j<k$, we have

$\displaystyle\alpha^{\star}_{q,k}\leq\alpha^{\star}_{q,j}$ (21)

when the associated moments agree up to order $j$.

###### Proof.

When the associated moments up to order $j<k$ agree, clearly we have
$\mathcal{P}^{q}_{k}\subseteq\mathcal{P}^{q}_{j}$, since additional moment
constraints restrict the set of distributions. For a fixed threshold
$\alpha^{\star}_{q,j}$ tuned for the false alarm rate $\mathcal{A}$, it
follows that

$\displaystyle\sup_{P_{q}\in\mathcal{P}^{q}_{k}}P_{q}(q_{t}>\alpha^{\star}_{q,j})\leq\sup_{P_{q}\in\mathcal{P}^{q}_{j}}P_{q}(q_{t}>\alpha^{\star}_{q,j})=\mathcal{A}.$
(22)

Thus, to achieve a desired false alarm rate $\mathcal{A}$ under the constraint
that $P_{q}\in\mathcal{P}^{q}_{k}$, the threshold $\alpha^{\star}_{q,k}$ must
be non-increasing in $k$, satisfying
$\alpha^{\star}_{q,k}\leq\alpha^{\star}_{q,j}$ for $j<k$. ∎

######  Theorem 4.

Let $\epsilon>0$, $k\in\mathbb{Z}_{\geq 2}$ and consider $q_{t}$ with its
first $k$ moments of distribution $(M^{1}_{q},M^{2}_{q},\dots,M^{k}_{q})$ (we
let $M^{0}_{q}=1$) defined on $\mathbb{R}_{+}$ being known (or estimated from
the given data) and the associated anomaly detector threshold $\alpha_{q,k}$
introduced in (15), which is intended to achieve a desired false alarm rate
$\mathcal{A}$. Suppose that $\alpha_{q,k}$ is obtained by solving the
following bisection algorithm (knowing that the desired
$\alpha^{\star}_{q,k}\in[\alpha_{l},\alpha_{u}]$ and
$\alpha_{u}=\alpha^{\star}_{q,k-1}$ by Lemma 3)

while  _$\alpha_{u}-\alpha_{l} >\epsilon$_ do

$\alpha^{\star}_{q,k}\leftarrow(\alpha_{u}+\alpha_{l})/2$

$p_{min}\leftarrow$ optimal value of (23) with
$\alpha_{q,k}=\alpha^{\star}_{q,k}$

if  _$p_{min} >\mathcal{A}$_ then

$\alpha_{l}\leftarrow\alpha^{\star}_{q,k}$

else

$\alpha_{u}\leftarrow\alpha^{\star}_{q,k}$

Algorithm 1 Bisection Subroutine

where the SDP that gives the tight upper bound on
$P_{q}(q_{t}\geq\alpha_{q,k})$ in the third line is

min $\displaystyle\sum^{k}_{r=0}y_{r}M^{r}_{q}$ (23) s.t.
$\displaystyle\sum_{i,j:i+j=2l-1}x_{ij}=0,$ $\displaystyle l=1,\dots,k,$
$\displaystyle(y_{0}-1)+\sum^{k}_{r=1}y_{r}\alpha^{r}_{q,k}=x_{00},$
$\displaystyle\sum^{k}_{r=l}y_{r}\binom{r}{l}\alpha^{r-l}_{q,k}=\sum_{i,j:i+j=2l}x_{ij},$
$\displaystyle l=1,\dots,k,$ $\displaystyle X\succeq 0,$
$\displaystyle\sum_{i,j:i+j=2l-1}z_{ij}=0,$ $\displaystyle l=1,\dots,k,$
$\displaystyle\sum^{l}_{r=0}y_{r}\binom{k-r}{l-r}\alpha^{r-l}_{q,k}=\sum_{i,j:i+j=2l}z_{ij},$
$\displaystyle l=0,\dots,k,$ $\displaystyle Z\succeq 0.$

with variables $X,Z\in\mathbb{R}^{(k+1)\times(k+1)},y_{r},r=0,1,\dots,k$. Then
with this optimal detector threshold $\alpha^{\star}_{q,k}$, the false alarm
rate under the worst-case distribution of the $q_{t}$ specified by the
$k$-moments-based ambiguity set (13) is at most $\mathcal{A}$.

###### Proof.

To get an upper bound on required probability in (19), it suffices to check
for the polynomials defined in (20) to be non-negative in their respective
sets. That is, the polynomial $p(q)$ satisfies $p(q)-1\geq 0,\forall
q\in\mathcal{S}_{k}$ if and only if there exists a positive semidefinite
matrix $X\in\mathbb{R}^{(k+1)\times(k+1)}$ that satisfies the first four
constraints defined in (23). Similarly, the same polynomial $p(q)$ satisfies
$p(q)\geq 0,\forall q\in[0,\alpha_{q,k}]$ if and only if there exists a
positive semidefinite matrix $Z\in\mathbb{R}^{(k+1)\times(k+1)}$ that
satisfies the last three constraints defined in (23). Then using the results
given in Theorem 3.2 of [17], for a fixed threshold $\alpha_{q,k}$, the
optimal value of the semidefinite program in (23) yields the worst-case false
alarm rate generated by the worst-case $P_{q}\in\mathcal{P}^{q}_{k}$ in (13).
To ensure that the optimal value of (23) is equal to the desired false alarm
rate $\mathcal{A}$, it suffices to tune the threshold $\alpha_{q,k}$ using the
bisection algorithm 1, to get $\alpha^{\star}_{q,k}$. ∎

###  III-C Discussion

For $k=1,2,3$, the solution to (23) can be obtained in closed-form using
Theorem 3.3 of [17] with the corresponding moments data without explicitly
solving the SDP in (23). For $k=1,2$, we recover the Markov bound and a
strictly improved Chebyshev bound respectively as the solution of (23).
Assuming that
$\hat{\delta}_{k}>0,\alpha^{\star}_{q,k}=(1+\hat{\delta}_{k})M^{1}_{q}$, the
resulting closed form solutions for the optimal detector threshold
$\alpha^{\star}_{q,k}$ for a desired false alarm rate $\mathcal{A}$ with
$k=1,2$ are summarized in Table I with
$C^{2}_{M}=\frac{M^{2}_{q}-(M^{1}_{q})^{2}}{(M^{1}_{q})^{2}}$ and these values
can be used as $\alpha_{u}$ towards computing $\alpha^{\star}_{q,k}$ for any
$k>2$. We omit the expression for the threshold $\alpha^{\star}_{q,3}$ for the
sake of brevity.

$k$ |  $\alpha^{\star}_{q,k}$ value
---|---
1 | $\frac{M^{1}_{q}}{\mathcal{A}}$
2 | $\left(1+\sqrt{\frac{1-\mathcal{A}}{\mathcal{A}}}\,C_{M}\right)M^{1}_{q}$
TABLE I: Closed-form solutions for $\alpha^{\star}_{q,k}$ with $k=1,2.$

Further, the optimal threshold $\alpha^{\star}_{q,k}$ is tight for a given
$\mathcal{P}^{q}_{k}$ up to a tolerance $\epsilon>0$ specified by the
bisection in Algorithm 1.

The exact detector threshold (which we call $\alpha^{\star}_{q}$)
corresponding to $P_{q}$ is obtained only asymptotically (or equivalently when
true $P_{q}$ is known exactly). That is,

$\displaystyle\alpha^{\star}_{q}=\begin{cases}\lim_{k\rightarrow\infty}\alpha^{\star}_{q,k},&\text{if
all moments are known},\\\ F^{-1}_{q}(\mathcal{A}),&\text{if true
$P_{q},F_{q}$ are known}.\end{cases}$ (24)

However, while determining how many moments are needed to get a close
approximation of the exact threshold is difficult in general, we find that
significant improvements can be obtained from a small number of moments. Given
$k$ moments of $q_{t}$, the complexity of this higher-order moment based
approach involves solving the SDP given by (23) with variables
$X,Z\in\mathbb{R}^{(k+1)\times(k+1)},y_{r},r=0,1,\dots,k$. For large $k$, it
is advisable to use, e.g., the Legendre polynomial basis instead of the
standard polynomial basis for obtaining the moment-based polynomial in (23) as
the former has nice orthogonal properties that improve numerical stability.

##  IV Attack-Reachable Set Bounds

With the attacker assumed to have perfect knowledge of the system dynamics,
the Kalman filter, control inputs, measurements along with read and write
access to all the sensors at each time step, we show that the volume of the
attack-reachable set shrinks with the tightened threshold. We define a zero-
alarm attack, which generates attack sequences so that no alarms are raised
during attack. With $\Sigma^{\frac{1}{2}}_{r}$ being the symmetric square root
of $\Sigma_{r}$ and an attack input $\bar{\delta}_{t}$ such that
$\bar{\delta}^{\top}_{t}\bar{\delta}_{t}\leq\alpha$ (here
$\alpha=\alpha^{\star}_{q,k}$), the attack sequence
$\delta_{t}=-Ce_{t}-v_{t}+\Sigma^{\frac{1}{2}}_{r}\bar{\delta}_{t}$ guarantees
that no alarm is raised. Then, using a static estimator feedback
$u_{t}=K\hat{x}_{t}$ with $\bar{\delta}_{t}$, the evolution of the system
dynamics with the joint state $\xi_{t}=[x_{t},e_{t}]^{\top}$ with input
$\zeta_{t}=[w_{t},\bar{\delta}_{t}]^{\top}$ is studied. We define a reachable
set of interest, driven by the ellipsoidally bounded inputs $w_{t}$ and
$\bar{\delta}_{t}$, as

$\mathcal{R}_{x}=\left\\{x_{t}=[I_{n},\,0_{n\times n}]\xi_{t}\ \left|\
\begin{aligned} &\xi_{t+1}=\hat{A}\xi_{t}+\hat{B}\zeta_{t},\\\
&\xi_{1}=\mathbf{0},\ \bar{\delta}^{\top}_{t}\bar{\delta}_{t}\leq\alpha,\\\
&w_{t}^{\top}\Sigma_{w}^{-1}w_{t}\leq\bar{w},\ \forall
t\in\mathbb{N}\end{aligned}\right.\right\\},$ (25)

where the noise threshold $\bar{w}$ obtained using (18) satisfies,

$\sup_{P_{w}\in\mathcal{P}^{w}_{k}}P_{w}(w^{\top}_{t}\Sigma^{-1}_{w}w_{t}\leq\bar{w})=1-\mathcal{A},$
(26)

and $\hat{A}=\begin{bmatrix}A+BK&-BK\\\
0&A\end{bmatrix},\hat{B}=\begin{bmatrix}I&0\\\
I&-L\Sigma_{r}^{1/2}\end{bmatrix}$. Using the geometric approach presented in
[22] with $H_{i}=A_{cl}^{i}-A^{i}$, the reachable set of states is the
Minkowski sum of the following ellipsoidal bound

$\mathcal{R}_{x,t}(\bar{w},\alpha)=\bigoplus_{i=0}^{t-2}\mathcal{E}\left(\bar{w}A^{i}\Sigma_{w}{A^{i}}^{\top}\right)\oplus\mathcal{E}\left(\alpha
H_{i}L\Sigma_{r}L^{\top}H_{i}^{\top}\right).$ (27)

Specifically, theorem 1 of [22] provides us the exact boundary of Minkowski
sum in (27) using an analytical formula. The following corollary highlights
the effects of the tightened threshold on the size of the reachable set.

######  Corollary 5.

The volume of the reachable set $\mathcal{R}_{x,t}(\bar{w},\alpha)$ shrinks
with the tightened threshold $\alpha^{\star}_{q,k}$ for all $k\geq 1$.
Further, $\forall k_{1},k_{2}\geq 1$ and $k_{1}<k_{2}$, the following
inclusion holds
$\mathcal{R}_{x,t}(\bar{w},\alpha^{\star}_{q,k_{1}})\supseteq\mathcal{R}_{x,t}(\bar{w},\alpha^{\star}_{q,k_{2}})$.

###### Proof.

Apply (27) with $k_{1},k_{2}\geq 1$ & $k_{1}<k_{2}$ along with
$\alpha^{\star}_{q,k_{1}}>\alpha^{\star}_{q,k_{2}}$ (by Lemma 3), to see that
volume$(\mathcal{R}_{x,t}(\bar{w},\alpha^{\star}_{q,k_{2}}))\leq\text{volume}(\mathcal{R}_{x,t}(\bar{w},\alpha^{\star}_{q,k_{1}}))$.
Inclusion follows from Minkowski sum with
$\alpha^{\star}_{q,k_{1}}>\alpha^{\star}_{q,k_{2}}$. ∎

##  V Numerical Simulation

We consider an empirical system under study with the detector tuned to a false
alarm rate $\mathcal{A}=0.05$ ($5$%). We demonstrate here simulation results
when the uncertainties are zero-mean Gaussian. We compare the size of the
reachable set boundary computed using (27) for thresholds
$\alpha^{\star}_{q,1},\alpha^{\star}_{q,2},\alpha^{\star}_{q,4}$. We assume
that the modeler is unaware of the functional form of the uncertainties and
has to arrive at a detector threshold satisfying the desired false alarm rate.
The SDP in (23) was solved with $\epsilon=10^{-4}$ using SOSToolbox in Matlab
with SeDuMi solver.

$\displaystyle A$ $\displaystyle=\begin{bmatrix}0.84&0.23\\\
-0.47&0.12\end{bmatrix},B=\begin{bmatrix}0.07&-0.32\\\
0.23&0.58\end{bmatrix},C=\begin{bmatrix}1&0\\\ 2&1\end{bmatrix}$
$\displaystyle K$ $\displaystyle=\begin{bmatrix}1.404&-1.402\\\
1.842&1.008\end{bmatrix},L=\begin{bmatrix}0.0276&0.0448\\\
-0.01998&-0.0290\end{bmatrix},$ $\displaystyle\Sigma_{v}$
$\displaystyle=6I_{p},\Sigma_{w}=\begin{bmatrix}0.0225&-0.0055\\\
-0.0055&0.0100\end{bmatrix}.$

Figure 2: The reachable sets associated with the detector thresholds
$\alpha^{\star}_{q,1},\alpha^{\star}_{q,2},\alpha^{\star}_{q,4},\alpha_{\chi^{2}}$
are shown in orange, green, blue and magenta colors respectively. It is
evident from the size of the reachable set that given a desired false alarm
rate $\mathcal{A}$, the knowledge of higher-order moments results in a
tightened detector threshold and thereby restricting the attacker’s ability to
launch a larger attack. Figure 3: The moment based polynomials in orange,
green and blue bounding their respective indicator functions
$\mathbf{1}_{(\alpha^{\star}_{q,1},\infty)},\mathbf{1}_{(\alpha^{\star}_{q,2},\infty)},\mathbf{1}_{(\alpha^{\star}_{q,4},\infty)}$
in shaded orange, green and blue colors are shown here. Clearly,
$\alpha^{\star}_{q,1}$ is very conservative and with $k=4$, the threshold
$\alpha^{\star}_{q,4}$ starts getting closer to the true threshold
$\alpha_{\chi^{2}}$ given by (16).

When the noises $w_{t}$ and $v_{t}$ are truly Gaussian, it is evident from
Fig. 2 that the reachable set corresponding to the thresholds
$\alpha^{\star}_{q,1}=\alpha^{\star}_{r,2}=40$ is conservative and ensures
that the false alarm does not exceed 5% but this also provides the attacker
with the ability to launch a larger attack. However, with the knowledge of
additional moments, the detector thresholds $\alpha^{\star}_{q,2}=10.7684$,
$\alpha^{\star}_{q,4}=9.1315$ get tightened with
$\alpha^{\star}_{q,4}\leq\alpha^{\star}_{q,2}\leq\alpha^{\star}_{q,1}$ as
shown in Fig. 3. Further, this threshold tightening limits the attacker’s
ability to launch a larger attack which is depicted through the reachable sets
corresponding to the thresholds
$\alpha^{\star}_{q,1},\alpha^{\star}_{q,2},\alpha^{\star}_{q,4}$ as shown in
Fig. 2. Subsequently, the false alarm rate corresponding to the threshold
$\alpha_{\chi^{2}}$ was 5% as expected and with the thresholds
$\alpha^{\star}_{q,4},\alpha^{\star}_{q,2},\alpha^{\star}_{q,1}$, it dropped
to $1\%,0.45\%,0\%$ respectively. When noises are truly multi-variate
Laplacian (which has heavier tails than normal distribution with same mean and
covariance), it resulted in thresholds $\alpha^{\star}_{q,1}=39.83$,
$\alpha^{\star}_{q,2}=17.23,\alpha^{\star}_{q,4}=16.54$ with false alarm rates
$0\%,0.9\%,1\%$ respectively. Thus, no matter what distributions satisfying
(3), (4) govern the noises $w_{t}$, $v_{t}$ respectively, the inclusion of
higher-order moments restricts the attacker’s potential impact through a
tightened detector threshold.

##  VI Conclusion & Future Outlook

We have proposed a distributionally robust approach to form the $k$ moments
based ambiguity set for the detection measure data and used it to tune the
anomaly detectors for a desired false alarm rate. We found a detector
threshold which guaranteed that the false alarm rate did not exceed a desired
value using a semidefinite program. We have demonstrated the effectiveness of
our proposed approach with a numerical example. Further, our approach using
higher-order moments restricted the attacker’s potential impact. Future works
include addressing the problems associated with the data-driven formulation
with moment estimation uncertainty and securing nonlinear cyberphysical
systems with distributionally robust unscented Kalman filter based state
estimation.

## References

  * [1] C. Murguia and J. Ruths, “On model-based detectors for linear time-invariant stochastic systems under sensor attacks,” _IET Control Theory & Applications_, vol. 13, no. 8, pp. 1051–1061, 2019.
  * [2] D. Umsonst and H. Sandberg, “Anomaly detector metrics for sensor data attacks in control systems,” in _2018 Annual American Control Conference (ACC)_. IEEE, 2018, pp. 153–158.
  * [3] V. Renganathan, N. Hashemi, J. Ruths, and T. H. Summers, “Distributionally robust tuning of anomaly detectors in cyber-physical systems with stealthy attacks,” in _2020 American Control Conference (ACC)_. IEEE, 2020, pp. 1247–1252.
  * [4] M. J. Khojasteh, A. Khina, M. Franceschetti, and T. Javidi, “Learning-based attacks in cyber-physical systems,” _IEEE Transactions on Control of Network Systems_ , 2020.
  * [5] D. Li and S. Martínez, “High-confidence attack detection via wasserstein-metric computations,” _IEEE Control Systems Letters_ , vol. 5, no. 2, pp. 379–384, 2021.
  * [6] K. Zetter, _Countdown to Zero Day: Stuxnet and the launch of the world’s first digital weapon_. Broadway books, 2014\.
  * [7] A. Cherepanov, “Win32/industroyer, a new threat for industrial control systems,” _White Paper. ESET_ , 2017.
  * [8] W. Wiesemann, D. Kuhn, and M. Sim, “Distributionally robust convex optimization,” _Operations Research_ , vol. 62, no. 6, pp. 1358–1376, 2014\.
  * [9] J. Goh and M. Sim, “Distributionally robust optimization and its tractable approximations,” _Operations research_ , vol. 58, no. 4-part-1, pp. 902–917, 2010.
  * [10] P. M. Esfahani and D. Kuhn, “Data-driven distributionally robust optimization using the wasserstein metric: Performance guarantees and tractable reformulations,” _Mathematical Programming_ , vol. 171, no. 1-2, pp. 115–166, 2018.
  * [11] N. Hashemi, C. Murguia, and J. Ruths, “A comparison of stealthy sensor attacks on control systems,” in _2018 Annual American Control Conference (ACC)_ , June 2018, pp. 973–979.
  * [12] N. Hashemi and J. Ruths, “Generalized chi-squared detector for lti systems with non-gaussian noise,” in _2019 American Control Conference (ACC)_. IEEE, 2019, pp. 404–410.
  * [13] J.-B. Lasserre, _Moments, positive polynomials and their applications_. World Scientific, 2010, vol. 1.
  * [14] J. B. Lasserre, “Global optimization with polynomials and the problem of moments,” _SIAM Journal on optimization_ , vol. 11, no. 3, pp. 796–817, 2001\.
  * [15] R. E. Curto, L. A. Fialkow, and H. M. Möller, “The extremal truncated moment problem,” _Integral Equations and Operator Theory_ , vol. 60, no. 2, pp. 177–200, 2008.
  * [16] J. W. Helton and J. Nie, “A semidefinite approach for truncated k-moment problems,” _Foundations of Computational Mathematics_ , vol. 12, no. 6, pp. 851–881, 2012.
  * [17] D. Bertsimas and I. Popescu, “Optimal inequalities in probability theory: A convex optimization approach,” _SIAM Journal on Optimization_ , vol. 15, no. 3, pp. 780–804, 2005.
  * [18] J. Shawe-Taylor and N. Cristianini, “Estimating the moments of a random vector with applications,” in _19th Colloque sur le traitement du signal et des images, FRA, 2003_. GRETSI, Groupe d’Etudes du Traitement du Signal et des Images, 2003.
  * [19] E. Delage and Y. Ye, “Distributionally robust optimization under moment uncertainty with application to data-driven problems,” _Operations research_ , vol. 58, no. 3, pp. 595–612, 2010.
  * [20] L. Vandenberghe, S. Boyd, and K. Comanor, “Generalized chebyshev bounds via semidefinite programming,” _SIAM review_ , vol. 49, no. 1, pp. 52–64, 2007\.
  * [21] X. Chen, “A New Generalization of Chebyshev Inequality for Random Vectors,” _preprint arXiv:0707.0805_ , Jul. 2007.
  * [22] N. Hashemi and J. Ruths, “Co-design for security and performance: Geometric tools,” _preprint arXiv:2006.08739_ , 2020.


# The International Pulsar Timing Array checklist for the detection of
nanohertz gravitational waves

Bruce Allen  Max Planck Institute for Gravitational Physics, Leibniz
Universität Hannover, Callinstrasse 38, D-30167 Hannover, Germany  Sanjeev
Dhurandhar  Inter University Centre for Astronomy & Astrophysics, Ganeshkhind,
Pune - 411 007, India  Yashwant Gupta  National Centre for Radio Astrophysics,
Pune University Campus, Pune 411007, India  Maura McLaughlin  West Virginia
University Department of Physics and Astronomy, Morgantown, WV, 26501, USA
Priyamvada Natarajan  Department of Astronomy, 52 Hillhouse Avenue, New Haven,
CT 06511  Black Hole Initiative, 20 Garden Street, Cambridge, MA 02138  Ryan
M. Shannon  Centre for Astrophyics and Supercomputing, Swinburne University of
Technology, Hawthorn, VIC, 3122, Australia  OzGrav: The ARC Centre of
Excellence for Gravitational Wave Discovery  Eric Thrane  School of Physics
and Astronomy, Monash University, Clayton VIC 3800, Australia  OzGrav: The ARC
Centre of Excellence for Gravitational Wave Discovery, Clayton VIC 3800,
Australia  Alberto Vecchio  School of Physics and Astronomy & Institute for
Gravitational Wave Astronomy, University of Birmingham, Birmingham, B15 2TT,
UK

###### Abstract

Pulsar timing arrays (PTAs) provide a way to detect gravitational waves at
nanohertz frequencies. In this band, the most likely signals are stochastic,
with a power spectrum that rises steeply at lower frequencies. Indeed, the
observation of a common red noise process in pulsar-timing data suggests that
the first credible detection of nanohertz-frequency gravitational waves could
take place within the next few years. The detection process is complicated by
the nature of the signals and the noise: the first observational claims will
be statistical inferences drawn at the threshold of detectability. To
demonstrate that gravitational waves are creating some of the noise in the
pulsar-timing data sets, observations must exhibit the Hellings and Downs
curve—the angular correlation function associated with gravitational waves—as
well as demonstrating that there are no other reasonable explanations. To
ensure that detection claims are credible, the International Pulsar Timing
Array (IPTA) has a formal process to vet results prior to publication. This
includes internal sharing of data and processing pipelines between different
PTAs, enabling independent cross-checks and validation of results. To oversee
and validate any detection claim, the IPTA has also created an eight-member
Detection Committee (DC) which includes four independent external members.
IPTA members will only publish their results after a formal review process has
concluded. This document is the initial DC checklist, describing some of the
conditions that should be fulfilled by a credible detection. At the present
time none of the PTAs have a detection claim; therefore this document serves
as a road map for the future.

##  1 Introduction

A nanohertz-frequency stochastic background of gravitational waves creates
corresponding low-frequency fluctuations in pulsar timing residuals (Sazhin,
1978; Detweiler, 1979; Foster & Backer, 1990). Such fluctuations may first be
inferred with the same spectrum in all PTA pulsars Hobbs et al. (2009). If
this common-spectrum process arises from a gravitational-wave background, it
will have statistically consistent amplitude and spectral shape in every
pulsar. However, the intrinsic rotation noise of pulsars, and the noise
generated by pulse propagation through the interstellar medium (ISM) are not
fully understood Groth (1975); Shannon & Cordes (2010); Melatos & Link (2014).
This means that the observation of a common (red) noise process is not, by
itself, compelling evidence of gravitational waves.

One way to distinguish gravitational waves from other sources of pulsar timing
fluctutations was proposed by Hellings & Downs (1983). They show that low-
frequency gravitational waves create a pattern of angular correlations between
pulsars in different parts of the sky. These correlations have quadrupolar
signature described by the Hellings and Downs curve (Hellings & Downs, 1983).
A clear observation of this pattern of angular correlations is one way to
distinguish a gravitational-wave background from intrinsic pulsar spin noise,
interstellar-medium effects , observatory clock errors, ephemeris errors, and
other sources of noise Tiburzi et al. (2016).

For several decades, PTAs have placed increasingly stringent constraints on
the amplitude of the nanohertz-frequency gravitational-wave background through
spectral analyses of pulsar timing measurements. Recently, three PTA
collaborations—the European Pulsar Timing Array (EPTA, Kramer & Champion,
2013), the North American Nanohertz Observatory for Gravitational Waves
(NANOGrav, McLaughlin, 2013), and the Parkes Pulsar Timing Array (PPTA,
Manchester et al., 2013) all reported the detection of a common-spectrum
process in their current datasets (Arzoumanian et al., 2020; Chen et al.,
2021; Goncharov et al., 2021). Here, “common spectrum” means that the
fluctuations in different pulsar arrival times are described by the same
spectrum, but does _not_ imply that the fluctuations have a common source (in
which case they would be correlated between different pulsars). Indeed, this
same common-spectrum process is also detected in the most recent International
Pulsar Timing Array (IPTA) dataset, consisting of data from these three
regional collaborations (Antoniadis et al., 2022).

The measured amplitude and spectral index of this common-spectrum process is
consistent in all four datasets at the 2$\sigma$ level (Antoniadis et al.,
2022) It is also consistent with some theoretical predictions for a
gravitational-wave background from a cosmological population of supermassive
black hole binaries (e.g., Begelman et al., 1980; Rajagopal & Romani, 1995;
Sesana, 2013), though, given the current level of understanding of the
uncertainties in the modeling of the growth and evolution of the supermassive
black hole population,a wide range of predictions have been made (e.g., Sesana
et al., 2008; Ravi et al., 2012; Zhu et al., 2019). However, angular
correlations with the characteristic quadrupolar signature have not yet been
published. Therefore, so far, it has not been possible to conclude that the
common-spectrum process is due to an astrophysical/cosmological gravitational-
wave background. However, if the currently observed common process is due to
gravitational waves, simulations predict that the spatial correlations could
be detectable at the 4–5$\sigma$ level in the NANOGrav 15-yr dataset (Pol et
al., 2021) and similar significance detections may be soon possible in other
PTA data sets.

In anticipation of evidence for a gravitational-wave background in one or
multiple PTA datasets currently under analysis, the IPTA collaboration (which
includes the EPTA, NANOGrav, PPTA and the Indian Pulsar Timing Array) convened
a Detection Committee comprising of members from all four PTAs, and four
independent experts. The Detection Committee is tasked with developing a
detection checklist to help verify candidate gravitational-wave signals. This
document is the result of the work of that Committee. At the present time none
of the PTAs have a detection claim. This checklist has not been applied to the
submitted papers, which therefore have not been vetted by the Detection
Committee.

##  2 Detection checklist

The detection checklist consists of three subsections. In 2.1, we describe
checklist items designed to ensure that Hellings-Downs correlations are
statistically significant. In 2.2, we describe checklist items that provide
consistency checks to ensure that the signal and noise model are well
specified, and that the signal cannot be easily explained by systematic error.
Finally, in 2.3, we describe checklist items that ensure that results are
cross-checked and vetted by independent experts.

###  2.1 Statistical significance

  * $\square$

The Hellings-Downs signal is evident in the data with $5\sigma$ significance.
There are different ways to show this. We provide two examples: one Bayesian
and one frequentist.

Bayesian formulation. Construct the $5\sigma$, highest posterior density
credible interval for the gravitational-wave amplitude $A$. Show that this
credible interval excludes $A=0$. The posterior is constructed using only
cross-power (no auto-power) so that this detection statement is not influenced
by the presence of (quasi-) common red noise, which may or may not be due to
gravitational waves. When constructing this posterior, analysts should
marginalize over all relevant sources of uncertainty including pulsar noise
models and astrophysical uncertainty in the signal model.

Frequentist formulation. Construct an estimator for the amplitude of the
Hellings-Downs amplitude $Y_{\text{HD}}$ with associated uncertainty
$\sigma_{Y}$ using only point estimates for the angular correlation function
from pairs of distinct pulsars and the associated error bars. The following is
a frequentist detection statistic for the Hellings-Downs signal:

$\displaystyle\rho_{\text{HD}}=Y_{\text{HD}}/\sigma_{Y}.$ (1)

Calculate the null distribution of $\rho_{\text{HD}}$ under the assumption
that the signal and noise models are correctly specified and that no
correlation is present. The observed value of $\rho_{\text{HD}}$ occurs with a
probability of $p\leq 3\times 10^{-7}$ in the null distribution ($\geq
5\sigma$ confidence).

  * $\square$

The statistical significance computed in the previous bullet is consistent
with estimates from bootstrap methods. At least $N=1000$ (preferably more)
quasi-independent noise realisations from phase-scrambling and/or sky-
scrambling have been analysed.111These boot-strap noise realisations should
take into account the relative quality of each pulsar as per Cornish & Sampson
(2016). Evidence needs to be provided to demonstrate that the required number
of quasi-independent realisations has been achieved. Repeating the previous
test, but with boot-strap noise realizations, yields a detection that is
consistent with a $\geq 5\sigma$ detection. There are no instances of boot-
strap noise with false-positive detections that are more significant than the
signal present in the data.222In the Bayesian formulation, there are no boot-
strap realizations where $A=0$ is excluded with higher credibility than it is
in the data. In the frequentist formulation, there are no boot-strap
realizations with larger values of $\rho_{\text{HD}}$.

###  2.2 Consistency checks

The data must be consistent with the signal and noise models. The checklist
items in this section are designed to ensure that the angular correlation
function is consistent with a Hellings-Downs signal and clearly quadrupolar in
nature. It also needs to be demonstrated that the significance is not overly
reliant on a small number of pulsars

  * $\square$

The angular correlation function is consistent with the Hellings-Downs curve.
Compute a $p$-value under the null hypothesis that the measured angular
correlation function is consistent with the Hellings-Downs curve for an
isotropic stochastic background. It can be calculated using a $\chi^{2}$
statistic, which sums over pulsar pairs (or angular bins):

$\displaystyle\chi^{2}(h)=$
$\displaystyle(y_{\alpha}-h^{2}\mu_{\alpha})^{\dagger}\,C^{-1}_{\alpha\beta}(h)\,(y_{\beta}-h^{2}\mu_{\beta}).$
(2)

Repeated indices imply summation. Here, $y_{\alpha}$ is the cross correlation
for pulsar pair $\alpha$, $\mu_{\alpha}$ is the Hellings-Downs curve evaluated
at the angular separation for pair $\alpha$, $h$ is the gravitational-wave
background amplitude, and $C_{\alpha\beta}$ is the covariance matrix. The
value of $h^{2}$ may be set “externally,” i.e., based on autocorrelation
measurements, or “internallly” by minimizing $\chi^{2}$, which reduces the
number of degrees of freedom by one. The $\chi^{2}$ value is used to calculate
a $p$-value. If the data are consistent with the Hellings-Downs hypothesis, we
expect $p\approx 50\%$. A small $p$-value $<0.02$ indicates that the model
does not provide an adequate fit.

In principle, the covariance matrix in Eq. 2 incorporates the effects of
measurement uncertainty arising from pulsar and measurement noise and cosmic
variance arising from the random amplitudes, phases and locations of
gravitational-wave sources; see Allen & Romano (2022). If it is included, the
cosmic variance contributions should be consistent with studies published in
peer-reviewed journals.

It is useful to contrast this checklist item (related to goodness of fit) with
the checklist items in 2.1 related to significance. We provide examples in
Table 1 showing different scenarios. The first rows illustrate two different
failure modes while the final row provides an example of a detection claim
that passes.

Example situation |  Significance: Credibility with which we exclude $A=0$ |  Consistency: $p$-value under the null hypothesis (that the data are described by the model) |  pass/fail
---|---|---|---
Significant, but misspecified |  $1.32\times 10^{-7}$ |  $0.003$ |  fail
Not sufficiently significant, but adequately specified |  $3.25\times 10^{-3}$ |  0.67 |  fail
Significant, adequately specified |  $1.32\times 10^{-7}$ |  $0.67$ |  pass
Table 1: Example situations illustrating how a detection claim is expected to
pass or fail based on the statistical significance of the signal and the
consistency of the data with the model. A failure is marked with red while a
pass is marked with green.

  * $\square$

The signal is clearly quadrupolar. In particular, the pure Hellings-Downs
correlation is preferred over a model consisting of a pure monopolar
correlation and/or a pure dipolar correlation with a Bayes factor of $>100$.
If “quadrupole + monopole” or “quadrupole + dipole” hypotheses are
significantly preferred over the pure Hellings-Downs hypothesis (with Bayes
factor $\gtrsim 100$), then there are likely still artifacts in the data that
are not yet correctly modeled (e.g., clock errors, ephemeris errors, etc.). In
such cases, more work may be required before we can establish a detection.

  * $\square$

The signal is present in more than a small number of pulsars. In particular,
the number of effective pulsar pairs $n_{\text{eff}}$ is at least $20$. The
effective number of pulsar pairs is

$\displaystyle n_{\text{eff}}=$
$\displaystyle\frac{\left(\sum_{k=1}^{m}w_{k}\right)^{2}}{\sum_{k=1}^{m}w_{k}^{2}}$
$\displaystyle=$
$\displaystyle\frac{\left(\sum_{k=1}^{m}\widehat{\sigma}_{k}^{-2}\right)^{2}}{\sum_{k=1}^{m}\widehat{\sigma}_{k}^{-4}}.$
(3)

Here, $w_{k}$ is the “weight” of pulsar pair $k$ while $\widehat{\sigma}_{i}$
is the angular correlation function uncertainty for pulsar pair $k$. The total
number of pulsar pairs is $m$. (We use the hat to differentiate the
uncertainty associated with a pulsar pair $\widehat{\sigma}_{k}$ from the
uncertainty associated with a cosine angular separation bin $\sigma_{i}$.)
Since some pairs are more important than others, $n_{\text{eff}}<m$. If
$n_{\text{eff}}<20$, then the measurement is dominated by a small number of
pulsars such that the Hellings-Downs curve will not be clearly evident through
visual inspection of the angular correlation function. This check ensures that
the Hellings-Downs correlation is visible when the data are binned. At least
seven well-timed pulsars are required to fulfill this requirement. This item
is similar in spirit to the “drop-out factors,” which pulsar timing arrays use
to quantify the relative importance of different pulsars.

  * $\square$

The result is consistent with previously published analyses by the same PTA.
The inferred gravitational-wave background is consistent with previous upper-
limit papers written using a subset of the currently available data. If the
inferred gravitational-wave background is inconsistent with previously
published papers, there is a convincing explanation to account for the
discrepancy.

###  2.3 Independent vetting of results

In addition to the checklist items described above, it is useful for the
analysis to be independently verified by domain experts. The following
checklist items are designed to facilitate this.

  * $\square$

The detection team will make available to the other PTAs everything required
to reproduce the detection. This includes:

    * •

Standard tempo2-format data including TOA tim files and ephemerides par files.

    * •

The noise model in either Enterprise or temponest format.

    * •

A technical note detailing all the assumptions needed to reproduce the
detection.

    * •

Weekly “office hours” to field questions from other analysis teams.

The Detection Committee may request additional code or documentation in order
to support the reproduction and verification of their main results.

The process by which NANOGrav, EPTA, and PPTA coordinate their papers is
referred to as “the 3P+ framework.” As per 3P+ rules, shared data may not be
used for any published papers without permission; it is _only_ for checking.
Analysis teams will provide data for these checks as soon as practical as part
of the 3P+ process. We recommend that the IPTA data combination working group
coordinate the sharing of data products and expertise between the constituent
PTAs. We note that the data combination working group has started some of this
work in preparation for Data Release 3.

Members of the IPTA will have at least six weeks from a detection/evidence
claim by any PTA to inspect the data. The data products include TOAs, pulsar
ephemerides, preferred noise models, and, if possible, intermediate data
products such as posterior chains.

During this time, independent teams are asked to interact with the Detection
Committee and to submit brief reports (and supporting analysis materials, as
far as practical) to the Detection Committee stating whether or not they are
able to find an alternative explanation that could explain the data at least
as well as the gravitational-wave hypothesis. The independent teams are
encouraged to share preliminary findings early with the Detection Committee.
These reports (and supporting materials) will be made available to the IPTA
via the IPTA 3P+ Committee. The Detection Committee may request the detection
team to provide a reply to one or more reports. The Detection Committee will
share all reports it receives with the 3P+ Committee. The Detection Committee
may amend the checklist based on developments during this period.

###  2.4 Initiating the detection procedure

The procedure for a PTA to initiate a detection/evidence claim is as follows:

  1. 1.

Prepare a paper draft and send it to the Detection Committee.

  2. 2.

Prepare a “response to the detection checklist” document and send it to the
Detection Committee. This technical note should address every item in the
detection checklist provided above. The reply should be relatively succinct
for readability; less than one page per checklist item. If multiple pages of
supporting material are required for a single checklist item, some of it can
be placed in an appendix.

  3. 3.

Present the detection case to the Detection Committee. The presentation will
be scheduled approximately two weeks after the Detection Committee has
receives both the paper draft and the response to the detection checklist.
This will provide the Detection Committee sufficient time to read both
documents and formulate their questions.The Detection Committee may request a
follow-up meeting for additional questions and/or to discuss additional
investigations that may be deemed necessary.

  4. 4.

The Detection Committee issues a recommendation stating whether or not they
endorse the detection claim. The Detection Committee will aim to provide this
recommendation within one month of the final presentation and/or completion of
any additional investigations. The recommendation will be written as a
technical note. If the detection claim is not endorsed, the Detection
Committee will specify the aspect that was unconvincing.

The Detection Committee shall endeavour to reach a consensus opinion. However,
if that proves impossible, the Committee will draft a report based on the
majority opinion. The perspective of the dissenting minority will be included
as an appendix. In the event of a tie vote, the majority opinion shall be
determined by the Detection Committee Chair. This report shall be made
publicly available.

##  Appendix A Future work

The checklist above reflects the short-time scale over which the first
paper(s) from the current analyses may become available and the wish of PTA
teams to submit them quickly for publication. Taking these factors into
account the IPTA Detection Committee has distilled down a list of absolutely
necessary checks, which are presented above, but has discussed a much longer
list of additional checks that it would strongly encourage PTAs to carry out
(partially or in full). Many suggestions have also been received by the
Detection Committee upon circulation of the detection checklist for comments
from the IPTA. A revised and updated checklist may well incorporate more
detailed checks (e.g., marginalizing over different solar system ephemerides,
for instance).

## References

  * Allen & Romano (2022) Allen, B., & Romano, J. D. 2022.  https://arxiv.org/abs/2208.07230
  * Antoniadis et al. (2022) Antoniadis, J., et al. 2022, MNRAS, 510, 4873
  * Arzoumanian et al. (2020) Arzoumanian, Z., et al. 2020, ApJ, 905, L34
  * Begelman et al. (1980) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307, doi: 10.1038/287307a0
  * Chen et al. (2021) Chen, S., et al. 2021, MNRAS, 508, 4970
  * Cornish & Sampson (2016) Cornish, N. J., & Sampson, L. 2016, Phys. Rev. D, 93, 104047
  * Detweiler (1979) Detweiler, S. 1979, ApJ, 234, 1100, doi: 10.1086/157593
  * Foster & Backer (1990) Foster, R. S., & Backer, D. C. 1990, ApJ, 361, 300, doi: 10.1086/169195
  * Goncharov et al. (2021) Goncharov, B., Shannon, R. M., Reardon, D. J., et al. 2021, ApJ, 917, L19
  * Groth (1975) Groth, E. J. 1975, ApJS, 29, 453, doi: 10.1086/190354
  * Hellings & Downs (1983) Hellings, R. W., & Downs, G. S. 1983, ApJ, 265, L39
  * Hobbs et al. (2009) Hobbs, G., Jenet, F., Lee, K. J., et al. 2009, MNRAS, 394, 1945, doi: 10.1111/j.1365-2966.2009.14391.x
  * Kramer & Champion (2013) Kramer, M., & Champion, D. J. 2013, Classical and Quantum Gravity, 30, 224009, doi: 10.1088/0264-9381/30/22/224009
  * Manchester et al. (2013) Manchester, R. N., Hobbs, G., Bailes, M., et al. 2013, PASA, 30, e017, doi: 10.1017/pasa.2012.017
  * McLaughlin (2013) McLaughlin, M. A. 2013, Classical and Quantum Gravity, 30, 224008, doi: 10.1088/0264-9381/30/22/224008
  * Melatos & Link (2014) Melatos, A., & Link, B. 2014, MNRAS, 437, 21, doi: 10.1093/mnras/stt1828
  * Pol et al. (2021) Pol, N. S., et al. 2021, ApJ, 911, L34
  * Rajagopal & Romani (1995) Rajagopal, M., & Romani, R. W. 1995, ApJ, 446, 543, doi: 10.1086/175813
  * Ravi et al. (2012) Ravi, V., Wyithe, J. S. B., Hobbs, G., et al. 2012, ApJ, 761, 84, doi: 10.1088/0004-637X/761/2/84
  * Sazhin (1978) Sazhin, M. V. 1978, Soviet Ast., 22, 36.  https://articles.adsabs.harvard.edu/pdf/1978SvA....22...36S
  * Sesana (2013) Sesana, A. 2013, MNRAS, 433, L1, doi: 10.1093/mnrasl/slt034
  * Sesana et al. (2008) Sesana, A., Vecchio, A., & Colacino, C. N. 2008, MNRAS, 390, 192, doi: 10.1111/j.1365-2966.2008.13682.x
  * Shannon & Cordes (2010) Shannon, R. M., & Cordes, J. M. 2010, ApJ, 725, 1607, doi: 10.1088/0004-637X/725/2/1607
  * Tiburzi et al. (2016) Tiburzi, C., Hobbs, G., Kerr, M., et al. 2016, MNRAS, 455, 4339, doi: 10.1093/mnras/stv2143
  * Zhu et al. (2019) Zhu, X.-J., Cui, W., & Thrane, E. 2019, MNRAS, 482, 2588, doi: 10.1093/mnras/sty2849


# No Efficient Disjunction or Conjunction of Switch-Lists

Stefan Mengel  Univ. Artois, CNRS, Centre de Recherche en Informatique de Lens
(CRIL)

###### Abstract

It is shown that disjunction of two switch-lists can blow up the
representation size exponentially. Since switch-lists can be negated without
any increase in size, this shows that conjunction of switch-lists also leads
to an exponential blow-up in general.

##  1 Introduction

Switch-lists are a representation language for Boolean functions introduced in
[7], strongly related to interval representations [6]. The idea is to write
the values of a Boolean function $f$ on all lexicographically ordered inputs
in a value table. Then, to encode $f$, it suffices to remember the value of
$f$ on the first input and the inputs at which the value of $f$ changes from
that of its predecessor. The resulting data structure is called a _switch-
list_ representation of $f$. Clearly switch list representations can be far
more succinct than truth tables, e.g. for constant functions.

To systematically understand the properties of switch-lists beyond this,
Chromý and Čepek [1] analyzed them in the context of the so-called _knowledge
compilation map_. This framework, introduced in the ground-breaking work of
Darwiche and Marquis [2] gives a list of standard properties which should be
analyzed for languages used in the area of knowledge compilation along
different axes: succinctness, queries and transformations. The idea of the
knowledge compilation map has had a huge influence and the approach of [2] is
widely applied in knowledge compilation, see e.g. [5, 3, 4] for a very small
sample.

Chromý and Čepek [1] analyzed switch-lists along the properties of the
knowledge compilation map and got a nearly complete picture. It turns out that
switch-lists, while being generally much more succinct than truth tables, have
many of their good properties. In particular, all of the queries in [2] (e.g.
consistency, entailement and counting) can be answered in polynomial time on
switch-lists and nearly all of the transformation can be performed
efficiently. The only exception is that [1] leaves open if switch-lists are
closed under bounded disjunction and bounded conjunction, i.e., given two
Boolean functions $f_{1}$ and $f_{2}$ represented by switch-lists, can one
compute a switch-list representation of $f_{1}\lor f_{2}$, resp. $f_{1}\land
f_{2}$, in polynomial time. It is shown here that this is not the case: there
are Boolean functions $f_{1}$, $f_{2}$ such that any switch list
representation of $f_{1}\lor f_{2}$ is exponentially larger than those of
$f_{1}$ and $f_{2}$. This completes the analysis of switch-lists along the
criteria of the knowledge compilation map and shows that (bounded) disjunction
and conjunction are the only “bad” transformations of switch-lists, as there
is no hope for a polynomial-time procedure in this case.

##  2 Preliminaries

Let $f$ be a Boolean function in the $n$ variables $\\{x_{1},\ldots,x_{n}\\}$.
Fix an order $\pi$ of $\\{1,\ldots,n\\}$. Then, the assignment
$a:\\{x_{1},\ldots,x_{n}\\}\rightarrow\\{0,1\\}$ can be identified with the
number $b(a)\in\\{0,\ldots,2^{n}-1\\}$ by identifying $a$ with
$b(a):=\sum_{i=1}^{n}a(x_{\pi(i)})2^{i-1}$. This allows to write $a\prec
a^{\prime}$ if and only if $b(a)<b(a^{\prime})$. The intuition behind all this
is that the assignments are written in lexicographical order with respect to
$\pi$ and then $a\prec a^{\prime}$ if and only if $a$ appears before
$a^{\prime}$.

A _switch_ of the function $f$ with respect to $\pi$ is a number
$b\in\\{1,\ldots,2^{n}-1\\}$ such that $f(b)\neq f(b-1)$ (note that here the
identification of numbers and assignments to $\\{x_{1},\ldots,x_{n}\\}$
depending on the order $\pi$ is used). The _switch-list_ representation of $f$
with respect to $\pi$ consist of the value $f(0)$ and an ordered list of all
switches of $f$ with respect to $\pi$. Note that, for fixed $\pi$ the switch-
list representation uniquely determines $f$ and $f$ uniquely determines the
switch-list representation.

The _size_ of a switch-list representation is defined as $n$ times the number
of switches which corresponds roughly to the natural encoding size111We do not
take into account the size of an encoding of $\pi$ in this since it is the
same for all switchlists in $n$ variables and thus would only complicate the
notion without giving any insights.. Note that the size depends strongly on
the order $\pi$.

Following Darwiche and Marquis [2], switch-lists are said to satisfy bounded
disjunction (resp. bounded conjunction) if there is a polynomial-time
algorithm that, given two switch-list representations of functions
$f_{1},f_{2}$, computes a switch-list representation of $f_{1}\lor f_{2}$
(resp. $f_{1}\land f_{2}$). Chromý and Čepek [1] also considered the
restricted version of bounded disjunction (resp. conjunction) in which one
assumes that the involved functions $f_{1},f_{2}$ depend on the same set of
variables.

##  3 The Proof

Let $n\in\mathbb{N}$ be even. Consider the functions
$f_{1}(x_{1},\ldots,x_{n}):=\Big{(}\bigwedge_{i=1}^{n/2}x_{i}\Big{)}\lor\Big{(}\bigwedge_{i=1}^{n}\neg
x_{i}\Big{)}$ and
$f_{2}(x_{1},\ldots,x_{n}):=\Big{(}\bigwedge_{i=n/2+1}^{n}x_{i}\Big{)}\lor\Big{(}\bigwedge_{i=1}^{n}\neg
x_{i}\Big{)}$.

######  Observation 1.

There are switch-list representations for $f_{1}$ and $f_{2}$ with at most two
switches.

###### Proof.

Only give the argument for $f_{1}$ is given as that for $f_{2}$ is completely
analogous. Fix any order $\pi$ in which the variables $x_{1},\ldots,x_{n/2}$
come before those in $x_{n/2+1},\ldots,x_{n}$. An assignment is a model of
$f_{1}$ if and only if it maps all variables to $0$ or it maps
$x_{1},\ldots,x_{n/2}$ to $1$. So all models different from $0$ lie in the
interval $[\sum_{j=n/2+1}^{n}2^{j-1},\sum_{j=1}^{n}2^{j-1}]$. Note that this
interval lies at the end of the order of all assignments. So for these models,
$f_{1}$ only has one switch at the beginning of the interval. To represent
$f_{1}$ with a switch-list one only needs one additional switch directly after
$0$. ∎

######  Proposition 1.

The function $f_{1}\lor f_{2}$ needs at least $2^{n/2+1}-3$ switches in any
switch-list representation.

###### Proof.

Let $X_{1}:=\\{x_{1},\ldots,x_{n/2}\\}$ and
$X_{2}:=\\{x_{n/2+1},\ldots,x_{n}\\}$. Fix any variable order $\pi$ of
$X_{1}\cup X_{2}$ and let $\preceq$ denote the lexicographical order with
respect to $\pi$. The last variable of $\pi$ is either in $X_{1}$ or in
$X_{2}$. Without loss of generality, assume that it is in $X_{2}$ and that the
last variable in $\pi$ is $x_{n}$.

For every assignment $a$ to $X_{1}$, two extensions $e_{0}(a)$ and $e_{1}(a)$
to $X_{1}\cup X_{2}$ are constructed as follows: on $X_{1}$, the assignments
$e_{0}(a)$ and $e_{1}(a)$ are both identical to $a$; all variables in
$X_{2}\setminus\\{x_{n}\\}$ are assigned $1$ and $x_{n}$ is assigned $0$ in
$e_{0}(a)$ and $1$ in $e_{1}(a)$. Let $\pi_{1}$ be the order $\pi$ restricted
to $X_{1}$ and let $\preceq_{1}$ be the order of the assignments to $X_{1}$
with respect to $\pi_{1}$. Then for two assignments $e_{i}(a_{1})$ and
$e_{j}(a_{2})$ it holds that $e_{i}(a_{1})\prec e_{j}(a_{2})$ if and only if
$a_{1}\prec_{1}a_{2}$; or $a_{1}=a_{2}$ and $i<j$. Note that none of the
assignments of the form $e_{i}(a)$ is the constant $0$-assignment, so
$e_{i}(a)$ satisfies $f_{1}\lor f_{2}$ if and only if it satisfies
$\left(\bigwedge_{i=1}^{n/2}x_{i}\right)\lor\left(\bigwedge_{i=n/2+1}^{n}x_{i}\right)$.

Now let $a_{1},\ldots,a_{2^{n/2}-1}$ be the assignments to $X_{1}$ different
from constant $1$-assignment given in the order $\preceq_{1}$. Then the
resulting sequence

$\displaystyle
e_{0}(a_{1}),e_{1}(a_{1}),\ldots,e_{0}(a_{2^{n/2}-1}),e_{1}(a_{2^{n/2}-1})$
(1)

is in lexicographical order as well. Note that because none of the $a_{i}$ is
the constant $1$-assignment, it holds that for every $i\in[2^{n/2}-1]$ that
$e_{1}(a_{i})$ is a model of $f_{1}\lor f_{2}$ while $e_{0}(a_{i})$ is not.
Thus there must be switches between each pair of consecutive elements of the
sequence (1). So there must be at least $2\times(2^{n/2}-1)-1=2^{n/2+1}-3$
switches in the switch-list representation of $f_{1}\lor f_{2}$ with respect
to the order $\pi$. ∎

The main result of this paper follows directly.

######  Theorem 1.

Switch-lists satisfy neither bounded disjunction nor bounded conjunction. This
remains true when the functions to be disjoined (resp. conjoined) are on the
same set of variables.

###### Proof.

For disjunction, this follows directly from Observation 1 and Proposition 1
since the outcome of any polynomial-time algorithm would in particular be of
polynomial size.

For conjunction, let us define $f^{\prime}_{1}=\neg f_{1}$ and
$f^{\prime}_{2}=\neg f_{2}$. Observe that a switch-list of $f$ can be negated
in constant time by simply flipping the value $f(0)$ (keeping the same
permutation of variables). Clearly $f^{\prime}_{1}\wedge f^{\prime}_{2}=\neg
f_{1}\wedge\neg f_{2}=\neg(f_{1}\vee f_{2})$ and the lower bound for
$f_{1}\vee f_{2}$ from Proposition 1 is of course valid also for
$\neg(f_{1}\vee f_{2})$. This gives us an identical lower bound for the size
of any switch list representing $f^{\prime}_{1}\wedge f^{\prime}_{2}$. ∎

##  4 Conclusion

I was shown that switch-lists neither satisfy bounded disjunction nor bounded
conjunction. This even remains true if both inputs depend on the same set of
variables. This completes the analysis of switch-lists in the framework of the
knowledge compilation map.

Let us remark that for practical applicability of switch-lists, this is rather
bad news. Many classical approaches to practical knowledge compilation use so-
called bottom-up compilation: given a conjunction of clauses, or more
generally constraints, $F=\bigwedge_{i=1}^{m}C_{i}$, one first computes
representations $R(C_{i})$ of individual constraints $C_{i}$. Then one uses
efficient conjunction to iteratively conjoin the $R(C_{i})$ to get a
representation of $F$. Since conjunction of even two switch-lists is hard in
general, this approach is ruled out by our results.

To better understand when switch-lists are useful, it would be interesting to
find classes of functions for which small switch-list representations can be
computed efficiently, either theoretically or with heuristic approaches.

### Acknowledgement.

The author would like to thank Ondřej Čepek for helpful comments on an earlier
version of this paper.

## References

  * [1] Milos Chromý and Ondrej Čepek.  Properties of switch-list representations of boolean functions.  J. Artif. Intell. Res., 69:501–529, 2020.
  * [2] Adnan Darwiche and Pierre Marquis.  A knowledge compilation map.  J. Artif. Intell. Res., 17:229–264, 2002.
  * [3] Hélène Fargier and Pierre Marquis.  Disjunctive closures for knowledge compilation.  Artif. Intell., 216:129–162, 2014.
  * [4] Hélène Fargier, Pierre Marquis, and Alexandre Niveau.  Towards a knowledge compilation map for heterogeneous representation languages.  In Francesca Rossi, editor, IJCAI 2013, Proceedings of the 23rd International Joint Conference on Artificial Intelligence, Beijing, China, August 3-9, 2013, pages 877–883. IJCAI/AAAI, 2013.
  * [5] Knot Pipatsrisawat and Adnan Darwiche.  New compilation languages based on structured decomposability.  In Dieter Fox and Carla P. Gomes, editors, Proceedings of the Twenty-Third AAAI Conference on Artificial Intelligence, AAAI 2008, Chicago, Illinois, USA, July 13-17, 2008, pages 517–522. AAAI Press, 2008\.
  * [6] Baruch Schieber, Daniel Geist, and Ayal Zaks.  Computing the minimum DNF representation of boolean functions defined by intervals.  Discret. Appl. Math., 149(1-3):154–173, 2005.
  * [7] Ondřej Čepek and Radek Hušek.  Recognition of tractable dnfs representable by a constant number of intervals.  Discret. Optim., 23:1–19, 2017.