# Framing Algorithmic Recourse for Anomaly Detection

Debanjan Datta<EMAIL_ADDRESS>Virginia TechArlingtonVirginiaUnited States ,
Feng Chen<EMAIL_ADDRESS>University of Texas, DallasDallasTexasUnited
States and  Naren Ramakrishnan<EMAIL_ADDRESS>Virginia
TechArlingtonVirginiaUnited States

(2022)

###### Abstract.

The problem of algorithmic recourse has been explored for supervised machine
learning models, to provide more interpretable, transparent and robust
outcomes from decision support systems. An unexplored area is that of
algorithmic recourse for anomaly detection, specifically for tabular data with
only discrete feature values. Here the problem is to present a set of
counterfactuals that are deemed normal by the underlying anomaly detection
model so that applications can utilize this information for explanation
purposes or to recommend countermeasures. We present an approach—Context
preserving Algorithmic Recourse for Anomalies in Tabular data (CARAT), that is
effective, scalable, and agnostic to the underlying anomaly detection model.
CARAT uses a transformer based encoder-decoder model to explain an anomaly by
finding features with low likelihood. Subsequently semantically coherent
counterfactuals are generated by modifying the highlighted features, using the
overall context of features in the anomalous instance(s). Extensive
experiments help demonstrate the efficacy of CARAT.

Anomaly Detection; Algorithmic Recourse; Deep Learning

††copyright: acmcopyright††conference: Proceedings of the 28th ACM SIGKDD
Conference on Knowledge Discovery and Data Mining; August 14–18, 2022;
Washington DC, DC, USA††booktitle: Proceedings of the 28th ACM SIGKDD
Conference on Knowledge Discovery and Data Mining (KDD ’22), August 14–18,
2022, Washington DC, DC, USA ††journalyear: 2022††doi: 10.1145/3534678.3539344

##  1\. Introduction

Algorithmic recourse can be defined as a a set of actions or changes that can
change the outcome for a data instance with respect to a machine learning
model, typically from an unfavorable outcome to a favorable one (Joshi et al.,
2019). This is an important and challenging task with practical applicability
in domains such as healthcare, hiring, insurance, and commerce that
incorporate machine learning models into decision support systems (Prosperi et
al., 2020; Karimi et al., 2020b). Algorithmic recourse is closely related to
explainability, specifically counterfactual explanations that are important to
improve fairness, transparency, and trust in output of machine learning (ML)
models. Indeed the most cited and intuitive explanation of algorithmic
recourse presents an example how to change input features of bank loan
application decided by a black-box ML algorithm to obtain a favorable outcome
(Karimi et al., 2021).

Although the primary focus of algorithmic recourse has been in supervised
learning contexts (Mothilal et al., 2020), specifically classification based
scenarios, it is also applicable in other scenarios. In this work, we address
the research of how to frame algorithmic recourse for outcomes of unsupervised
anomaly detection. Specifically, we seek to obtain a set of actions to modify
the feature values of a data instance deemed anomalous by a black-box anomaly
detection model such that it is no longer anomalous. A motivating example
would be the case of a shipment transaction that is flagged as suspicious or
illegal by a monitoring system employing anomaly detection, and our exploring
what needs to be modified in this transaction to no longer merit that outcome.
An entity such as a trading company might seek to address its future shipment
patterns, by adjusting routes, products or suppliers to avoid getting flagged
as potentially fraudulent – thus motivating the problem of algorithmic
recourse for anomaly detection.

Algorithmic recourse for anomaly detection has some factors that
differentiates it from the classification based scenario due to the underlying
ML model w.r.t. which one tries to achieve a different but favorable outcome.
While classification models are supervised, anomaly detection models are
mostly unsupervised, and archetypes of anomalies are difficult to determine
and are application scenario dependant. Prior works consider tabular data for
algorithmic recourse in the context of classification (Karimi et al., 2020b;
Rawal and Lakkaraju, 2020), and use comparatively simpler datasets where
features are mostly real-valued. We explore the scenario where features are
strictly categorical with high dimensionality (cardinality), such as found in
real world data from commerce, communication and shipping (Cao et al., 2018;
Datta et al., 2020). Concepts such as proximity in the context of
counterfactuals are simpler to define for real-valued data. Moreover, metrics
used in classification specific algorithmic recourse do not directly translate
to the scenario of anomaly detection. Our key contributions in this work are:

  1. (i)

A novel formulation for the unexplored problem of algorithmic recourse for
unsupervised anomaly detection.

  2. (ii)

A novel approach CARAT to generate counterfactuals for anomalies in tabular
data with categorical features. CARAT is demonstrated to be effective,
scalable and agnostic to the underlying anomaly detection model.

  3. (iii)

A new set of metrics that can effectively quantify the quality of the
generated counterfactuals w.r.t. multiple objectives.

  4. (iv)

Empirical results on multiple real world shipment datasets along with a case
study highlighting the practical utility of our approach.

##  2\. Preliminaries

Tabular data with strictly categorical attributes can be formally represented
in terms of domains and entities (Datta et al., 2020). A domain or attribute
or categorical feature is defined as a set of elements sharing a common
property, e.g. Port. A domain consists of a set of entities which are the set
of possible values for the categorical variable, e.g. Port: { Baltimore, New
York, …}. Context (Datta et al., 2020) is defined as the reference group of
entities with which an entity occurs, implying an entity can be present in
multiple contexts. A data instance (record) is anomalous if it contains
unexpected co-occurrence among two or more of its entities (Das and Schneider,
2007; Hu et al., 2016; Datta et al., 2020).

###### Definition 0 (Anomalous Record).

An anomalous record is a record where certain domains have entity values that
are not consistent with the remaining entity values, termed as the context,
with respect to the expected data distribution.

Explanation for a model typically refers to an attempt to convey the internal
state or logic of an algorithm that leads to a decision (Wachter et al.,
2017). Closely related to the idea of explanations are counterfactuals.
Counterfactuals are hypothetical examples that demonstrate to an user how a
different and desired prediction can be obtained. Algorithmic recourse has
been defined as an actionable set of changes that can be made to change the
prediction of a system with respect to a data instance from an unfavourable
one to a desirable one (Joshi et al., 2019). The idea is to change one or more
of the feature values of the input in an feasible manner in order to produce a
favorable outcome. Algorithmic recourse has been explored in the context of
mostly classification problems, with a generalized binary outcome scenario.
Algorithmic recourse for anomaly detection is an important yet mostly
unexplored problem. In this work, the hypothetical instances that are the
result of algorithmic recourse on a data instance are referred to as
counterfactuals or recourse candidates. It is important to note that while
counterfactual explanations provide explanations through contrasting examples,
algorithmic recourse refers to the set of actions that provides the desired
outcome.

While nominal points are assumed to be generated from an underlying data
distribution $\mathcal{D}$, anomalies can be assumed to be generated from a
different distribution $\mathcal{D^{\prime}}$. It can be hypothesized that an
anomaly $\textbf{x}_{a}\sim\mathcal{D^{\prime}}$, is generated from some
$\textbf{x}_{n}\sim\mathcal{D}$ through some transformation function set
$\mathcal{F}$, such that
$\textbf{x}_{a}=\mathcal{F}(\textbf{x}_{n})\sim\mathcal{D^{\prime}}$. A
simplifying view of $\mathcal{F}$ can be a process of feature value
perturbation or corruption. Therefore, we can also hypothesize that there
exists some arbitrary function set $\mathcal{G}$, such that
$\mathcal{G}(\textbf{x}_{a})\sim\mathcal{D}$ and possibly
$\mathcal{G}(\textbf{x}_{a})\neq\textbf{x}_{n}$ — which is emulated through
algorithmic recourse.

###### Definition 0 (Algorithmic Recourse for Anomaly Detection).

Algorithmic recourse for anomaly detection can be defined as a set of
actionable changes on an anomalous data instance, such that it is no longer
considered an anomaly with respect to the underlying anomaly detection model.

Specifically, we consider the research question that given a row of tabular
data, with strictly categorical values, which is deemed anomalous by an
anomaly detection model $\mathcal{M}_{AD}$ — how can we generate a set of
hypothetical records such which would be deemed normal by $\mathcal{M}_{AD}$.
In this setting, without loss of generality we consider $\mathcal{M}_{AD}$ to
be  (i) trained using a training set which is assumed to be clean (Chen et
al., 2016; Datta et al., 2020),  (ii) a likelihood based model that produces
real-valued scores,  (iii) a queriable black box model

###### Problem Description 0.

Given data instance $\textbf{x}_{a}$ which is deemed anomalous by a given
anomaly detection model $\mathcal{M}_{AD}$, the objective is to generate a set
of counterfactuals $Y_{cf}$ such that $\textbf{x}_{cf}\in Y_{cf}$ is not an
anomaly according to $\mathcal{M}_{AD}$.

Since in the case of unsupervised anomaly detection an application or dataset
specific threshold is often used which is difficult to determine, we can relax
the definition of recourse to a obtain a set of counterfactuals $Y$ such that
$\textbf{x}_{cf}\in Y_{cf}$ are ranked lower by $\mathcal{M}_{AD}$ in terms of
anomaly score. There are multiple objectives that require optimization to
obtain counterfactuals that satisfy different criterion (Karimi et al., 2020a)
such as sparsity, diversity, prolixity (Keane and Smyth, 2020), proximity to
the anomalous record, low cost to the end user, along with feasibility,
actionability and non-discriminatory nature which depend on application
scenario. Since we address a general scenario without apriori application
specific knowledge, some of these problem specific objectives such as user
specific cost or feasibility are not applicable. We consider the key criterion
such as validity, diversity and sparsity and discuss them on evaluation
metrics in Section 5.

(a) Architecture of the Encoder and the Decoder-R used in the first phase of
pretraining the Encoder to learn the representation the entities and capture
overall context of the record (Section 4.1.1).

(b) Architecture of Decoder-P, which takes both the embeddings and the
contextual representation of entities from the encoder, and predicts
likelihood of each entity in the record (Section 4.1.2)

Figure 1. Architecture of the Explainer model in CARAT, comprising of the
encoder and decoder-P, which captures likelihood of each entity in the input
record.

##  3\. Related Work

In this section, we discuss the prior literature that explores the concepts of
algorithmic recourse, counterfactuals and anomaly explanation, which are
relevant in this discourse.

Explainability in machine learning models has gained burgeoning research focus
over the last decade due to the need for building trust and achieving
transparency in decision support systems that often employ black-box models.
Post-hoc explanations through feature importance has been proposed to explain
prediction of classifiers. LIME (Ribeiro et al., 2016) presents an approach to
obtaining explanations through locally approximating a model in the
neighborhood of of a prediction of a prediction. DeepLIFT (Shrikumar et al.,
2017) proposed a method to decompose the prediction of a neural network by
recursively calculating the contributions by individual neurons. SHAP
(Lundberg and Lee, 2017)present a unified framework that assigns each feature
an additive importance measure for a particular prediction, based on Shapley
Values. InterpretML (Nori et al., 2019) presents an unified framework for Ml
interpretability.

There has been recent work on explaining outcome for anomaly detection models
(Antwarg et al., 2021; Macha and Akoglu, 2018; Yepmo et al., 2022). ACE (Zhang
et al., 2019) proposes an approach for explaining anomalies in cybersecurity
datasets. While some some anomaly detection methods such as LODI (Dang et al.,
2013) and LOGP (Dang et al., 2014) provide feature importance to explain
anomalies by design, most methods employ a post-hoc explanation approach.
DIFFI (Carletti et al., 2020) provides explanations for outputs from an
Isolation Forest. Explanation through gradient based approaches have been
proposed in anomaly detection based methods on neural networks (Amarasinghe et
al., 2018; Nguyen et al., 2019; Kauffmann et al., 2020).

Karimi et al. (Karimi et al., 2020b) presents a comprehensive survey on
algorithmic recourse. Ustun et al. (Ustun et al., 2019) introduced the notion
of actionable recourse, that ensures that the counterfactuals are obtained
through appropriate feature value modification. DiCE (Mothilal et al., 2020)
presents an a framework for generating and evaluating a diverse set of
counterfactuals based on determinantal point processes. Neural network model
based approaches for generating counterfactuals have also been proposed
(Pawelczyk et al., 2020; Mahajan et al., 2019; Chapman-Rounds et al., 2021).
Causal reasoning has also been explored towards algorithmic recourse (Karimi
et al., 2020a; Prosperi et al., 2020; Karimi et al., 2021; Crupi et al.,
2021). Approaches for recourse based on heuristics, specifically using genetic
algorithms have also been explored (Sharma et al., 2019; Barredo-Arrieta and
Del Ser, 2020; Dandl et al., 2020). To our knowledge, only one method RCEAA
(Haldar et al., 2021) has been proposed towards recourse in anomaly detection
based on autoencoders with real valued inputs.

##  4\. Algorithmic Recourse Through Modeling Context

Algorithmic recourse consists of two steps:  (i) Understanding what is causing
a data instance to be an anomaly,  (ii) How to define a set of actions to
modify the feature values in order to remedy the unfavorable outcome.  We
propose CARAT: Context preserving Algorithmic Recourse for Anomalies in
Tabular Data that decomposes the task into these two sequential logical steps
and address them. CARAT comprises of a model based approach to identify the
presence of entities that causes the record to be anomalous, and an algorithm
to modify those feature values in the record for recourse.

###  4.1. Explainer Model

Given a record or data instance with categorical features, the tuple of
entities is anomalous when one or more of the entities are out-of-context with
respect to the remaining entities (Das and Schneider, 2007) with unexpected
co-occurrence patterns. We use a Transformer (Vaswani et al., 2017) based
architecture to jointly model the context of the entities of records.
Transformers have been extensively utilized in other applications on text,
image and tabular data (Huang et al., 2020). A record can be considered as a
sequence of entities, without a predefined ordering of domains or any semantic
interpretation of the relative ordering. Transformer based architectures are
appropriate for tabular data with categorical features since  (i) they can
handle large cardinality values for each category and are scalable  (ii) can
provide contextual representations of entities  (iii) can model context with a
prespecified ordering of domains  and do not consider any relative ordering
among the domains (categories). We adopt an encoder-decoder architecture,
similar to language models (Devlin et al., 2019) with the objective to predict
the likelihood of each entity in a given record with possible corruptions. The
predicted likelihood for each entity is conditioned on the context—implicitly
capturing the pair-wise and higher order co-occurrence patterns among
entities.

####  4.1.1. Pretrained Row Encoder

The encoder has a transformer based architecture and consists of multiple
layers. Sequential architectures are not effective for rows of tabular data
where the relative ordering of entities (and domains) do not have any semantic
interpretation. To handle domains with large number of entities, a $m$
parallel domain specific embedding layers are used with same dimensionality.
To inform the model which domain an entity embedding vector belongs to, we
utilize positional encoding (Devlin et al., 2019) vectors which is
concatenated to each of the entity embedding vectors. The tuple of vectors is
then passed to the subsequent transformer block comprising of multiple layers
of transformers.

To train the encoder such that it learns contextual representation for each
entity in a record, we require a corresponding decoder and training objective
which we design as follows. We refer to this decoder as decoder-R, which is
used for pretraining the encoder and not in the final objective. Decoder-R
comprises of multiple fully connected layers and is trained to reconstruct
data. In order to aid the network to retain information and reconstruct it
accurately, the contextual entity embeddings from the encoder layer are
augmented with positional vectors through concatenation, after the first fully
connected layer. Note that both the encoder and decoder-R utilize positional
encoding to indicate domain and help the model reconstruct the entity for a
given domain using the contextual embedding. The remainder of the decoder-R
consists of $m$ parallel dense layers with GELU activation, with the last
layer being softmax to obtain the index of the entity for a specific domain.
Note that while the first transformation layer is shared for entity embedding
of all domains, the latter layers are domain specific. The encoder and
decoder-R are jointly trained, using a reconstruction based objective, similar
to Masked Language Model where we randomly perturb or remove entities from
records and train the model to predict the correct one from the partial
context. The trained encoder captures the shallow embedding in the first layer
as well as the contextual representation of entities in the record.

####  4.1.2. Entity Likelihood Prediction Model

The decoder-P is designed to predict the likelihood of each entity in a record
as output—using the outputs from the pretrained encoder as it’s input. The
input to decoder-P consists of  (i) The embedding representation for the
$j^{th}$ entity $e^{j}$, obtained from the first embedding layer of the
encoder ($x_{0}^{j}\in R^{d}$)  (ii) The contextual representation of the
entity $e^{j}$, obtained from the last layer of the encoder, $z^{j}$.  Domain
specific positional encoding vector ($p^{j}\in R^{d}$) is concatenated with
$x_{0}^{j}$ to obtain $x^{j}\in R^{2d}$. We want to capture the semantic
coherence and interaction between $x^{j}$, and the contextual representation
of the entity $z^{j}$. To accomplish this we use a Bilinear layer. The output
of this Bilinear layer is fed to a dense network with multiple hidden layers,
and finally a sigmoid activation function to obtain a likelihood of whether an
entity should occur in the given record. We utilize a simple 2-layered
architecture with ReLU activation for this domain specific dense layer. Binary
cross-entropy loss is used to train decoder-P, keeping the weights of the
pretrained encoder fixed. The training of decoder-R differs from decoder-P due
to the divergent objective. We generate labelled samples from the training
data, where we perturb samples with a probability $1-\alpha$. For $\alpha$
fraction of samples, the model is given unchanged records from the training
set to enable it to recognize expected patterns and predict higher likelihood
scores for co-occurring entities. For the remaining samples, we randomly
perturb one or more of its entities and task the decoder to recognize which of
the entities have been perturbed. It is important to note that the objective
of the explainer model is not anomaly detection, but to predict the
likelihoods of individual entities in a record.

Figure 2. Finding counterfactuals through replacing a lower likelihood entity
in a record with another entity based on semantic similarity with respect to
the other entities.

###  4.2. Generating Counterfactuals

Input : Anomalous record $\textbf{x}_{a}$; Explainer Model $\mathcal{M}_{t}$,
Anomaly Detection model $\mathcal{M}_{AD}$, Set of metapaths
$\mathrm{MP}:\\{mp_{1},mp_{2}..mp_{q}\\}$, K: No. of counterfactuals

Output : Set of conterfactuals CF for $\textbf{x}_{a}$

$\\{d_{mod}\\}$ $\leftarrow$ Set of domains $d_{j}$ where $\mathcal{M}_{t}$:
$P(e_{j}^{a}\in\textbf{x}_{a})<0.5$

If $\\{d_{mod}\\}=\phi$ : $\\{d_{mod}\\}$ $\leftarrow d_{j}$ where
$min_{j}(P(e_{j}^{a}\in\textbf{x}_{a}))$ ;

//Ensure at least a single domain is modified

for  _$d_{i}\in\\{d_{mod}\\}$ _ do

$C_{i}=\\{\\}$;

//Candidate entities to replace $e_{i}^{a}$

for  _$mp_{j}\in\mathrm{MP}$_ do

$Q=\\{\\}$

if  _$d_{i}\in mp_{j}$ _ then

$S\leftarrow Nbr(mp_{j},d_{i})$, such that $S\cap\\{d_{mod}\\}=\phi$

$Q\leftarrow Q\cup S$

for  _each $d_{q}\in Q$_ do

$P\leftarrow K\,Nearest\,Neighbors$ for $e_{q}^{a}\in d_{i}$ ;

//Find entities for $d_{i}$ similar to other entities in $\textbf{x}_{a}$

$C_{i}=C_{i}\cup P$ ;

//Update candidate entities

$CF\leftarrow\\{\\}$

for  _combinations of ( $C_{i},d_{i}\in\\{d_{mod}\\}$) for $\forall
d_{i}\in\\{d_{mod}\\}$_ do

$\textbf{x}_{cf}\leftarrow$ Replace $e_{i}^{a}\in d_{i}$ in $\textbf{x}_{a}$
with $c_{i}\in C_{i}$ ;

$CF$ $\leftarrow$ $CF$ $\cup$ $\textbf{x}_{cf}$

Score $\textbf{x}_{cf}\in CF$ using $\mathcal{M}_{AD}$ ;

$CF$ $\leftarrow$ Least K anomalous records in $CF$ ;

Algorithm 1 Counterfactual generation in CARAT

Records in tabular data with categorical features can be considered as tuple
of entities, with data specific inherent relationships between the domains
(attributes). For instance in the case of shipment records, products being
shipped are closely related to the company trading them and their origin. Many
real-life applications involve tabular data which can be represented as a
heterogeneous graph or Heterogeneous Information Network (HIN). A HIN (Sun et
al., 2011) is formally defined as a graph $\mathrm{G}=(V,E)$ with a object
type mapping function $\phi:V\rightarrow\mathcal{A}$ and edge type mapping
function $\psi:E\rightarrow\mathcal{R}$. Here $v\in V$ are the nodes
representing entities, $\phi(v)\in\mathcal{A}$ are the domains, $e\in E$ are
the edges representing co-occurrence between entities and
$\psi(e)\in\mathcal{R}$. A metapath (Sun et al., 2011) or metapath schema is
an abstract path defined on the graph network schema of a heterogeneous
network that describes a composite relation $\mathbf{R}=R_{1}oR_{2}\ldots
R_{l}$ between nodes of type $A_{1},A_{2}\ldots A_{l+1}$, capturing
relationships between entities of different domains. There can exist multiple
metapaths, and we consider $\mathcal{R}$ and thus the metapaths to be
symmetric in our problem setting. Metapaths have been utilized in similarity
search in complex data, and to find patterns through capturing relevant entity
relationships (Cao et al., 2018). Recent approaches on knowledge graph
embeddings (KGE) (Wang et al., 2017) have demonstrated their effectiveness in
capturing the semantic relationships between objects of different types in
knowledge graphs which are HINs. Many approaches for KGE consider symmetric
relationships as in our case, and is more generally applicable. We choose one
such model DistMult (Yang et al., 2014) to obtain KGE for the entities in our
data. DistMult uses both node and edge embeddings to predict semantically
similar nodes, since it models relationships between entities in form of
$<e_{a},R,e_{b}>$.

In generating counterfactuals for an anomalous record, we intend to replace
the entity $e_{j}$(or entities) which is predicted to have low likelihood by
the explainer model, given the context comprising of the other entities in the
record. The intuition is to replace such entities with other entities (of the
corresponding domain) which are semantically similar to the other entities in
the record. Let $\textbf{x}_{a}$ be the anomalous record and let entity
$e_{j}^{p}$ in domain $d_{j}$ be selected for replacement. In this task, we
utilize the associated HIN constructed from the data, along with the set of
metapaths $MP=\\{mp_{1},mp_{2}\ldots mp_{q}\\}$ that are defined using domain
knowledge. Here metapath $mp_{i}$ is of the form $\\{d_{a},d_{b}\ldots
d_{i}\\}$. Thus, candidates to replace $e_{j}^{p}$ are selected using the
metapaths that contain $d_{j}$. Let us consider one such metapath $mp_{p}$
such that $d_{j}\in mp_{p}$, with relations of $(d_{i},d_{j})$ and
$(d_{j},d_{k})$. Let the respective entities in $\textbf{x}_{a}$ for $d_{i}$
and $d_{k}$ be $e_{i}^{a}$ and $e_{k}^{a}$. In a generated counterfactual
$\textbf{x}_{cf}$, the entity $e_{j}^{cf}$ that replaces $e_{j}^{a}$ should
ideally be semantically similar to $e_{i}^{a}$ and $e_{k}^{a}$. KGE can be
effectively used for this task. This idea is described in Figure 2. We find
$K$ nearest entities to $e_{i}^{a}$ and $e_{k}^{a}$, belonging to domain
$d_{j}$. Note that it is possible that $d_{i}$ or $d_{k}$ is null based on the
schema of $mp_{p}$. We replace the entities in the domains with low likelihood
with all combinations of the candidate replacements for the respective domains
to obtain the set of candidate counterfactuals, of which $K$ least anomalous
are chosen. The steps are summarized in Algorithm 1.

##  5\. Evaluation Metrics

Evaluation metrics are crucial to understanding the performance of
counterfactual generation methods, more so due to the fact that generated
counterfactuals have multiple objectives and associated trade-offs. We discuss
some of the metrics proposed in prior literature, and their limitations in the
current problem setting. Further, we propose a set of new metrics that are
more appropriate.

###  5.1. Existing Metrics for Counterfactuals

Recourse Correctness or validity (Mothilal et al., 2020) captures the ratio of
counterfactuals that are accurate in terms obtaining the desired outcome from
the blackbox prediction model. For unsupervised anomaly detection since
$\mathcal{M}_{AD}$ provides a real valued likelihood (or anomaly score), a
direct prediction (decision value) is unavailable. Recourse Coverage (Rawal
and Lakkaraju, 2020) refers to quantification of the criterion that the
algorithmic recourse provided covers as many instances as possible. Distance
(Crupi et al., 2021; Karimi et al., 2020b) or proximity measures the mean
feature-wise distance between the original data instance and the set of
recourse candidates. Distance is often calculated separately for categorical
and continuous attributes. For continuous attributes $l_{p}$ norms (Dhurandhar
et al., 2018) or their combinations are used whereas for categorical(discrete)
variables overlap measure (Chandola et al., 2007) $\mathbbm{1}(x_{i}=x_{j})$
has been used. With purely categorical attributes, this measure however fails
to convey any information other than merely how many of the attributes are
different in the counterfactual.

Cost (Crupi et al., 2021; Karimi et al., 2020a) refers to the cost incurred in
changing a particular feature value in a recourse candidate. Prior works have
utilized $l_{p}$ norms to quantify this criteria, for real-valued features. In
our problem scenario, this metric is directly not applicable without any
external real-world constraints which can help quantify the difference in cost
in changing $x_{i}$ to $x_{j}$ vs. $x_{k}$. Diversity (Mothilal et al., 2020)
refers to the feature-wise distances between the set of recourse candidates.
Diversity encourages sufficient variation among the set of recourse candidates
so that it increases the chance of finding a feasible solution. However, it
has been noted that in certain cases diversity as an objective correlates
poorly with user cost (Yadav et al., 2021). Sparsity (Mothilal et al., 2020)
refers to the number of features that are different in the recourse
candidates, with respect to the original data instance.

###  5.2. Proposed Metrics

We propose a new set of metrics based on previously defined metrics, which are
more suited to our problem setting.

Sparsity-Index: Sparsity is an important objective along with diversity that
encourages minimal change is made to a data instance in terms of features. To
capture this notion, we define Sparsity Index for tabular data with
categorical features. Let $\textbf{x}_{a}$ be the anomalous record, and
$d_{j}$ be the $j^{th}$ domain or feature.

(1) $Sparsity\,Index=\frac{1}{|Y|}\Sigma_{\textbf{x}\in
Y}\frac{1}{1+\Sigma_{d_{j}}\mathbbm{1}(x_{j}\neq x_{j}^{a})}$

The values of Sparsity Index $\in[0.5,1)$, with the low value corresponding to
modification of all feature values and the maximum value corresponding to
none.

Coherence: We define coherence as measure to quantify the consistency of the
counterfactuals similar to density consistency (Karimi et al., 2020b). Let
$\mathcal{D}_{p}$ be the set of domains which are modified in $\textbf{x}_{a}$
to obtain a counterfactual $\textbf{x}_{cf}$, $\mathcal{D}_{r}$ be the
remaining domains. Let $e_{j}$ be the entity in $\textbf{x}_{cf}$ for domain
$j$. Coherence measures the mean probability of co-occurrence of the entities
$e_{i}\in D_{p}$ with $e_{j}\in D_{r}$. Maximizing coherence implies
$e_{j}^{r}$ in $\textbf{x}_{a}$ is replaced with a candidate entity
$e_{j}^{cf}$ in $\textbf{x}_{cf}$ which has a high probability of co-
occurrence given the context of other entities of $D_{r}$ in $\textbf{x}_{a}$,
and leads to plausible counterfactuals.

(2) $coherence=\Sigma_{e_{i}\in
D_{p}}\frac{1}{|\mathcal{D}_{r}|}\Sigma_{e_{j}\in\mathcal{D}_{r}}P(e_{i},e_{j})$

Conditional Correctness: This metric quantifies the validity of the
counterfactuals, conditional upon the underlying anomaly detection model
$\mathcal{M}_{AD}$ which has a scoring function $score_{\mathcal{M}}()$. Let
$\mathcal{D}_{k}$ be a set of randomly chosen data instances from the training
and testing set. Let $\textbf{x}_{a}$ be the anomalous record and let the rank
of $\textbf{x}_{a}$ in $\textbf{x}_{a}\cap\mathcal{D}_{k}$ be $r$, sorted by
$score_{\mathcal{M}}()$ with appropriate order. Without loss of
generalization, we can assume a higher score indicates a more normal or
nominal data instance and a low score indicates anomalousness. For $x\in Y$,
where Y is the set of counterfactuals, conditional correctness can be defined
as

(3) $CC=\frac{1}{|Y|}\Sigma_{x\in Y}\mathbbm{1}(Rank(x)-r>0)$

This implies that $x$ is ranked lower in terms of being an anomaly, since
higher ranked data instances are more anomalous. Ranking is a more suitable
approach to designing a metric than utilizing thresholds which are data and
application dependant an is difficult to determine. The relative ordering of
records are important in this setting, since test instances are sorted based
on $score_{\mathcal{M}}()$.

Feature Accuracy: The concept of anomaly in tabular data with categorical
variables has been described as one or more attributes being out of context
with respect to the others, as discussed in Section 2. Therefore, it is
important to accurately measure how well can a model identify which of the
domain values should be modified in $\textbf{x}_{a}$, and relates to the
explanation aspect of algorithmic recourse. This requires having a Gold
Standard (ground truth) knowledge where we know which domain values (features)
have been corrupted and the entities for those domains are out of context. Let
$dom$ be the set of domains (features) with $m$ domains. Let $\mathbbm{q()}$
be a binary valued function that has value $1$ if the domain value is changed
in counterfactual $\textbf{x}_{cf}$ from $\textbf{x}_{a}$
$(r_{j}\in\textbf{x}_{a}\rightarrow x_{j}\in\textbf{x}_{cf})$ is an actual
cause of the anomaly or if a domain value remains unchanged if it was not a
cause of the anomaly.

(4) $FA=\frac{1}{|Y|}\Sigma_{\textbf{x}_{cf}}\frac{1}{m}\Sigma_{j\in
dom}\mathbbm{1}(\mathbbm{q}(r^{j}\rightarrow x^{j})=1)$

Heterogeneity: Although diversity is an important objective for recourse
candidates, existing diversity metrics like Count Diversity (Mothilal et al.,
2020) are inadequate. A trivial random modification of all feature values will
maximize such metrics for our setting with strictly categorical features,
where distance between discrete feature values is computed using overlap
measure. Two factors are important here:  (i) the variation among the entities
that are proposed to replace original entity in $\textbf{x}_{a}$ and  (ii) the
correct domain’s value is modified or not.  We require the Gold Standard
(ground truth) to determine whether the counterfactual modifies the a correct
domain’s value. In our experiments, use of synthetic anomalies enables
calculation of this metric. Between any two pair of recourse candidates,
heterogeneity encourages dissimilarity while taking into account if both the
pair of counterfactuals modify the correct domain’s value. Let $m$ be the
number of domains, $K$ be the size of the set of counterfactuals $Y$, and
$\mathbbm{1}(\mathbbm{q})$ be 1 if the correct domain has been modified.

(5)
$\centering\begin{multlined}H=\frac{1}{K^{2}m}\Sigma_{i=1}^{K-1}\Sigma_{j=i+1}^{K}\Sigma_{l=1}^{m}w_{ij}^{l}\mathbbm{1}(x_{i}^{l}\neq
x_{j}^{l})\\\ w_{ij}^{l}=\mathbbm{1}(\mathbbm{q}(r^{l}\rightarrow
x_{i}^{l})=1)*\mathbbm{1}(\mathbbm{q}(r^{l}\rightarrow
x_{j}^{l})=1)\end{multlined}H=\frac{1}{K^{2}m}\Sigma_{i=1}^{K-1}\Sigma_{j=i+1}^{K}\Sigma_{l=1}^{m}w_{ij}^{l}\mathbbm{1}(x_{i}^{l}\neq
x_{j}^{l})\\\ w_{ij}^{l}=\mathbbm{1}(\mathbbm{q}(r^{l}\rightarrow
x_{i}^{l})=1)*\mathbbm{1}(\mathbbm{q}(r^{l}\rightarrow
x_{j}^{l})=1)\@add@centering$

Table 1. Details of the datasets used for evaluation. Dataset | Source |  Total entity count |  Domain Count |  Train size
---|---|---|---|---
Dataset-1 | US Import | 6353 | 8 | 38291
Dataset-2 | US import | 6151 | 8 | 35177
Dataset-3 | US import | 7340 | 8 | 43495
Dataset-4 | Colombia Export | 4008 | 5 | 16758
Dataset-5 | Ecuador Export | 3198 | 7 | 13956

##  6\. Empirical Evaluation

The key objective here is to obtain counterfactuals for for anomalies in
tabular data with categorical features. We consider MEAD (Datta et al., 2020)
and APE (Chen et al., 2016) as the base anomaly detection models suited to
categorical tabular data. For our comparative evaluation against baselines we
use MEAD. For an objective and quantifiable analysis of the performance of our
approach with possible alternatives, we perform extensive experiments to
capture the varied desiderata in terms of the metrics defined in Section 5.2.
Further, we analyze the computational cost as well as the stability of the
proposed approach.

Table 2. Comparison of performance of baselines and our approach for the
adopted metrics.

(a) Feature Accuracy Dataset | Replace-m | FIMAP | RCEAA | Xformer-R | CARAT
---|---|---|---|---|---
Dataset-1 | $0.8638\pm 0.0858$ | $0.1897\pm 0.0615$ | $0.2328\pm 0.0414$ | $0.9803\pm 0.0527$ | $0.9822\pm 0.0509$
Dataset-2 | $0.8581\pm 0.0930$ | $0.1899\pm 0.0617$ | $0.2585\pm 0.0488$ | $0.9771\pm 0.0641$ | $0.9813\pm 0.0583$
Dataset-3 | $0.8591\pm 0.0876$ | $0.1898\pm 0.0615$ | $0.2381\pm 0.0580$ | $0.9769\pm 0.0633$ | $0.9786\pm 0.0630$
Dataset-4 | $0.7167\pm 0.1152$ | $0.2995\pm 0.0962$ | $0.3094\pm 0.0988$ | $0.9111\pm 0.1388$ | $0.9164\pm 0.1370$
Dataset-5 | $0.5658\pm 0.1566$ | $0.2181\pm 0.0698$ | $0.2761\pm 0.0542$ | $0.9577\pm 0.0855$ | $0.9601\pm 0.0831$

(b) Heterogenity Dataset | Replace-m | FIMAP | RCEAA | Xformer-R | CARAT
---|---|---|---|---|---
Dataset-1 | $0.3589\pm 0.3056$ | $0.9737\pm 0.0261$ | $0.6625\pm 0.1098$ | $0.9693\pm 0.0821$ | $0.9017\pm 0.1244$
Dataset-2 | $0.3728\pm 0.3273$ | $0.9716\pm 0.0278$ | $0.6646\pm 0.1499$ | $0.9579\pm 0.1311$ | $0.8930\pm 0.1522$
Dataset-3 | $0.3472\pm 0.3164$ | $0.9737\pm 0.0262$ | $0.6503\pm 0.1236$ | $0.9662\pm 0.1016$ | $0.8978\pm 0.1319$
Dataset-4 | $0.0634\pm 0.1216$ | $0.9145\pm 0.0965$ | $0.6863\pm 0.2824$ | $0.8164\pm 0.2333$ | $0.7805\pm 0.2264$
Dataset-5 | $0.1122\pm 0.2041$ | $0.9570\pm 0.0441$ | $0.6176\pm 0.1957$ | $0.8879\pm 0.1957$ | $0.8378\pm 0.2003$

(c) Coherence Dataset | Replace-m | FIMAP | RCEAA | Xformer-R | CARAT
---|---|---|---|---|---
Dataset-1 | $0.0012\pm 0.0005$ | $0.0000\pm 0.0000$ | $0.0002\pm 0.0001$ | $0.0004\pm 0.0003$ | $0.0025\pm 0.0022$
Dataset-2 | $0.0008\pm 0.0004$ | $0.0000\pm 0.0000$ | $0.0002\pm 0.0001$ | $0.0003\pm 0.0002$ | $0.0020\pm 0.0018$
Dataset-3 | $0.0013\pm 0.0006$ | $0.0000\pm 0.0000$ | $0.0002\pm 0.0001$ | $0.0004\pm 0.0003$ | $0.0030\pm 0.0028$
Dataset-4 | $0.0003\pm 0.0004$ | $0.0000\pm 0.0000$ | $0.0001\pm 0.0001$ | $0.0007\pm 0.0012$ | $0.0018\pm 0.0027$
Dataset-5 | $0.0007\pm 0.0007$ | $0.0000\pm 0.0001$ | $0.0005\pm 0.0005$ | $0.0011\pm 0.0014$ | $0.0037\pm 0.0047$

(d) Sparsity-Index Dataset | Replace-m | FIMAP | RCEAA | Xformer-R | CARAT
---|---|---|---|---|---
Dataset-1 | $0.8889\pm 0.0000$ | $0.5015\pm 0.0010$ | $0.5355\pm 0.0089$ | $0.8389\pm 0.0524$ | $0.8381\pm 0.0523$
Dataset-2 | $0.8889\pm 0.0000$ | $0.5016\pm 0.0010$ | $0.5391\pm 0.0089$ | $0.8381\pm 0.0523$ | $0.8572\pm 0.0463$
Dataset-3 | $0.8889\pm 0.0000$ | $0.5015\pm 0.0010$ | $0.5321\pm 0.0042$ | $0.8366\pm 0.0532$ | $0.8358\pm 0.0531$
Dataset-4 | $0.8750\pm 0.0000$ | $0.5049\pm 0.0022$ | $0.5445\pm 0.0052$ | $0.7804\pm 0.0705$ | $0.7825\pm 0.0678$
Dataset-5 | $0.8333\pm 0.0000$ | $0.5027\pm 0.0014$ | $0.5468\pm 0.0039$ | $0.8305\pm 0.0586$ | $0.8300\pm 0.0585$

(e) Conditional Correctness Dataset | Replace-m | FIMAP | RCEAA | Xformer-R | CARAT
---|---|---|---|---|---
Dataset-1 | $1.0000\pm 0.0000$ | $0.7404\pm 0.1388$ | $0.6230\pm 0.1609$ | $0.7389\pm 0.1947$ | $0.9774\pm 0.0631$
Dataset-2 | $1.0000\pm 0.0000$ | $0.6646\pm 0.1556$ | $0.4990\pm 0.2092$ | $0.7414\pm 0.2077$ | $0.9733\pm 0.0800$
Dataset-3 | $1.0000\pm 0.0000$ | $0.7544\pm 0.1342$ | $0.6090\pm 0.1489$ | $0.7634\pm 0.2061$ | $0.9727\pm 0.0679$
Dataset-4 | $1.0000\pm 0.0000$ | $0.3666\pm 0.2223$ | $0.3770\pm 0.2712$ | $0.5572\pm 0.2725$ | $0.8411\pm 0.2429$
Dataset-5 | $1.0000\pm 0.0000$ | $0.5122\pm 0.1867$ | $0.4240\pm 0.1744$ | $0.5858\pm 0.2708$ | $0.8695\pm 0.2023$

(f) Summary (normalized mean) of all metrics across all datasets, for the baselines and our proposed approach CARAT. Replace-m |  FIMAP | RCEAA |  Xformer-R | CARAT
---|---|---|---|---
0.6150 |  0.2410 | 0.1657 |  0.6752 | 0.9183

Table 3. Comparison of metrics for generated counterfactuals using different anomaly detection models with CARAT. | Sparsity Index | Conditional Corr. | Coherence
---|---|---|---
Dataset | APE | MEAD | APE | MEAD | APE | MEAD
Dataset-1 |  $0.8764$ $\pm 0.0362$ |  $0.8730$ $\pm 0.0391$ |  $0.7323$ $\pm 0.3012$ |  $0.7237$ $\pm 0.3283$ |  $0.0004$ $\pm 0.0004$ |  $0.0004$ $\pm 0.0004$
Dataset-2 |  $0.8870$ $\pm 0.0178$ |  $0.8695$ $\pm 0.0404$ |  $0.6624$ $\pm 0.3003$ |  $0.6853$ $\pm 0.3473$ |  $0.0003$ $\pm 0.0003$ |  $0.0004$ $\pm 0.0004$
Dataset-3 |  $0.8810$ $\pm 0.0370$ |  $0.8753$ $\pm 0.0416$ |  $0.7517$ $\pm 0.2462$ |  $0.6326$ $\pm 0.3721$ |  $0.0002$ $\pm 0.0002$ |  $0.0003$ $\pm 0.0003$
Dataset-4 |  $0.8415$ $\pm 0.0058$ |  $0.8194$ $\pm 0.0517$ |  $0.6422$ $\pm 0.2703$ |  $0.5414$ $\pm 0.3877$ |  $0.0004$ $\pm 0.0003$ |  $0.0003$ $\pm 0.0004$
Dataset-5 |  $0.8720$ $\pm 0.0238$ |  $0.8654$ $\pm 0.0355$ |  $0.7694$ $\pm 0.2150$ |  $0.4601$ $\pm 0.3657$ |  $0.0009$ $\pm 0.0008$ |  $0.0015$ $\pm 0.0014$

###  6.1. Datasets

The datasets used in for the experimental and evaluation setup are real world
proprietary datasets of shipping records obtained from Panjiva Inc (Panjiva,
2019). Specifically we use 5 datasets with no overlap, constructed from a
larger corpus of records. We consider records of a time period as the training
set and subsequent time as the test set. We fix the test set size to 5000. The
dataset details are described in Table 1. The training set of each dataset is
used to train the AD model as well as the KGE model. Since we do not have
ground truth data of anomalies, we use synthetic anomalies generated from the
test set of each dataset following prior work (Chen et al., 2016), and allows
us to analyze the results using the ground truth knowledge of what caused the
record to be an anomaly.

###  6.2. Competing Baseline Methods

The area of algorithmic recourse specifically for anomaly detection is
unexplored, and to the best of our knowledge only one prior work RCEAA (Haldar
et al., 2021) exists on this. Prior work on algorithmic recourse deals with
classification scenarios and they are not directly applicable to our setting.
Moreover, as previously noted, most prior work deals with real valued or mixed
valued data where the cardinality of categorical variables are significantly
lower. Also, most prior approaches convert discrete variables to binary
vectors through one-hot encoding and treat them as real valued vectors. We
choose the following baselines for comparison:

Replace-m: This approach generates an initial candidate set of all possible
records by replacing the entity values in $m$ domains simultaneously, using
all possible combinations. The records in this candidate set are scored by the
given anomaly detection model, and $K$ top scored (least anomalous) records
are considered as set of counterfactuals. We set $m=1$ due to computational
limitations.

FIMAP (Chapman-Rounds et al., 2021): FIMAP is a model based approach for
generating counterfactuals through adverserial perturbations using a
perturbation network, for a classification setting with known labels. To train
the proxy classifier, a set of synthetic anomalies (assigned label $y=0$) and
normal instances (assigned label $y=1$). The perturbation network, which
generates counterfactuals, is trained by providing synthetic anomalies and
passing the perturbed data instance to the pretrained proxy classifier, to
obtain the desired label ($y=1$).

RCEAA (Haldar et al., 2021): RCEAA uses an optimization based objective to
exactly calculate a set of $K$ counterfactuals. Since the optimization
requires real valued inputs, we adopt real-value relaxation on the one-hot
encoded discrete representation of $\textbf{x}_{a}$ and use soft-threshold
approach to obtain discrete outputs.

Xformer-R: This method utilizes the explainer model to identify entities in an
anomaly with low likelihood. Counterfactuals are generated by replacing the
entity in the identified domains $\textbf{x}_{a}$ with entity values that are
sampled uniformly from the domain.

###  6.3. Results

We present the results for the metrics discussed in Section 5.2. For
conditional correctness, we sample a set of records containing both known
synthetic anomalies and normal data instances. It is important to note that no
single metric quantifies the different desiderata of the generated
counterfactuals. Beginning with feature accuracy, which is reported in Table
2(a) we see the approaches based on Transformer based explainer (Xformer-R and
CARAT) have significantly better performance compared to the others. This
demonstrates that the explainer can effectively identify entities in records
which do not conform to expected co-occurrence patterns. For heterogeneity,
the results are presented in Table 2(b), while CARAT performs well, but FIMAP
has somewhat better performance. This can be explained by the fact that
counterfactuals generated by FIMAP modify most of the entity values—which
violates the sparsity objective. Next we consider coherence, which that
captures how semantically similar the replaced entities are to the remaining
ones in the counterfactuals generated. As reported in Table 2(c), we find
CARAT performs significantly better. This implies the generated
counterfactuals are consistent with the underlying data distribution.
Considering sparsity, FIMAP and RECEAA have significantly lower values as
shown in Table 2(d), since the counterfactuals have multiple feature values
modified from the given anomaly. Replace-m has a high sparsity since a single
entity is modified (m=1). Xformer-R and CARAT have similar performance in
terms of sparsity. Lastly, for conditional correctness reported in Table 2(e)
we find Replace-m has perfect score due to performing exhaustive search for
least anomalous records. CARAT shows competitive performance here, better than
other baselines. Since no single metric comprehensively captures the requisite
objectives that we are trying to maximize in generating counterfactuals—we
summarize the model performances across all the metrics. For each approach, we
first obtain the average of the values across all datasets and then normalize
them. Then we perform an unweighted average across all of these normalized
metrics values to find a single performance value. The results are reported in
Table 2(f), which shows CARAT has a significant overall advantage.

Figure 3. Computational time of compared approaches. Our proposed approach
CARAT shows a clear advantage.

###  6.4. Stability

The process of algorithmic recourse is inherently dependant on the underlying
anomaly detection model which finds the anomalous data instances. Thus it is
important to understand the variation in performance of our proposed approach
in the context of the underlying AD model $\mathcal{M}_{AD}$. Here
$\mathcal{M}_{AD}$ is a black-box model, and assumed to perfectly capture the
underlying data distribution to find data instances that are true anomalies.
We choose two embedding based algorithms for tabular data with strictly
categorical features as $\mathcal{M}_{AD}$, APE (Chen et al., 2016) and MEAD
(Datta et al., 2020). We apply APE and MEAD on test sets of each of the
datasets, consider $1\%$ of the lowest scored records as anomalies, and apply
CARAT with the respective AD models on the corresponding anomalies. Comparison
of the applicable metrics defined in Section 5.2 for the generated
counterfactuals across multiple metrics are reported in Table 3. We observe
that similar performance is obtained in both cases, demonstrating the
stability of our proposed approach.

###  6.5. Computational Cost

One of the major challenges of generating recourse is the computational
complexity. We consider the computational cost in the execution phase, after
all applicable pretraining and set up has been completed. An approach with
with exponential computational complexity would be simply infeasible in most
practical scenarios. For instance the Replace-m approach outlined in Section
6.2, the complexity is $O(\Sigma_{1}^{l}\prod d_{i},d_{i+1}..d_{i+m})$, where
$m$ domains are chosen to be modified at most, $d_{i},d_{i+1}\dots d_{i+m}$
are the cardinalities of the $m$ domains, and $l={|d|\choose m}$. RCEAA
(Haldar et al., 2021) also suffers from high computational cost, as shown in
Figure 3. The major bottlenecks in such an optimization based approach are:
(i) expensive loss function  (ii) grid search for hyperparameters  (iii)
operations on a very high dimensional vector due to one-hot encoding.  Our
proposed approach is computationally efficient, since the explanation phase
uses only a pretrained model with linear time complexity and the
counterfactual generation phase that requires finding nearest neighbors is
sped up using an indexing library (Johnson et al., 2019).

Figure 4. Example of an detected anomaly (first row) and the generated
counterfactuals using CARAT.

###  6.6. Case Study of an Anomaly

We perform a short case study based on one of the anomalies detected from test
set of Dataset-1. The detected anomaly and a set of counterfactuals are shown
in Figure 4. Our explainer model finds that the entity in domain HS Code in
the anomalous record $\textbf{x}_{a}$ has low likelihood of occurrence in its
context. This is supported by the empirical data, as we find goods represented
by the entity $<$HS Code:7211$>$ was previously traded by the neither the
consignee (buyer) nor the shipper. Additionally, $<$HS Code:7211$>$ was not
previously transported through the ports of lading and unlading in the record.
Our explainer only presents a low score for this entity, and not the others in
the record although their context, which contains $<$HS Code:7211$>$, is
altered. This demonstrates that the model can accurately predict likelihood
from partially correct contextual information—which is the case for anomalous
records. Let us look at the counterfactuals, where $<$HS Code:7211$>$ is
modified to alternate values. $<$HS Code:7211$>$ refers to products of iron or
non-alloy steel. Firstly $<$HS Code:7210$>$, $<$HS Code:7225$>$ and $<$HS
Code:7219$>$ refer to products very similar to $<$HS Code:7211$>$,
specifically rolled or non-rolled steel or alloy products. $<$HS Code: 2818$>$
and $<$HS Code:7304$>$ are metal products as well. So the counterfactuals are
essentially suggesting the buyer to obtain alternate products from the
supplier (shipper). This can be explained based on the data, the particular
type of goods $<$HS Code:7211$>$ are generally not sourced from the given
origin and the shipper evidenced by any similar prior occurrences. A practical
explanation might relate to industrial production and proficiency patterns,
since certain regions specialize in specific products. From prior records it
is further observed that the consignee buys metal and alloy goods, such as
with HS Codes 7210, 7323, 7304, 7225 and 7310. Thus CARAT provides meaningful
counterfactuals for detected anomalies.

##  7\. Conclusion and Future Work

In this work we address the previously unexplored problem of algorithmic
recourse for anomaly detection in tabular data with categorical features with
high cardinality. We propose a novel deep learning based approach CARAT to
find counterfactuals for detected anomalies. We also define a set of relevant
metrics that can enable effective evaluation of counterfactuals in such a
problem setting. The scalability and efficacy of our model in terms of the
applicable metrics as well computational cost is demonstrated through
extensive experiments. However the current research leads to further research
questions. While we consider multiple objectives to optimize in the process of
algorithmic recourse there are application scenario specific constraints that
are not considered. One of the questions that has practical implications and
requires domain knowledge is the actual cost incurred by the user. Another
aspect is feasibility and actionability of counterfactuals, which are often
dependant on extrinsic factors that need to explicitly considered and
incorporated into the algorithmic recourse for anomalies. Thus there are
multiple continuing research directions which are a natural progression of the
problem we address in this work.

###### Acknowledgements.

This work was supported in part by US NSF grants CCF-1918770, NRT DGE-1545362,
and OAC-1835660 to NR, and IIS-1954376 and IIS-1815696 to FC.

## References

  * (1)
  * Amarasinghe et al. (2018) Kasun Amarasinghe, Kevin Kenney, and Milos Manic. 2018\.  Toward explainable deep neural network based anomaly detection. In _2018 11th ICHSI_. IEEE.
  * Antwarg et al. (2021) Liat Antwarg, Ronnie Mindlin Miller, Bracha Shapira, and Lior Rokach. 2021.  Explaining anomalies detected by autoencoders using Shapley Additive Explanations.  _Expert Systems with Applications_ 186 (2021), 115736.
  * Barredo-Arrieta and Del Ser (2020) Alejandro Barredo-Arrieta and Javier Del Ser. 2020.  Plausible counterfactuals: Auditing deep learning classifiers with realistic adversarial examples. In _2020 IJCNN_. IEEE.
  * Cao et al. (2018) Bokai Cao et al. 2018\.  Collective fraud detection capturing inter-transaction dependency. In _KDD 2017 Workshop on Anomaly Detection in Finance_. PMLR.
  * Carletti et al. (2020) Mattia Carletti, Matteo Terzi, and Gian Antonio Susto. 2020\.  Interpretable Anomaly Detection with DIFFI: Depth-based Isolation Forest Feature Importance.  _arXiv preprint arXiv:2007.11117_ (2020).
  * Chandola et al. (2007) Varun Chandola, Shyam Boriah, and Vipin Kumar. 2007\.  Similarity Measures for Categorical Data–A Comparative Study.  (2007).
  * Chapman-Rounds et al. (2021) Matt Chapman-Rounds et al. 2021\.  FIMAP: Feature Importance by Minimal Adversarial Perturbation. In _AAAI_ , Vol. 35.
  * Chen et al. (2016) Ting Chen et al. 2016\.  Entity embedding-based anomaly detection for heterogeneous categorical events. In _IJCAI_.
  * Crupi et al. (2021) Riccardo Crupi et al. 2021\.  Counterfactual Explanations as Interventions in Latent Space.  _arXiv e-prints_ (2021), arXiv–2106.
  * Dandl et al. (2020) Susanne Dandl et al. 2020\.  Multi-objective counterfactual explanations. In _International Conference on Parallel Problem Solving from Nature_. Springer, 448–469.
  * Dang et al. (2013) Xuan Hong Dang et al. 2013\.  Local outlier detection with interpretation. In _ECML PKDD_. Springer, 304–320.
  * Dang et al. (2014) Xuan Hong Dang et al. 2014\.  Discriminative features for identifying and interpreting outliers. In _IEEE ICDE 2014_. IEEE, 88–99.
  * Das and Schneider (2007) Kaustav Das and Jeff Schneider. 2007.  Detecting anomalous records in categorical datasets. In _13th ACM SIGKDD_. 220–229.
  * Datta et al. (2020) Debanjan Datta et al. 2020\.  Detecting Suspicious Timber Trades. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 34. 13248–13254.
  * Devlin et al. (2019) Jacob Devlin et al. 2019\.  BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding. In _NAACL_. ACL.
  * Dhurandhar et al. (2018) Amit Dhurandhar et al. 2018\.  Explanations based on the missing: Towards contrastive explanations with pertinent negatives.  _NeurIPS_ 31 (2018).
  * Haldar et al. (2021) Swastik Haldar et al. 2021\.  Reliable Counterfactual Explanations for Autoencoder Based Anomalies. In _8th ACM IKDD CODS and 26th COMAD_.
  * Hu et al. (2016) Renjun Hu, Charu C Aggarwal, Shuai Ma, and Jinpeng Huai. 2016\.  An embedding approach to anomaly detection. In _ICDE_.
  * Huang et al. (2020) Xin Huang et al. 2020\.  Tabtransformer: Tabular data modeling using contextual embeddings.  _arXiv preprint arXiv:2012.06678_ (2020).
  * Johnson et al. (2019) Jeff Johnson, Matthijs Douze, and Hervé Jégou. 2019\.  Billion-scale similarity search with gpus.  _IEEE Transactions on Big Data_ 7, 3 (2019), 535–547.
  * Joshi et al. (2019) Shalmali Joshi et al. 2019\.  Towards realistic individual recourse and actionable explanations in black-box decision making systems.  _arXiv preprint arXiv:1907.09615_ (2019).
  * Karimi et al. (2020a) Amir-Hossein Karimi et al. 2020a.  Model-agnostic counterfactual explanations for consequential decisions. In _AISTATS_. PMLR, 895–905.
  * Karimi et al. (2020b) Amir-Hossein Karimi et al. 2020b.  A survey of algorithmic recourse: definitions, formulations, solutions, and prospects.  _CoRR_ abs/2010.04050 (2020).
  * Karimi et al. (2021) Amir-Hossein Karimi, Bernhard Schölkopf, and Isabel Valera. 2021.  Algorithmic recourse: from counterfactual explanations to interventions. In _ACM FaccT_.
  * Kauffmann et al. (2020) Jacob Kauffmann et al. 2020\.  Towards explaining anomalies: a deep Taylor decomposition of one-class models.  _Pattern Recognition_ 101 (2020), 107198.
  * Keane and Smyth (2020) Mark T Keane and Barry Smyth. 2020.  Good counterfactuals and where to find them: A case-based technique for generating counterfactuals for explainable AI (XAI). In _International Conference on Case-Based Reasoning_. Springer, 163–178.
  * Lundberg and Lee (2017) Scott M Lundberg and Su-In Lee. 2017.  A Unified Approach to Interpreting Model Predictions. In _NeurIPS 2017_ , Vol. 30.
  * Macha and Akoglu (2018) Meghanath Macha and Leman Akoglu. 2018.  Explaining anomalies in groups with characterizing subspace rules.  _DMKD 2018_ 32 (2018).
  * Mahajan et al. (2019) Divyat Mahajan et al. 2019\.  Preserving causal constraints in counterfactual explanations for machine learning classifiers.  _arXiv preprint arXiv:1912.03277_ (2019).
  * Mothilal et al. (2020) Ramaravind K Mothilal et al. 2020\.  Explaining machine learning classifiers through diverse counterfactual explanations. In _ACM FAT 2020_.
  * Nguyen et al. (2019) Quoc Phong Nguyen et al. 2019\.  Gee: A gradient-based explainable variational autoencoder for network anomaly detection. In _2019 IEEE CNS_. 91–99.
  * Nori et al. (2019) Harsha Nori et al. 2019\.  Interpretml: A unified framework for machine learning interpretability.  _arXiv preprint arXiv:1909.09223_ (2019).
  * Panjiva (2019) Panjiva. 2019.  Panjiva Trade Data.  https://panjiva.com.
  * Pawelczyk et al. (2020) Martin Pawelczyk et al. 2020\.  Learning model-agnostic counterfactual explanations for tabular data. In _The Web Conference 2020_. 3126–3132.
  * Prosperi et al. (2020) Mattia Prosperi et al. 2020\.  Causal inference and counterfactual prediction in machine learning for actionable healthcare.  _NMI_ 2, 7 (2020).
  * Rawal and Lakkaraju (2020) Kaivalya Rawal and Himabindu Lakkaraju. 2020.  Beyond Individualized Recourse: Interpretable and Interactive Summaries of Actionable Recourses.  _NeurIPS_ (2020).
  * Ribeiro et al. (2016) Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. 2016\.  ” Why should i trust you?” Explaining the predictions of any classifier. In _22nd ACM SIGKDD_.
  * Sharma et al. (2019) Shubham Sharma, Jette Henderson, and Joydeep Ghosh. 2019\.  Certifai: Counterfactual explanations for robustness, transparency, interpretability, and fairness of artificial intelligence models.  _arXiv preprint arXiv:1905.07857_ (2019).
  * Shrikumar et al. (2017) Avanti Shrikumar, Peyton Greenside, and Anshul Kundaje. 2017\.  Learning important features through propagating activation differences. In _International conference on machine learning_. PMLR, 3145–3153.
  * Sun et al. (2011) Yizhou Sun et al. 2011\.  Pathsim: Meta path-based top-k similarity search in heterogeneous information networks.  _Proceedings of the VLDB Endowment_ 4, 11 (2011), 992–1003.
  * Ustun et al. (2019) Berk Ustun, Alexander Spangher, and Yang Liu. 2019\.  Actionable recourse in linear classification. In _ACM FAT_. 10–19.
  * Vaswani et al. (2017) Ashish Vaswani et al. 2017\.  Attention is all you need. In _NeurIPS_. 5998–6008.
  * Wachter et al. (2017) Sandra Wachter et al. 2017\.  Counterfactual explanations without opening the black box: Automated decisions and the GDPR.  _Harv. JL & Tech._ 31 (2017), 841.
  * Wang et al. (2017) Quan Wang et al. 2017\.  Knowledge graph embedding: A survey of approaches and applications.  _IEEE TKDE_ 12 (2017).
  * Yadav et al. (2021) Prateek Yadav et al. 2021\.  Low-Cost Algorithmic Recourse for Users With Uncertain Cost Functions.  _arXiv preprint arXiv:2111.01235_ (2021).
  * Yang et al. (2014) Bishan Yang et al. 2014\.  Embedding entities and relations for learning and inference in knowledge bases.  _ICLR_.
  * Yepmo et al. (2022) Véronne Yepmo, Grégory Smits, and Olivier Pivert. 2022\.  Anomaly explanation: A review.  _Data & Knowledge Engineering_ 137 (2022), 101946.
  * Zhang et al. (2019) Xiao Zhang et al. 2019\.  ACE–an anomaly contribution explainer for cyber-security applications. In _2019 IEEE Big Data_. IEEE, 1991–2000.

##  Appendix A Dataset Background

The datasets used in the empirical evaluation are from shipping domain and are
proprietary due to security and legal reasons. We discuss some of the
attributes of this real world data and their interpretations.

HS Code or Harmonized Tariff Schedule Codes are globally standardized codes
that define what type of goods are being transported. Carrier is the
transporting entity that operates between ports. The ports of lading and
unlading are the points where the cargo is laden onto the transporting vessel
or vehicle and and unladen from it. We have received help of collaborating
domain experts who deal with shipping data to help us understand the data
characteristics and the relationships between attributes. The original data
has many attributes which contain redundant information, and we select only
meaningful attributes from the raw data. Also we remove rows with missing
values and perform standard data cleaning to obtain our datasets.

The metapaths that describe these relationships are shown in Table 4(c). These
are designed with the knowledge of the structure of supply chains that are
captured in this Bill of Lading corpus.

Table 4. Metapaths for the datasets belonging to the three data sources, viz.
US Import, Colombia Export and Ecuador Export.

(a) Schema of the metapaths describing the relationships between attributes of
the data for US Import  Shipment Origin $\leftrightarrow$ HS Code
$\leftrightarrow$Port Of Lading
---
Shipment Destination $\leftrightarrow$ HS Code $\leftrightarrow$ Port Of
Unlading
Port Of Lading $\leftrightarrow$ HS Code $\leftrightarrow$ Carrier
HS Code $\leftrightarrow$ Carrier $\leftrightarrow$Port Of Unlading
Shipper $\leftrightarrow$ Shipment Origin $\leftrightarrow$ Port Of Lading
Consignee $\leftrightarrow$ Shipment Destination$\leftrightarrow$Port Of
Unlading
Consignee $\leftrightarrow$ Carrier $\leftrightarrow$ Shipment Destination
Shipper $\leftrightarrow$ Carrier $\leftrightarrow$ Shipment Origin
Consignee $\leftrightarrow$ Carrier $\leftrightarrow$ Port Of Unlading
Shipper $\leftrightarrow$ Carrier $\leftrightarrow$ Port Of Lading

(b) Schema of the metapaths describing the relationships between attributes of
the data for exports from Colombia.  Shipper $\leftrightarrow$ Shipment Origin
$\leftrightarrow$ HS Code
---
Consignee $\leftrightarrow$ Shipment Destination $\leftrightarrow$ HS Code
Shipment Destination $\leftrightarrow$ HS Code $\leftrightarrow$ Shipment
Origin
Shipper $\leftrightarrow$ HS Code $\leftrightarrow$ Consignee

(c) Schema of the metapaths describing the relationships between attributes of
the data for exports from Ecuador.  Shipment Destination,Goods Shipped,Port Of
Unlading
---
Shipper $\leftrightarrow$ Goods Shipped $\leftrightarrow$ Shipment Origin
Goods Shipped $\leftrightarrow$ Carrier $\leftrightarrow$ Port Of Unlading
Consignee $\leftrightarrow$ Shipment Destination $\leftrightarrow$ Port Of
Unlading
Consignee $\leftrightarrow$ Carrier $\leftrightarrow$ Shipment Destination
Shipper $\leftrightarrow$ Carrier $\leftrightarrow$ Shipment Origin
Consignee $\leftrightarrow$ Carrier $\leftrightarrow$ Port Of Unlading

##  Appendix B Experimental Setup Details

###  B.1. Hardware and Libraries

We provide the implementation details to faithfully reproduce the results
obtained. All implementation is done in Python 3.9, and uses standard
libraries such as Numpy, Pandas and scikit-learn. For optimization and neural
network based models, PyTorch (version 1.10) is used. All data preprocessing,
training and evaluation presented in this work are performed on a 40-core
machine, with a single GPU and distributed training required. To train our
Knowledge Graph Embedding model, we use the library StellarGraph, which
provides an implementation of DistMult.

###  B.2. Experimental Settings and Hyperparameters

####  B.2.1. Anomaly Detection Model

Anomaly detection for tabular data with strictly categorical features,
especially where the attributes have high dimensionality (cardinality) is a
challenging task. We choose Multi-relational Embedding based Anomaly Detection
(Datta et al., 2020) as the base anomaly detection model for our experiments.
MEAD uses an additive model based on shallow embedding, where the likelihood
of a record is a function of the magnitude of transformed sum of the entity
embeddings. We use an embedding size of $32$ for anomaly detection models our
experiments.

####  B.2.2. Synthetic Anomalies

Synthetic anomalies are generated using the approach followed in prior works
(Chen et al., 2016; Datta et al., 2020). For each record, randomly one or more
feature values are perturbed i.e. replaced with a random but valid feature
value for the categorical attribute. Since our data has at most $8$
categorical attribute we limit the number of perturbations to 2. In generating
counterfactuals we consider a balanced mix for all cases.

####  B.2.3. CARAT Explainer Model Details

For the explainer in CARAT presented in Section 4 we an entity embedding
dimension of 64. The encoder employs 4 layers of transformer blocks with 8
heads for multi-headed self-attention. The fully connected layers Decoder-R
has 3 layers, with $256$,$128$ and $64$. The fully connected layers Decoder-P
has 3 layers, with $32$ and $16$. We use the same architecture across all
datasets.

In pretraining the encoder with decoder-R, the training objective is similar
to Masked Language Model but not identical. Specifically we replace
approximately $20\%$ of entities in each records to mask, and $20\%$ of
entities are perturbed by replacement with a randomly sample entity from the
same domain. In the second phase of training the decoder-P, $\alpha$ — the
fraction of records which are not changed is set to $0.3$.

Both the pre-training phase of the encoder and the final explainer
architecture are trained for 250 epochs with a batch size of 512 and learning
rate of $0.0005$. All optimization for our model and the baselines are
performed using Adam.

####  B.2.4. CARAT KGE Details

The knowledge graph embedding model adopted here is DistMult. We use an
embedding size of $100$. The training batch size used is 1024, and we train
the model for 300 epochs. We use both the node and edge embeddings to find
entities that are similar to a target entity. Since DistMult uses
$<$head,rel,tail$>$ format to calculate similarity, we perform nearest
neighbor search for the tail entity using precomputed embedding of head nodes
and relation type.

###  B.3. Empirical Evaluation Setup

We generate counterfactuals for synthetic anomalies. For each approach, we use
a set of 400 anomalies and we generate 50 counterfactuals for each anomaly. We
use the same set of anomalies for all approaches to perform a fair comparative
evaluation. For RCEAA we are able generate counterfactuals for 40 anomalies,
due to the excessively long execution time as explained in Section 6.5.

###  B.4. Additional Detail on Competing Baselines

For the baseline models, we adapt the models to the current problem setting in
an appropriate manner. We utilize the hyperparameters provided in the original
work and do not perform significant hyperparameter tuning for these
approaches. Similarly, we do not perform significant hyperparameter tuning for
our model to tune performance as well since the objective is to demonstrate
the validity of our approach in a general setting.

####  B.4.1. RCEAA

In our implementation of RCEAA (Haldar et al., 2021), we adapt the original
approach. Firstly, we replace the autoencoder based anomaly detection model
with our likelihood based model. In the loss function, we use $10^{th}$
percentile scores of training data as the threshold so that the generated
counterfactuals have a low anomaly score (higher likelihood) according to our
anomaly model. We set $\omega$ to 2. Additionally in place of euclidean
distance, we use cosine distance since the dimensionality of the vectors is
high. In order to reduce the computational complexity due to grid search of
hyperparameters, we set the upper and lower bounds of $\lambda_{1}$ and
$\lambda_{2}$ to $0.25$ and $0.75$, and adopt step size of $0.25$. The
training epochs are increased from $15$ mentioned in the paper, to $20$. We
report the results with the lowest optimization loss.

####  B.4.2. FIMAP

For FIMAP, which is an approach for generating counterfactuals in a
classification based setting, we adapt it to our problem setting
appropriately. We adopt a more complex neural network architecture that
follows that archetype proposed in the original work. Specifically, an
embedding projection layer of dimensionality 32 is used for each categorical
variable, whose output is concatenated and fed to a fully connected network
for both proxy classifier network and perturbation network. The fully
connected network for proxy classifier has size $(256,256])$, following the
original work that uses a 3 layered neural network. The fully connected
network for perturbation model has size $(256,128,128)$ and uses dropout of
$0.2$. The value of $\tau$ in Gumbel softmax set to $0.5$. We use $10000$ data
points to train the proxy classifier for each dataset, with samples from
training set and synthetic anomalies. The batch size used is 512, and the
networks are trained for 100 epochs, with early stopping.

##  Appendix C Code and Data Link

Please find the code for this work at :

https://github.com/ddatta-DAC/Algorithmic_Recourse_AnomalyDetection

##  Appendix D Ethical Implications of Algorithmic Recourse for Anomaly
Detection

Since algorithmic recourse provides an approach towards generating
counterfactuals that are not deemed anomalous by an anomaly detection model
which may be part of a decision support system, it raises an obvious ethical
question. Is algorithmic recourse adverserial to anomaly detection i.e.
intended to help nefarious actors attempting to evade detection?

That is not our motivation here. First, data and anomaly detection models are
expected to be secured and adverserial agents would have no access to them in
order to circumvent the detection process. Recourse for anomaly detection is
intended to help the decision making process. It can enable organizations like
enforcement agencies in our case study, that make use of anomaly detection
systems to better handle false positive cases. Verified agents, whose
transactions may be erroneously flagged, could be intimated of the issue and
may be provided alternatives or countermeasures. The issue of false positives
exist since it is not feasible to always incorporate the application specific
notions of anomaly into anomaly detection models. This effectively aids the
decision support system user as well the agents whose data is being assessed.
In practical scenarios, systems that utilize algorithmic recourse do require
transparency and human surveillance to ensure they are used as intended. There
is a minimal risk, as in any system, that unscrupulous personnel who are
insiders might utilize algorithmic recourse to provide alternatives to
nefarious agents enabling them avoid detection. This however can be eliminated
through correct operational and access protocols where such a system is
deployed.


# Synchronization Conditions for Nonlinear Oscillator Networks

Sanjeev Kumar Pandey, Shaunak Sen, Indra Narayan Kar
This work was not supported by any organization. The authors are with the
Department of Electrical Engineering, IIT Delhi, New Delhi, 110016,
India(e-mail:{sanjeev; shaunak.sen; ink} @ee.iitd.ac.in).

###### Abstract

Understanding conditions for the synchronization of a network of
interconnected oscillators is a challenging problem. Typically, only
sufficient conditions are reported for the synchronization problem. Here, we
adopted the Lyapunov-Floquet theory and the Master Stability Function approach
in order to derive the synchronization conditions for a set of coupled
nonlinear oscillators. We found that the positivity of the coupling constant
is a necessary and sufficient condition for synchronizing linearly full-state
coupled identical oscillators. Moreover, in the case of partial state
coupling, the asymptotic convergence of volume in state space is ensured by a
positive coupling constant. The numerical calculation of the Master Stability
Function for a benchmark two-dimensional oscillator validates the
synchronization corresponding to the positive coupling. The results are
illustrated using numerical simulations and experimentation on benchmark
oscillators.

##  I Introduction

Synchronization of oscillators is a phenomenon that occurs in a wide variety
of natural and engineered contexts. Examples include synchronized oscillations
in biology [1, 2, 3, 4], neuroscience [5, 6, 7, 8], and electronics [9, 10,
11, 12, 13, 14]. Understanding conditions for synchronization to occur would
benefit their analysis and design.

A significant advance in obtaining conditions for synchronization was the
Master Stability Function approach [15]. Through a computation of the Floquet
multipliers of nodes of a network of oscillators, conditions for
synchronization could be verified numerically. There have been several
attempts to obtain theoretical conditions for synchronization [16, 17, 18, 19,
20, 21, 22, 23], but these are often conservative and provide only sufficient
conditions. For example, the synchronization condition in the case of coupled
Van der Pol oscillators, obtained using a Lyapunov function approach and a
contraction theory-based approach, is that the critical coupling gain depends
on the system parameter and the connectivity graph.

In addition to numerical studies, electronic implementations of oscillator
networks have been investigated as more realistic models that often show rich
dynamics [12, 13, 11, 23, 24, 25]. These studies provide important foundations
for the analysis and design of synchronization. However, the problem of
obtaining tighter bounds on the coupling strength for an oscillator network is
generally unresolved.

Here, we aim to obtain the synchronizing condition for a set of coupled
nonlinear oscillators. To derive these conditions, we used the Lyapunov-
Floquet Theory and the Master Stability Function approach. The key finding is
that the positivity of the coupling constant is a necessary and sufficient
condition for the synchronization of a network of identical oscillators
connected linearly in a full-state fashion. A positive coupling constant
assures that the volume in the state space converges to zero asymptotically
for a set of oscillators with partial state linear coupling. Moreover, a
numerical computation of the Master Stability Function for a benchmark two-
dimensional oscillator confirm the synchronization behavior for positive
coupling. These findings are demonstrated on benchmark oscillators by
numerical simulations, LT SPICE, and electronic implementation.

##  II Mathematical Background

This section provides a brief overview of the relevant mathematical theory
that is used in this study. Consider a nonlinear dynamical system

$\dot{x}_{i}=f\left(x_{i}\right),\text{ }x_{i}\in\mathbb{R}^{m},$ (1)

with a limit cycle $x_{s}(t)$ with time period $T$. Linearizing this system
around $x=x_{s}(t)$ gives the linear time-periodic system

$\dot{y}=A(t)y,\text{ }y\in\mathbb{R}^{m},$ (2)

$A(t)=Df(x_{s}(t))$, where $Df$ is the Jacobian of $f$ at $x(t)=x_{s}(t)$.

###  II-A Floquet Theory

Floquet theory was developed to describe the behavior of a system of linear
differential equations with time-periodic coefficients [26]. The theory
introduces the concept of the Floquet multipliers ($\mu$), which are
eigenvalue-like entities that describe the exponential growth or decay of
solutions over one time period of the system. The fundamental solutions denote
the state transition by $\phi(t,t_{0})$. Define $R$ by $e^{RT}=\phi(T,0)$ and
$P^{-1}(t)=\phi(t,0)e^{-RT}$. Then $P^{-1}(t)$ is periodic with period $T$ as
$P^{-1}(t+T)=\phi(t+T,0)e^{-R(t+T)}=\phi(t+T,T).\phi(T,0).e^{-RT}.e^{-Rt}=\phi(t,0).e^{-Rt}=P^{-1}(t)$.
In the co-ordinates $z=P(t)y$, the system is the linear time invariant system
$\dot{z}=Rz$ as $\dot{P}(t)P^{-1}(t)+P(t)A(t)P^{-1}(t)=R.$ The eigenvalues of
$R$, called the Floquet exponents, determine the stability. Floquet
multipliers are the eigenvalues of $e^{RT}$. For asymptotic stability, one of
the Floquet multiplier has absolute value as $1$, representing perturbations
along the limit cycle, and others have moduli strictly less than $1$.

###  II-B Master Stability Function

Consider $n$ identical oscillators

$\dot{x}_{i}=f\left(x_{i}\right),\text{ }x_{i}\in\mathbb{R}^{m},\text{
}i\in\\{1,2,...,n\\},$ (3)

diffusively and identically coupled in a fully connected network

$\dot{x}_{i}=f\left(x_{i}\right)+K\sum_{j\in{N_{i}}}G_{ij}H\left(x_{j}\right),$
(4)

where $K$ is the coupling constant, $N_{i}$ is the neighborhood of the
oscillator $i$, $G=[a_{ij}]$ is the graph Laplacian and $H$ is an matrix
containing information on the coupled variables. For linear coupling, the
coupling between an oscillator $i$ and an oscillator $j\in N_{i}$ is
$K(x_{j}-x_{i})$. A simplified representation is

$\dot{X}=I_{n}\otimes f(X)+KG\otimes H(X),$ (5)

where $X=\left[x_{1}^{T},x_{2}^{T},\ldots,x_{n}^{T}\right]^{T}$ is a
$nm$-dimensional state vector, and $\otimes$ is the Kronecker product.
Linearisation around the state
$X_{s}=\left[x_{s}^{T},x_{s}^{T},\ldots,x_{s}^{T}\right]^{T}$ yields

$\dot{Y}=\left[I_{n}\otimes Df\left(x_{s}\right)-KG\otimes
DH\left(x_{s}\right)\right]Y,$ (6)

where $DH(x_{s})$ is the Jacobian of $H$ at $x(t)=x_{s}(t)$. Equation (6) can
be written in block matrix form by introducing a transformation
$\zeta=(Q\otimes I_{n})Y$, where $Q$ is the transformation matrix. On
diagonalising $G$ through the transformation $Q$, where
$Q^{-1}GQ=\text{diag}\left[0,\lambda_{2},\lambda_{3},...,\lambda_{n}\right]$,
$0<\lambda_{2}\leq\lambda_{3}\leq,...,\leq\lambda_{n}$. In transformed co-
ordinates, the dynamics are in block form, with each block representing an
eigenmode

$\dot{\zeta_{i}}=[DF\left(x_{s}\right)-K\lambda_{i}DH\left(x_{s}\right)]\zeta_{i},\text{
}i=1,2,...,n.$ (7)

$\zeta_{i}$ is referred as the synchronization mode. The dynamics
corresponding to the first eigenmode is $\dot{\zeta}_{1}=DF(x_{S})\zeta_{1}$,
which is the same as the linearised dynamics around the limit cycle. Assuming
this mode to be stable, the stability of the synchronized state is determined
by Floquet multipliers of all other eigenmodes. The Master Stability Function
($\mu_{max}(\lambda)$) is the largest non-unity Floquet multiplier of the
system matrix $[DF\left(x_{s}\right)-K\lambda DH\left(x_{s}\right)]$. The
necessary and sufficient condition for the stability of the limit cycle is
$\mu_{max}(\lambda)<1$.

##  III Condition For Synchronization

This section investigates the condition for synchronization of identical
oscillators coupled identically and linearly. Synchronization of two or more
identical oscillators is defined as the situation when the difference between
their corresponding states is asymptotically zero.

###### Theorem 1

A network of identical oscillators (4) coupled identically and linearly in
full-state fashion synchronizes if and only if $K>0$.

###### Proof:

The proof relies on the Lyapunov-Floquet transformation
$P(t)\in\mathbb{R}^{m\times m}$ [27] and the fact that for identical linear
full-state coupling, $DH=I.$ Applying the transformation $P(t)$ at each
eigenmode $\lambda_{i}$ (refer (7)) results in

$\dot{Z}_{i}=(R-KP(t)DH(x_{s})P^{-1}(t))Z_{i},\text{ }i=1,2,...,n,$ (8)
$\dot{Z}_{i}=(R-K\lambda_{i}I)Z_{i},\text{ }i=1,2,...,n,$ (9)

where $R$ is defined in subsection II-A. For $K=0$, system dynamics (7) become
uncoupled, and the matrix $R\in\mathbb{R}^{m\times m}$ has one eigenvalue at
$0$ and other $(m-1)$ eigenvalues with strictly negative real parts. For
coupled dynamics in (9), the eigenvalues are the roots of the polynomial,

$\det(sI-(R-K\lambda_{i}I))=0,$ $\Rightarrow\det((s+K\lambda_{i})I-R)=0.$

The roots of the above characteristic equation depend on the coupling strength
$K$ and the eigenmode $\lambda_{i}$. To prove sufficiency, first assume that
$K>0$ and define $s+K\lambda_{i}=\alpha\text{ }\Rightarrow
s=\alpha-K\lambda_{i}$. For $K=0,\text{ }s=\alpha$, and the eigenvalues are
same as those of $R$. One can conclude that the set of eigenvalues of
$(R-K\lambda_{i}I)$ are the same as the set of the eigenvalues of $R$ except
offset by $-K\lambda_{i}$. All eigenvalues of $R$ corresponding to each
synchronization modes ($\lambda_{i}$) have a strictly negative real part
except for one eigenvalue of the first eigenmode ($\lambda_{1}=0$), which is
zero.

To prove the necessary condition, we suppose that $K<0$. Define
$s+K\lambda_{i}=\alpha\text{ }\Rightarrow s=\alpha-K\lambda_{i}$. We use the
method of contradiction and assume that the synchronization is achieved if
$K<0$. The matrix $R\in\mathbb{R}^{m\times m}$ has one eigenvalue $0$ and
other $(m-1)$ eigenvalues with strictly negative real parts. However, as
$K<0$, the first eigenvalue of matrix $R$, which is zero, moves to the right
half of the complex plane for all eigenmode $\lambda_{i}$, ($i=2,...,n$). It
means that at least one of the eigenvalues of ($R-K\lambda_{i}I$)
corresponding to $i^{th}$ eigenmode has a positive real part. This contradicts
the supposition. ∎

###### Remark 1

Using the Lyapunov-Floquet transformation, the necessary and sufficient
condition on the coupling gain $K$ for the synchronization of $n$ oscillators
has been obtained to be $K>0$. This condition is non-conservative as compare
to existing results. For example, using a contraction theory-based approach, a
sufficient condition for the synchronization of $n$ coupled oscillators was
found to be $K>\frac{\alpha}{n}$, where $\alpha$ is an upper bound of a matrix
measure of the associated Jacobian [17]. Moreover, using a Lyapunov function-
based approach, a similar sufficient condition for the synchronization of $n$
coupled Van der Pol oscillators was found to be $K>\mu/n$, where $\mu$ is the
system parameter[23].

Theorem 1 is derived in the case that all states are coupled. When only some
states are coupled (referred to as a partially coupled system). In general,
$DH=diag[d_{1},d_{2},...,d_{n}]$ where $d_{i}\in\\{0,1\\}$ $i=1,2,...,n$ and
not all $d_{i}=1$ simultaneously. The above proof doesn’t work because the
matrix multiplication of $P(t)$ and $DH$ is not necessarily commutative.
However, a sufficient condition on the time evolution of the determinant of
the state transition matrix is obtained using Abel-Jacobi-Liouville (AJL)
identity [28].

###### Theorem 2

For the network of identical oscillators (4) coupled linearly in a partial
state fashion, $K>0$ implies
$\lim\limits_{t\to\infty}\operatorname{det}\phi(t,t_{0})=0$, where
$\phi(t,t_{0})$ is the state transition matrix of (7).

###### Proof:

The AJL identity for an uncoupled system (2) is

$\operatorname{det}\phi\left(t,t_{0}\right)=\exp\left[\int_{t_{0}}^{t}{\operatorname{tr}(DF(x_{s}(\tau))d\tau)}\right],$
(10)

where $\operatorname{det}\phi\left(t,t_{0}\right)$ is the state transition
matrix. The determinant of the transition matrix can be interpreted as a
measure of the volume in the phase space. For the linearized partial state
coupled system (7), the state transition matrix corresponding to each
eigenmode follows

$\operatorname{det}\phi\left(t,t_{0}\right)=\exp\left[\int_{t_{0}}^{t}{\operatorname{tr}(DF(x_{s}(\tau))-K\lambda_{i}DH(x_{s}))d\tau}\right],$
(11)
$=\exp\left[\int_{t_{0}}^{t}{(\operatorname{tr}(DF(x_{s}(\tau)))d\tau)}\right]\exp(-K\lambda_{i}(t-t_{0})).$
(12)

For positive coupling gain $K>0$, (12) ensures the convergence of
$\operatorname{det}\phi\left(t,t_{0}\right)$. In other words, for positive
values of $K$, $\operatorname{det}\phi\left(t,t_{0}\right)$ $\rightarrow 0$ as
$t\rightarrow\infty$. ∎

We further used a numerical approach to study the synchronization behavior of
a second-order partial-state coupled system (7). For this purpose, the Master
Stability Function was calculated numerically [15]. We found the maximum
Floquet multiplier is a decreasing function of the coupling strength ($K$).
This has been shown in Fig. 1. As the maximum Floquet multiplier is less than
one for $K>0$, the computation suggests synchronization.

Figure 1: Numerical computation of the Master Stability Function for a coupled
Van der Pol oscillator shows that the maximum Floquet multiplier decreases as
the coupling strength ($K$) increases. Calculating the Floquet multiplier
numerically required calculating the period and limit cycle. Solve the matrix
differential equation over one period given the Identity matrix as the initial
condition

##  IV Numerical Simulations

The present section uses numerical simulations to report the synchronization
behavior for both full-state and partial-state coupling. Here, oscillators
were coupled linearly. As case studies, two benchmark oscillators are
considered— the Van der Pol oscillator [2] and the repressilator [29, 30, 31].

Example 1: Van der Pol oscillator. The Van der Pol oscillator is a nonlinear
oscillator that exhibits limit cycle behavior and is frequently employed to
simulate self-sustained oscillations. This oscillator describes the behavior
of a nonlinearly damped oscillator and is modeled as

$\begin{split}\dot{x}_{1}&=x_{2},\\\
\dot{x}_{2}&=\mu(1-x_{1}^{2})x_{2}-x_{1},\end{split}$

where $\mu$ is a parameter. We investigated synchronization for a coupling
gain $K>0$. The graph Laplacian and the linearized coupling matrix of three
coupled oscillators with full-state coupling are, respectively

$G=\begin{pmatrix}2&-1&-1\\\ -1&2&-1\\\ -1&-1&2\end{pmatrix}\text{ and
}DH=\textbf{I}_{2\times 2}.$ (13)

We denote $x_{ij}$ as the $j^{th}$ state of the $i^{th}$ oscillator, where
$j\in\\{1,2\\}$ and $i\in\\{1,2,3\\}$. In the absence of coupling, all three
oscillators oscillate independently between $0-20$ sec ( Fig. 2(a)). After
$20$ sec, the coupling gain ($K=1$) is turned on, and the oscillators get
synchronized (Fig. 2(a)).

In partial-state coupling, the graph Laplacian matrix of coupled three
oscillators was the same as (13), but the linearized coupling matrix $DH$ was
$\begin{pmatrix}0&0\\\ 0&1\end{pmatrix},$ This means that only the second
state from each oscillator is coupled. As the coupling gain ($K=1$) is turned
on at $t=20$sec, all three oscillators get synchronized (Fig. 2(b)).

Figure 2: Numerical simulation of three coupled Van der Pol oscillators for
$\mu=1$, $\text{initial condition:}(x_{10}=[0,1],x_{20}=[2,3],x_{30}=[4,5])$
and $K=1$. (a) Synchronization with full-state coupling. (b) Synchronization
with partial-state coupling

Example 2: Repressilator.

This oscillator models a biomolecular oscillator with three proteins
inhibiting each other’s expression in a cyclic fashion. The mathematical model
of a repressilator is

$\begin{split}\frac{dm_{j}}{dt}&=-m_{j}+\frac{\alpha}{1+p_{j-1}^{n}}+\alpha_{0},\\\
\frac{dp_{j}}{dt}&=-\beta(p_{j}-m_{j}),\hskip 7.11317ptj=1,2,3,\end{split}$
(14)

where $p_{i}$ and $m_{i}$ are the concentrations of proteins and the
concentration of the corresponding mRNA, respectively. The parameters are
$\alpha$, $\alpha_{0}$, $\beta,$ and $n$, with $p_{0}\equiv p_{3}$. In
particular, $\frac{\alpha}{1+p_{i}^{n}}$ represents the nonlinear sigmoidal
function.

Three identical repressilators were coupled with full-state and partial-state
coupling. For full-state coupling, the graph Laplacian matrix is similar to
(13), but the linearized coupling matrix is $I_{3\times 3}$. We denote
$p_{ij}$ as the $j^{th}$ state of the $i^{th}$ oscillator, where
$j\in\\{1,2,3\\}$ and $i\in\\{1,2,3\\}$. Initially, all the repressilators
oscillate independently (Fig. 3). As the coupling gain ($K=1$) is turned on at
$t=20$sec, the oscillators acheived synchronization(Fig. 3(a)).

Similar behavior has been observed for the partial-state coupling. The graph
Laplacian matrix is same as (13) and the linearized coupling matrix is
$DH=Diag([0\hskip 5.69046pt1\hskip 5.69046pt0\hskip 5.69046pt1\hskip
5.69046pt0\hskip 5.69046pt1])$. As the coupling gain ($K=1$) is turned on at
$t=20$sec, oscillators get synchronized 3(b).

Figure 3: Numerical simulation of three coupled repressilator for different
parameters value $\alpha=1000$, $\alpha_{0}=1$, $\beta=5$, $n=2$,
$\text{initial
condition:}(p_{10}=[0,1,0,3,0,5],p_{20}=[0,7,0,9,0,11],p_{30}=[0,13,15,17,4,6])$
and $K=1$. (a) Synchronization with full-state coupling (b) Synchronization
with partial-state coupling.

##  V Experimental Results

To complement the numerical simulations, an electronic testbed was designed to
investigate the synchronization behaviour of nonlinear oscillators.

Symbol | Parameter | Value |  Units
---|---|---|---
$R_{i}$ | resistor | $1.0\pm 5\%$ |  $K\Omega$
$R_{j}$ | resistor | $1.0\pm 5\%$ |  $M\Omega$
$C_{i}$ | capacitor | $1.0\pm 10\%$ |  $\mu F$
op-amp | UA741CN |  |
Potentiometer | variable resistor | $100.0\pm 5\%$ |  $K\Omega$
$V_{1}$ and $V_{2}$ | voltage source | $\pm 12$ |  V
Analog multiplier | AD633JNZ |  |
TABLE I: Electronic components in the experiment

Example 1: Van der Pol oscillator. The oscillator circuit was realized using
operational amplifiers, analog multipliers, resistors, and capacitors [32].
Op-amps introduce nonlinearities, while analog multipliers capture the
quadratic terms in the Van der Pol equation. Resistors and capacitors are used
to tune the time constants and obtain the oscillator behavior. The electronic
circuit of the Van der Pol oscillator is given in Fig. 4.

Figure 4: Circuit diagram of a Van der Pol oscillator. (a) LT Spice design,
(b) experimental implementation.

We investigated the synchronization behavior of three identical Van der Pol
oscillators connected in all-to-all topology. The LT Spice simulation, using
the components value listed in Table 1, is shown in Fig. 5(a). As depicted in
the figure, the Van der Pol oscillators exhibited independent oscillations
between $0-50$ sec. After $50$ sec, the coupling is turned on, and all three
oscillators get synchronized (Fig. 5(a)).

Figure 5: Synchronization of three Van der Pol oscillators coupled with
different physical parameters is given in Table 1. (a) LT Spice simulation and
(b) experiment.

Similar behavior was observed experimentally (Fig. 5(b)), where the coupling
gain was set by varying a potentiometer.

Example 2: Repressilator. Each oscillator was implemented using three op-amps,
resistors, and capacitors in a negative feedback ring topology (Fig. 6) [13,
14]. The nonlinear sigmoidal [33] characteristic of the op-amp was well suited
to the nonlinearity in the repressilator equations. The output voltage of each
of the op-amps was analogous to protein concentration. A resistor-capacitor
network was used to implement the integration. The component values and
biasing voltages are listed in Table 1.

Figure 6: Circuit diagram of repressilator. (a) LT Spice design, (b) coupling
circuit, (c) experimental implementation

We found that the circuit exhibited oscillations in the LT Spice simulations.
As the output voltage of node two ($V_{2}$) increases, it results in the
reduction of the output voltage $V_{3}$. This is because of the influence at
the inverting input of the corresponding op-amps and the crossing of the
threshold voltage, which is $0$ V. Similarly, as $V_{3}$ decreases, $V_{1}$
will increase, which in turn will decrease $V_{2}$. This process leads to
oscillatory behavior. The phase difference between successive nodes was 120 ∘,
as expected for a 3-node ring oscillator. All three repressilators were
connected in an all-to-all topology and were synchronized once the coupling
was activated at $t=0.034$ sec (Fig. 7(a)). Before the coupling was activated,
all the oscillators were oscillating independently (Fig. 7(a)). Similar
behaviour was observed experimentally (Fig. 7(b)), where the coupling gain was
set by varying a potentiometer.

Figure 7: Synchronization of three coupled repressilators with different
physical parameters is given in Table 1. (a) LT Spice simulation and (b)
experiment.

##  VI Conclusion

We proved a necessary and sufficient condition for synchronizing a network of
identical oscillators coupled linearly and in a full-state fashion. The
derived condition is that the coupling constant is positive. These conditions
are derived using the Lyapunov-Floquet theory and the Master Stability
Function approach. Moreover, in the case of partial state coupling, a positive
coupling constant ensures that the volume in the state space converges to zero
asymptotically. Numerical computation of the Master Stability Function for a
second-order oscillator showed that a positive coupling ensures
synchronization. These results are illustrated on benchmark oscillators using
numerical simulations, LT SPICE, and electronic implementation.

## Acknowledgments

The author would like to thank Prof. S. Janardhanan and Prof. Deepak Patil of
the control and automation group for their constant support and guidance. I
also thank Dr. Abhilash Patel and Shivanagouda Biradar for their valuable time
and suggestions.

## References

  * [1] S. Strogatz, “Sync: The emerging science of spontaneous order,” 2004.
  * [2] A. Pikovsky, J. Kurths, M. Rosenblum, and J. Kurths, Synchronization: a universal concept in nonlinear sciences.  No. 12, Cambridge university press, 2003.
  * [3] B. Ermentrout, “An adaptive model for synchrony in the firefly pteroptyx malaccae,” Journal of Mathematical Biology, vol. 29, no. 6, pp. 571–585, 1991.
  * [4] S. H. Strogatz, “Nonlinear dynamics and chaos,” 1996.
  * [5] E. M. Izhikevich, Dynamical systems in neuroscience.  MIT press, 2007.
  * [6] C. D. Brody and J. Hopfield, “Simple networks for spike-timing-based computation, with application to olfactory processing,” Neuron, vol. 37, no. 5, pp. 843–852, 2003.
  * [7] N. Kopell, “We got rhythm: Dynamical systems of the nervous system,” Notices of the AMS, vol. 47, no. 1, pp. 6–16, 2000.
  * [8] W. Singer, “Neuronal synchrony: a versatile code for the definition of relations?,” Neuron, vol. 24, no. 1, pp. 49–65, 1999.
  * [9] T. L. Carroll and L. M. Pecora, “Synchronizing nonautonomous chaotic circuits,” IEEE Transactions on Circuits and Systems II: Analog and Digital Signal Processing, vol. 40, no. 10, pp. 646–650, 1993.
  * [10] Y. Tang, A. Mees, and L. Chua, “Synchronization and chaos,” IEEE Transactions on Circuits and Systems, vol. 30, no. 9, pp. 620–626, 1983.
  * [11] C. W. Wu and L. O. Chua, “Synchronization in an array of linearly coupled dynamical systems,” IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, vol. 42, no. 8, pp. 430–447, 1995.
  * [12] Z. Liu, J. Ma, G. Zhang, and Y. Zhang, “Synchronization control between two Chua’s circuits via capacitive coupling,” Applied Mathematics and Computation, vol. 360, pp. 94–106, 2019.
  * [13] J. M. Buldú, J. García-Ojalvo, A. Wagemakers, and M. A. Sanjuán, “Electronic design of synthetic genetic networks,” International Journal of Bifurcation and Chaos, vol. 17, no. 10, pp. 3507–3511, 2007.
  * [14] S. Y. Shafi, M. Arcak, M. Jovanović, and A. K. Packard, “Synchronization of diffusively-coupled limit cycle oscillators,” Automatica, vol. 49, no. 12, pp. 3613–3622, 2013.
  * [15] L. M. Pecora and T. L. Carroll, “Master stability functions for synchronized coupled systems,” Physical Review Letters, vol. 80, no. 10, pp. 2109–2112, 1998.
  * [16] W. Wang and J.-J. E. Slotine, “On partial contraction analysis for coupled nonlinear oscillators,” Biological cybernetics, vol. 92, no. 1, pp. 38–53, 2005.
  * [17] G. Russo and J. J. E. Slotine, “Global convergence of quorum-sensing networks,” Physical Review E, vol. 82, no. 4, p. 041919, 2010.
  * [18] L. Scardovi and R. Sepulchre, “Synchronization in networks of identical linear systems,” in 47th Conference on Decision and Control (CDC), pp. 546–551, IEEE, 2008.
  * [19] M. Barahona and L. M. Pecora, “Synchronization in small-world systems,” Physical review letters, vol. 89, no. 5, p. 054101, 2002.
  * [20] H. Gu, K. Liu, and J. Lü, “Adaptive pi control for synchronization of complex networks with stochastic coupling and nonlinear dynamics,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 67, no. 12, pp. 5268–5280, 2020.
  * [21] A. Jadbabaie, N. Motee, and M. Barahona, “On the stability of the kuramoto model of coupled nonlinear oscillators,” in Proceedings of the 2004 American Control Conference, vol. 5, pp. 4296–4301, IEEE, 2004.
  * [22] E. N. Davison, B. Dey, and N. E. Leonard, “Synchronization bound for networks of nonlinear oscillators,” in 54th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pp. 1110–1115, IEEE, 2016\.
  * [23] S. K. Joshi, S. Sen, and I. N. Kar, “Synchronization of coupled benchmark oscillators: analysis and experiments,” International Journal of Dynamics and Control, vol. 10, no. 2, pp. 577–597, 2022.
  * [24] A. Wagemakers, J. M. Buldú, J. García-Ojalvo, and M. A. Sanjuán, “Synchronization of electronic genetic networks,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 16, no. 1, p. 013127, 2006\.
  * [25] M. Lodi, M. Forti, and M. Storace, “Stability analysis of the synchronous solution in arrays of memristive Chua’s circuits,” IEEE Transactions on Circuits and Systems II: Express Briefs, no. 5, pp. 1694–1698, 2023.
  * [26] V. A. Yakubovich and V. M. Starzhinskii, Linear differential equations with periodic coefficients.  New York: Wiley, 1975.
  * [27] S. Sinha, R. Pandiyan, and J. Bibb, “Liapunov-floquet transformation: Computation and applications to periodic systems,” Journal of Vibration and Acoustics, vol. 118, pp. 209–219, 1996.
  * [28] R. W. Brockett, Finite dimensional linear systems.  SIAM, 2015.
  * [29] M. B. Elowitz and S. Leibler, “A synthetic oscillatory network of transcriptional regulators,” Nature, vol. 403, no. 6767, pp. 335–338, 2000\.
  * [30] E. H. Hellen, E. Volkov, J. Kurths, and S. K. Dana, “An electronic analog of synthetic genetic networks,” PLoS One, vol. 6, no. 8, p. e23286, 2011.
  * [31] E. H. Hellen, S. K. Dana, B. Zhurov, and E. Volkov, “Electronic implementation of a repressilator with quorum sensing feedback,” PLoS One, vol. 8, no. 5, p. e62997, 2013.
  * [32] J. K. Roberge, Operational amplifiers: theory and practice, vol. 197.  Wiley New York, 1975.
  * [33] R. Gesztelyi, J. Zsuga, A. Kemeny-Beke, B. Varga, B. Juhasz, and A. Tosaki, “The hill equation and the origin of quantitative pharmacology,” Archive for history of exact sciences, vol. 66, no. 4, pp. 427–438, 2012.