# Guarantees for Spontaneous Synchronization on Random Geometric Graphs.

Pedro Abdalla, Afonso S. Bandeira and Clara Invernizzi
Department of Mathematics, ETH Zürich

###### Abstract

The Kuramoto model is an important mathematical model to study the
synchronization phenomena in complex networks. The rich theory behind the
rigorous treatment of the synchronization phenomenon has been connected to
many branches of science. In particular, the study of the energy function
associated with the Kuramoto model is relevant in the context of non-convex
optimization for machine learning and data science.

In this work, we use tools from random matrix theory, random graphs, and
statistics to derive rigorous guarantees for the random geometric graph on the
sphere to synchronize with probability of success tending to one as the number
of nodes tends to infinity. To the best of our knowledge, this is the first
rigorous result for the Kuramoto model on random geometric graphs.

##  1 Introduction

Spontaneous synchronization is an ubiquitous phenomenon in nature. The
experimental investigation of this phenomenon dates back to the work of
Huygens in 1665, who observed that two pendulums held by the same wooden beam
tend to swing in anti-phase after a certain amount of time, regardless their
initial position. Similar phenomena have been observed in many branches of
science [26]. In a seminal work [18], Kuramoto introduced the first
mathematical tractable model to study synchronization in a system of all-to-
all phase oscillators. Natural generalizations for complex networks were
considered in the literature, but their analysis remains not fully understood.
We point out that different versions of the Kuramoto model have been studied
with applications in electric power networks, wireless networks, neuroscience
[1, 10, 24, 9], to name a few.

In this work, we study the following model: Consider a graph $G$ with
adjacency matrix $A$ with entries $a_{ij}\in\\{0,1\\}$. Each node $i$ is
associated with a phase $\theta_{i}\in(0,2\pi]$ that varies with the time $t$
and a natural frequency $\omega_{i}\in\mathbb{R}$. The system of differential
equations below represents the Kuramoto model associated with the graph $G$

$\forall
i\in[n]:\frac{d\theta_{i}(t)}{dt}=w_{i}-\frac{K}{n}\sum_{i,j=1}^{n}a_{ij}\sin(\theta_{i}-\theta_{j}),$

where $K$ is the coupling strength.

In order to give a more intuitive description of the model, we describe an
alternative viewpoint due to Dörfle, Chertkov, and Bullo [10]. For simplicity,
in this manuscript, we only discuss the homogeneous Kuramoto model, where
$w_{i}$ is zero for all $1\leq i\leq n$. The coupled oscillators system can be
interpreted as a system of particles that tends to minimize the following
energy function

$E_{\omega}(\theta):=\sum_{i,j}a_{ij}(1-\cos(\theta_{i}-\theta_{j})),$

and the standard Cauchy problem associated with this energy function can be
written as

$\forall i\in[n]:\frac{d\theta_{i}(t)}{dt}=-\frac{\partial
E_{w}(\theta_{1},\ldots,\theta_{n})}{\partial\theta_{i}},\quad\theta_{i}=\theta_{i}(0).$

From the equation above, the intuitive meaning of Kuramoto model is that the
evolution moves according to the gradient flow of the energy function, i.e,
the direction of maximum decrease of $E_{\omega}$. With this gradient flow
formulation, we study the relationship between the graph $G$ and the
synchronization phenomenon associated to the topology of $G$. Mathematically,
this means that we are interested in determining whether the energy function
$E_{\omega}(\theta)$ admits a local minimum that is not the global minimum,
also called a spurious minimum. It is easy to see that the global minimum
corresponds to $\theta_{1}=\ldots=\theta_{n}=0$ (up to a global shift). The
challenge here is that the function is not convex so it is necessary to do a
deeper analysis of the optimization landscape associated with the energy
function. In fact, there is a significant body of work dedicated to the study
of the optimization landscape of non-convex functions because such problems
often appear in machine learning and data science (e.g. [13, 14]), some works
have been directly connected to the study of energy functions such as the one
treated in this work [2, 20].

###  1.1 Kuramoto Model on random graphs

A trivial necessary condition for the oscillators to synchronize is that the
underlying graph is connected. However, this is far from being sufficient.
There are many examples of networks that do not synchronize even if all nodes
are connected over half of the other nodes. There are also known sparse
networks that synchronize [27, 20]. Moreover, it is known that the addition of
a single edge can create spurious local minima [20].

The motivation to study the synchronization phenomenon in random graphs is to
study the behavior of oscillators on a typical graph. In this spirit, the
standard Erdős-Rényi random graphs, in which an edge is placed between every
pair of vertices independently with probability $p$, were considered to
quantify the typical behaviour of networks in the homogeneous Kuramoto model,
i.e, when all $\omega_{i}$ are equal to zero [20, 17].

In this work, we also focus on the homogeneous case, and analyse the existence
of spurious local minima when the graph $G$ is drawn at random with some
geometric structure behind, the so called random geometric graphs on the
sphere. These are graphs with $n$ vertices sampled by generating $n$ points in
the unit sphere, and placing an edge between every pair of vertices with
Euclidean distance smaller than a certain threshold. The random geometric
graph on the sphere is an interesting model because it captures the property
that two nodes are more likely to connect if both are connected with a third
node. Similar models have been studied in the literature [12, 19, 8], but we
are not aware of any rigorous mathematical guarantees.

In contrast to the Erdős-Rényi random graph model studied in [20, 17], the
adjacency matrix of the random geometric graph does not enjoy independence
among entries or columns, which requires different tools in their analysis.
Our main contribution is to give a rigorous treatment for the random geometric
graphs by exploring its connections with random matrix theory and recent work
on hypothesis testing among random graph models [21]. Our main result
establishes that a random geometric graph with $n$ nodes drawn uniformly at
random from the unit sphere $\mathbb{S}^{d-1}$ synchronizes, with high
probability, as long as $d$ is sufficiently large compared to $n$. We state
here a simplified statement

######  Theorem 1.

(Simplified Statement of the Main Result Theorem 4) Let $G$ be a random
geometric graph with $n$ nodes on the unit sphere $\mathbb{S}^{d-1}$ such that
each edge is connected with probability $p$. Then if $np=\Omega(\log^{2}(n))$
and $d=\tilde{\Omega}(n^{2}p^{2})$, the graph $G$ synchronizes with high
probability. Here $\tilde{\Omega}(.)$ hides polylogarithmic factors.

A more detailed version and its proof is presented in Section 3. We end up
this section with a discussion about our main result. It shows that if
$p=\tilde{\Theta}\left(\frac{1}{n}\right)$ and
$d=\tilde{\Omega}\left(n^{2}p^{2}\right)=\tilde{\Omega}\left(1\right)$, then
the random geometric graph on the sphere $\mathbb{S}^{d-1}$ with connectivity
probability $p$ synchronizes with high probability. On the other hand, it is
natural to expect that if $p$ grows and the other parameters are fixed, then
$G$ also synchronizes with high probability. If we could formally prove that
this monotonicity holds, then our result would imply that if we fix any
$p=\tilde{\Omega}\left(\frac{1}{n}\right)$ and
$d=\tilde{\Omega}\left(1\right)$, then the graph synchronizes. This is
essentially the best we can hope for (up to logarithmic factors) because below
this threshold the graph will not be connected with non-vanishing probability.

The challenge is that the proof technique relies in comparisons between the
geometric and Erdős-Rényi random graph models which impose upper bounds on
$p$. While one expects that, on average, larger $p$ helps synchronization
(even if it might make a geometric graph more easily distinguishable from an
Erdős-Rényi one), we were not able to establish such a monotonicity result,
and we leave it as an open problem. An important difficulty is that adding
edges, in a worst case sense, may destroy the synchronization property.

### Notation

The standard unit sphere in $\mathbb{R}^{d}$ is denoted by $\mathbb{S}^{d-1}$
as usual. The matrix $J_{n}$ represents the $n\times n$ all ones matrix and
$I_{n}$ denotes the $n\times n$ identity matrix. For vectors $u,v$ of same
size, the inner product $\langle u,v\rangle$ corresponds to the standard inner
product, $\|v\|_{2}$ denotes the standard Euclidean norm. For a matrix A,
$\|A\|$ denotes the standard spectral norm. We consider the standard "Big-O"
notation (Landau notation) and we write $\tilde{\Omega}(.)$ to hide log
factors.

##  2 Preliminaries on the Kuramoto Model

We will establish spontaneous synchronization of the graphs drawn from our
random graph model by showing that they satisfy, with high probability,
explicit properties that imply synchronization. In what follows we describe
two such sets of properties. One approach, developed to analyse the landscape
of certain optimization problems arising from the Burer-Monteiro relaxation to
semidefinite programming (see [5, 6]), yields a weaker estimate than the one
on this paper, and is available in an earlier version of our work [15].

A second approach, based on very recent remarkable progress by Kassabov,
Strogatz and Townsend [17] on the problem of synchronization of Erdős-Rényi
graphs, provides better sharper guarantees. They established a sufficient
condition for the graph to synchronize based on the so called complex order
parameters $\rho_{k}:=\frac{1}{n}\sum_{j=1}^{n}e^{ik\theta_{j}}$, for
$k\in\mathbb{N}$. These parameters encode information about the landscape of
$E(\theta)$, for example, $\rho_{1}=1$ corresponds to the state of pure
synchrony. Their main result in this direction is the following theorem

######  Theorem 2.

[17, Theorem 6] Let $G$ be a graph with $n$ vertices and fix $p_{o}$ such that
$0<p_{o}<1$. Let $A$ and $L$ be the associated adjacency, Laplacian matrix
respectively, and define $\bm{\Delta}_{A}=A-p_{o}J_{n}$,
$\bm{\Delta}_{L}=L+p_{o}J_{n}-np_{o}I_{n}$. If the three conditions below are
satisfied

  1. 1.

$\frac{\|\bm{\Delta}_{A}\|}{np_{o}}<\frac{1}{12}$,

  2. 2.

$\frac{\|\bm{\Delta}_{L}\|}{np_{o}}<\frac{1}{4}$,

  3. 3.

$\frac{\pi/4}{\sin^{-1}\left(\frac{12\|\bm{\Delta}_{A}\|}{np_{o}}\right)}>\frac{\log(n/6)}{\log(\frac{np_{o}}{2\|\bm{\Delta}_{L}\|}-1)}+1$,

then the homogeneous Kuramoto model on the graph $G$ globally synchronizes.

In the next section, we prove that these conditions are satisfied for random
geometric graphs on the sphere in the regime stated in Theorem 4.

##  3 The Random Geometric Graph on the Sphere

We start by formally defining the random geometric graph on the sphere. We
also include the definition of the Erdős-Rényi random graph for the sake of
completeness.

######  Definition 1 (Erdős-Rényi Random Graph [11]).

Let $n\in\mathbb{N}$ and consider $p$ a fixed probability. The random graph
$G(n,p)$ is drawn as follows: For each pair $\\{i,j\\}\in[n]\times[n]$ with
$i\neq j$, an edge is placed with probability $p$ independently.

######  Definition 2 (Random Geometric Graph on the Sphere [22, 23]).

Let $n\in\mathbb{N}$ and $X_{1},...,X_{n}\in\mathbb{S}^{d-1}$ be points drawn
independently and uniformly at random on the unit sphere $\mathbb{S}^{d-1}$.
We define the random geometric graph $G(n,p,d)$ as the graph on the vertex set
$[n]$ such that for $i\neq j$, $\\{i,j\\}\in[n]\times[n]$ is an edge of the
graph if and only if $\langle X_{i},X_{j}\rangle\geq t_{p,d}$, where
$t_{p,d}\in[-1,1]$ is related to $p$ by

$p=\mathbb{P}\left(\langle X_{i},X_{j}\rangle\geq t_{p,d}\right).$

Equivalently, $\\{i,j\\}$ is an edge of the graph if and only if the Euclidean
distance between $X_{i}$ and $X_{j}$ satisfies
$\|X_{i}-X_{j}\|_{2}\leq\sqrt{2-2t_{p,d}}$.

A significant body of work in statistics [7, 4, 3, 21] is dedicated to the
following question: If we receive one sample from a random graph that is
either sampled from $G(n,p,d)$ or $G(n,p)$, is it possible to determine, with
high probability, from which graph it was sampled? In other words, can we do
hypothesis testing with the probability of error going to zero as $n$ goes to
infinity?

Clearly, if we cannot distinguish between the $G(n,p)$ and $G(n,p,d)$, then
the random graph $G(n,p,d)$ must synchronize with high probability since we
know that the $G(n,p)$ does. Recall that the definition of total variation
distance between two random graphs $G$ and $G^{\prime}$ with set of vertices
$V$ is defined by

$\operatorname{TV}(G,G^{\prime})=\max_{\mathcal{G}}|\mathbb{P}(G\in\mathcal{G})-\mathbb{P}(G^{\prime}\in\mathcal{G})|,$

where $\mathcal{G}$ is any set of graphs over $V$. In our context, the
relevant result is the following

######  Theorem 3.

[21, Theorem 1.2] For any fixed constant $\alpha>0$, if
$\frac{\alpha}{n}<p<\frac{1}{2}$ and $d=\tilde{\Omega}(n^{3}p^{2})$, then

$\lim_{n\rightarrow\infty}\operatorname{TV}\left(G(n,p),G(n,p,d)\right)=0.$

The latter result directly implies that, in the regime mentioned above for
$d$, the random graph $G(n,p,d)$ synchronizes as long as the $G(n,p)$ does.
However, it is not clear if the random graph $G(n,p,d)$ synchronizes, with
high probability, below this regime. Our result shows that it is indeed
possible.

######  Theorem 4.

(Main Result) Let $G$ be a random graph sampled from $G(n,p,d)$. For any
constant $\alpha>0$, there are absolute constants $C_{p},C_{d}>0$ such that if
$d\geq C_{d}(n^{2}p^{2}+\log^{4}(n))\log^{4}(n)$,
$\frac{\alpha}{n}<p<\frac{1}{2}$ and $np\geq C_{p}\log(n)^{2}$, then, with
probability at least $1-n^{-\Omega(1)}$, the graph $G$ synchronizes.

Our proof strategy is to verify that the three conditions in Theorem 2 are
satisfied for $G(n,p,d)$ in the regime above. We start by collecting a few
preliminary results for random graphs. The first one is a standard fact, but
we could not locate the exact statement in the literature.

######  Lemma 1.

Let $A$ be the adjacency matrix of a graph sampled from the $G(n,p)$ with
i.i.d Bernoulli $p$ random variables above the diagonal. Suppose that
$np>\log(n)$. Then there exist absolute constants $K,c>0$ such that, with
probability at least $1-e^{-t^{2}/c}$,

$\|A-\mathbb{E}A\|\leq K\sqrt{np}+t.$

We provide a short sketch of the proof.

###### Proof.

By Theorem 2.1 in [25], we know that $\mathbb{E}\|A-\mathbb{E}A\|$ is at most
a constant times the expectation of the maximum between the Euclidean norms of
the columns. Based on this fact, an easy computation gives that
$\mathbb{E}\|A-\mathbb{E}A\|\leq C\sqrt{np}$. To get the tail bound, we apply
Talagrand convex Lipschitz concentration inequality [28, Theorem 5.2.16]. ∎

The second preliminary proposition that we need is the following.

######  Proposition 1.

[21, Proposition 1.3] For any constant $\alpha>0$, there exist constants
$C_{d},C_{\varepsilon}>0$ such that if $\frac{\alpha}{n}\leq p\leq\frac{1}{2}$
and $d\geq C_{d}(n^{2}p^{2}+\log^{4}n)\log^{4}n$, then, for any
$\varepsilon\geq C_{\varepsilon}\sqrt{\frac{1}{d}(np+\log n)\log^{4}n}$, there
exists a probability space such that $G_{-}\sim G(n,p(1-\varepsilon))$, $G\sim
G(n,p,d)$ and $G_{+}\sim G(n,p(1+\varepsilon))$ satisfy, with probability at
least $1-n^{-\Omega(\log n)}$, $G_{-}\subseteq G\subseteq G_{+}$.

For simplicity, we assume that both random graphs have self-loops, i.e, every
node may be connected with itself. We remark that this assumption does not
change the energy function. We isolate the main technical proof in the
following lemma, it verifies the first condition of Theorem 2.

######  Lemma 2.

Let $A$ be the adjacency matrix of a graph sampled from $G(n,p,d)$. There
exists absolute constants $K,C_{\varepsilon},C_{d}>0$ such that if
$np>\log(n)^{2}$ and $d\geq C_{d}(n^{2}p^{2}+\log^{4}(n))\log^{4}(n)$, then,
with probability at least $1-n^{-\Omega(1)}$,

$\frac{\|\bm{\Delta}_{A}\|}{np}<\sqrt{\frac{2}{np}}\left(K+\max\left(4\frac{C_{\varepsilon}}{\sqrt{C_{d}}},24\right)\right).$

###### Proof.

We choose $C_{\varepsilon},C_{d}>0$ to be the constants in Proposition 1 and
$K>0$ as in Lemma 1. We fix

$\varepsilon:=\max\left(\frac{C_{\varepsilon}}{\sqrt{C_{d}}},6\right)\sqrt{\frac{2}{np}}\geq
C_{\varepsilon}\sqrt{\frac{1}{d}(np+\log n)\log^{4}n}.$ (1)

The inequality above follows easily from the choice of $d$ and $np$. The need
for the constant $6$ will become clear later. Notice that the hypothesis of
Proposition 1 is satisfied, so if $G^{-}\coloneqq G(n,p(1-\varepsilon))$ and
$G^{+}\coloneqq G(n,p(1+\varepsilon))$ are Erdős-Rényi random graphs, then
there exists a probability space such that $G^{-}\subseteq G\subseteq G^{+}$
with the required probability. Now we define $L\coloneqq D-A$, the Laplacian
matrix associated to $G$, and $L^{-}\coloneqq D^{-}-A^{-}$, $L^{+}\coloneqq
D^{+}-A^{+}$ the Laplacian matrices of $G^{-}$ and $G^{+}$ respectively.
Recall that the quadratic form of a Laplacian matrix can always be written as
$\frac{1}{2}\sum_{\\{i,j\\}\in E(G)}(v_{i}-v_{j})^{2}$. Hence the property
$G^{-}\subseteq G\subseteq G^{+}$ implies in the following monotone relation:
For all vectors $v$,

$v^{T}L^{-}v\leq v^{T}Lv\leq v^{T}L^{+}v.$ (2)

Therefore, for an arbitrary unit vector $v\in\mathbb{R}^{n}$, we have

$\displaystyle v^{T}(A-\mathbb{E}A)v$ $\displaystyle=v^{T}Av-
pv^{T}J_{n}v=v^{T}Dv-v^{T}Lv-pv^{T}J_{n}v$ $\displaystyle\leq
v^{T}D^{+}v-v^{T}L^{-}v-pv^{T}J_{n}v+v^{T}\mathbb{E}[A^{-}]v-v^{T}\mathbb{E}[A^{-}]v$
$\displaystyle=v^{T}(D^{+}-D^{-})v+v^{T}(A^{-}-\mathbb{E}[A^{-}])v-p\varepsilon
v^{T}J_{n}v$ $\displaystyle\leq
v^{T}(D^{+}-D^{-})v+v^{T}(A^{-}-\mathbb{E}[A^{-}])v,$

we used here the self-loop simplification, $\mathbb{E}A=pJ_{n}$ and
$\mathbb{E}[A^{-}]=p(1-\varepsilon)J_{n}$. In the second line we use property
(2) and the fact that $G\subseteq G^{+}$ implies $D_{ii}\leq D_{ii}^{+}$ for
every $i\in[n]$. The last line then follows because $v^{T}J_{n}v\geq 0$ for
every $v\in\mathbb{R}^{n}$. In an analogous way, we obtain

$v^{T}(A-\mathbb{E}A)v\geq
v^{T}(D^{-}-D^{+})v+v^{T}(A^{+}-\mathbb{E}[A^{+}])v.$

Putting everything together, we obtain

$\displaystyle\|A-\mathbb{E}A\|$
$\displaystyle=\sup_{\|v\|_{2}=1}|v^{T}(A-\mathbb{E}A)v|$ (3)
$\displaystyle\leq\|D^{+}-D^{-}\|+\max\left\\{\|A^{-}-\mathbb{E}[A^{-}]\|,\|A^{+}-\mathbb{E}[A^{+}]\|\right\\}.$

By Lemma 1, we obtain that

$\displaystyle\|A^{-}-\mathbb{E}[A^{-}]\|$ $\displaystyle\leq
K\sqrt{np(1-\varepsilon)},$ $\displaystyle\|A^{+}-\mathbb{E}[A^{+}]\|$
$\displaystyle\leq K\sqrt{np(1+\varepsilon)},$

with the desired probability and $K>0$ is the constant in Lemma 1. For the
first term in the right hand side of (3), we claim that $\|D^{+}-D^{-}\|\leq
4p\varepsilon n$ with the desired probability. Indeed, each entry in the
diagonal matrix $D^{+}-D^{-}$ is the sum of $n$ independent Bernoulli random
variables:

$(D^{+}-D^{-})_{ii}=\sum_{j=1}^{n}X_{ij}\quad\text{for}\quad
X_{ij}=\begin{cases}1\ \text{if}\ A^{+}_{ij}=1\ \text{and}\ A^{-}_{ij}=0,\\\
0\ \text{else},\end{cases}$

and

$\mathbb{P}(X_{ij}=1)=1-(\mathbb{P}(A^{+}_{ij}=0)+\mathbb{P}(A^{-}_{ij}=1))=1-(1-p(1+\varepsilon)+p(1-\varepsilon))=2p\varepsilon.$

The expectation of each entry of $D^{+}-D^{-}$ is
$\mathbb{E}[(D^{+}-D^{-})_{ii}]=2p\varepsilon n$. We now apply Chernoff’s
deviation inequality [28, 16] to get

$\mathbb{P}((D^{+}-D^{-})_{ii}\geq 4p\varepsilon
n)\leq\exp\left(-\frac{p}{3}\varepsilon n\right).$

By union bound, definition of $\varepsilon$ in (1) and the fact that
$np>\log^{2}(n)$, we obtain that the claim holds with probability at least

$1-\exp\left(-\frac{p\varepsilon n}{3}+\log n\right)\geq
1-\exp\left(-2\sqrt{2}\log n+\log n\right),$

the probability converges to one even if $n=2$. Now, for the range of $d$ in
the hypothesis of this lemma, the following holds with the desired
probability,

$\begin{split}\frac{\|\bm{\Delta}_{A}\|}{np}&\leq
K\frac{\max(\sqrt{1-\varepsilon},\sqrt{1+\varepsilon})}{\sqrt{np}}+\frac{4p\varepsilon
n}{np}\\\ &\leq K\sqrt{\frac{2}{np}}+4\varepsilon.\end{split}$

Recall the choice of $\varepsilon$ in the beginning to conclude the proof. ∎

We now present the proof of the main theorem.

###### Proof of Theorem 4.

As explained above, we verify that the three conditions listed in Theorem 2
are satisfied. Thanks to Lemma 2, we know that
$\frac{\|\bm{\Delta}_{A}\|}{np}<\frac{1}{12}$ with the required probability.
To verify the second condition, we apply triangular inequality to obtain

$\|\bm{\Delta}_{L}\|=\|D-A-\mathbb{E}D+\mathbb{E}A\|\leq\|D-\mathbb{E}D\|+\|A-\mathbb{E}A\|<\|D-\mathbb{E}D\|+\frac{1}{12}.$

To deal with the first term in the right hand side, observe that the entries
of $D-\mathbb{E}D$ are only nonzero on the diagonal and
$(D-\mathbb{E}D)_{ii}=\sum_{j=1}^{n}\mathds{1}_{\langle X_{i},X_{j}\rangle\geq
t_{p,d}}-p$. Conditioned on $X_{i}$, the distribution of
$(D-\mathbb{E}D)_{ii}$ is a centered binomial random variable and we apply
Chernoff’s deviation again to obtain, for every $\delta>0$,

$\mathbb{P}\left(\left|(D-\mathbb{E}D)_{ii}\right|\geq\delta np\right)\leq
2\exp\left(-\frac{\delta^{2}np}{3}\right).$

We choose $\delta=\frac{\log n}{\sqrt{np}}$ and union bound over all $i\in[n]$
to obtain that, with the required probability,

$\frac{\|\bm{\Delta}_{L}\|}{np}<\frac{\log
n\sqrt{np}}{np}+\frac{1}{12}\leq\frac{1}{\sqrt{C_{p}}}+\frac{1}{12}<\frac{1}{8},$

for a sufficiently large $C_{p}>0$. The reason why $\frac{1}{8}$ instead of
$\frac{1}{4}$ will become clear later. We proceed to verify the last
condition. It is easy to see that for a sufficiently large $C_{p}>0$, we have

$\frac{12\|\bm{\Delta}_{A}\|}{np}\leq\frac{1}{4\log n}.$

Since the function $\sin^{-1}(x)$ is non-increasing and satisfies
$\sin^{-1}(x)<2x$ for $x\geq 0$, we obtain that

$\sin^{-1}\left(\frac{12\|\bm{\Delta}_{A}\|}{np}\right)\leq\frac{1}{2\log n}.$
(4)

On the other hand, we know that $\frac{np}{2\|\bm{\Delta}_{L}\|}>4$ and then

$\log\left(\frac{np}{2\|\bm{\Delta}_{L}\|}-1\right)>\log(3).$ (5)

Combining (4) with (5), we get

$\frac{\pi/4}{\sin^{-1}\left(\frac{12\|\bm{\Delta}_{A}\|}{np}\right)}\geq\frac{\pi}{2}\log
n>\log\left(\frac{n}{6}-3\right)+1>\frac{\log(n/6)}{\log(\frac{np}{2\|\bm{\Delta}_{L}\|}-1)}+1.$

We conclude the proof by applying Theorem 2. ∎

## Acknowledgment

We would like to thank Alex Townsend and Tselil Schramm for insightful
discussions. In particular, we are indebted to Tselil Schramm for various
insights and clarifications with regards to properties of random geometric
graphs.

## References

  * [1] J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler.  The kuramoto model: A simple paradigm for synchronization phenomena.  Reviews of modern physics, 77(1):137, 2005.
  * [2] N. Boumal.  Nonconvex phase synchronization.  SIAM Journal on Optimization, 26(4):2355–2377, 2016.
  * [3] M. Brennan, G. Bresler, and D. Nagaraj.  Phase transitions for detecting latent geometry in random graphs.  Probability Theory and Related Fields, 178(3):1215–1289, 2020.
  * [4] S. Bubeck, J. Ding, R. Eldan, and M. Z. Rácz.  Testing for high-dimensional geometry in random graphs.  Random Structures & Algorithms, 49(3):503–532, 2016.
  * [5] S. Burer and R. D. Monteiro.  A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization.  Mathematical Programming, 95(2):329–357, 2003.
  * [6] S. Burer and R. D. Monteiro.  Local minima and convergence in low-rank semidefinite programming.  Mathematical programming, 103(3):427–444, 2005.
  * [7] L. Devroye, A. György, G. Lugosi, and F. Udina.  High-dimensional random geometric graphs and their clique number.  Electronic Journal of Probability, 16:2481–2508, 2011.
  * [8] A. Diaz-Guilera, J. Gómez-Gardenes, Y. Moreno, and M. Nekovee.  Synchronization in random geometric graphs.  International Journal of Bifurcation and Chaos, 19(02):687–693, 2009.
  * [9] F. Dörfler and F. Bullo.  Synchronization in complex networks of phase oscillators: A survey.  Automatica, 50(6):1539–1564, 2014.
  * [10] F. Dörfler, M. Chertkov, and F. Bullo.  Synchronization in complex oscillator networks and smart grids.  Proceedings of the National Academy of Sciences, 110(6):2005–2010, 2013.
  * [11] P. Erdős and A. Rényi.  On the evolution of random graphs.  Magyar Tud. Akad. Mat. Kutató Int. Közl., 5:17–61, 1960\.
  * [12] E. Estrada and G. Chen.  Synchronizability of random rectangular graphs.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(8):083107, 2015.
  * [13] R. Ge, J. D. Lee, and T. Ma.  Matrix completion has no spurious local minimum.  Advances in neural information processing systems, 29, 2016.
  * [14] R. Ge and T. Ma.  On the optimization landscape of tensor decompositions.  Advances in Neural Information Processing Systems, 30, 2017.
  * [15] C. Invernizzi.  Synchronization problems on random geometric graphs.  Master’s thesis, ETH Zürich (Available at request), 2022.
  * [16] S. Janson, A. Rucinski, and T. Luczak.  Random graphs.  John Wiley & Sons, 2011.
  * [17] M. Kassabov, S. H. Strogatz, and A. Townsend.  A global synchronization theorem for oscillators on a random graph.  arXiv preprint arXiv:2203.03152, 2022.
  * [18] Y. Kuramoto.  Self-entrainment of a population of coupled non-linear oscillators.  In International symposium on mathematical problems in theoretical physics, pages 420–422. Springer, 1975.
  * [19] I. Leyva, A. Navas, I. Sendiña-Nadal, J. M. Buldu, J. A. Almendral, and S. Boccaletti.  Synchronization waves in geometric networks.  Physical Review E, 84(6):065101, 2011.
  * [20] S. Ling, R. Xu, and A. S. Bandeira.  On the landscape of synchronization networks: A perspective from nonconvex optimization.  SIAM Journal on Optimization, 29(3):1879–1907, 2019.
  * [21] S. Liu, S. Mohanty, T. Schramm, and E. Yang.  Testing thresholds for high-dimensional sparse random geometric graphs.  arXiv preprint arXiv:2111.11316, 2021.
  * [22] G. Lugosi.  Lectures on combinatorial statistics.  47th Probability Summer School, Saint-Flour, pages 1–91, 2017.
  * [23] M. Penrose.  Random geometric graphs, volume 5.  OUP Oxford, 2003.
  * [24] F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths.  The kuramoto model in complex networks.  Physics Reports, 610:1–98, 2016.
  * [25] Y. Seginer.  The expected norm of random matrices.  Combinatorics, Probability and Computing, 9(2):149–166, 2000.
  * [26] S. Strogatz.  Sync: The emerging science of spontaneous order.  Penguin UK, 2004.
  * [27] A. Townsend, M. Stillman, and S. H. Strogatz.  Dense networks that do not synchronize and sparse ones that do.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(8):083142, 2020.
  * [28] R. Vershynin.  High-dimensional probability: An introduction with applications in data science, volume 47.  Cambridge university press, 2018.


# Cascading failures: dynamics, stability and control

Stefanny Ramirez, Maaike Odijk
Jan C. Willems Center for Systems and Control,
Engineering and Technology Institute Groningen,
Faculty of Science and Engineering, University of Groningen
Nijenborgh 4, 9747 AG Groningen, The Netherlands
<EMAIL_ADDRESS>
Dario Bauso
Jan C. Willems Center for Systems and Control,
Engineering and Technology Institute Groningen,
Faculty of Science and Engineering, University of Groningen
Nijenborgh 4, 9747 AG Groningen, The Netherlands
Dipartimento di Ingegneria, Università di Palermo
viale delle Scienze 90128 Palermo, Italy
<EMAIL_ADDRESS>

###### Abstract

We develop a dynamic model of cascading failures in a financial network
whereby cross-holdings are viewed as feedback, external assets investments as
inputs and failure penalties as static nonlinearities. We provide sufficient
milder and stronger conditions for the system to be a positive one, and study
equilibrium points and stability. Stability implies absence of cascades and
convergence of market values to constant values. We provide a constructive
method for control design to obtain stabilizing market investments in the form
of feedback-feedforward control inputs.

_K_ eywords Cascading failures $\cdot$ financial networks $\cdot$ positive
systems $\cdot$ stability analysis $\cdot$ systemic risks

##  1 Introduction.

Cascading failures in financial networks is the phenomenon by which a little
shock at one node propagates throughout the network. Nodes represent
organizations and edges represent interdependent cross-holding. The shocks are
financial failures and take place whenever the market value of an organization
falls below a threshold. The market value of an organization depends on their
investments on external assets and on shares of other organizations.

###  1.1 Main contribution.

We highlight the following main contributions. First, we extend the model in
[10] and [11] to a dynamic setting whereby market values change with time
based on the current investments and current cross-holdings. The model has a
feedback structure (cross-holdings) with external inputs (assets investments)
and static nonlinearities (discontinuous failure penalties).

Second, we provide a detailed analysis within the context of positive systems.
In particular, we provide sufficient conditions for the system to be a
positive one and provide a physical interpretation in terms of cascading
failure network resilience. Under such conditions, the system admits steady-
state solutions (equilibrium points).

Third, we analyze conditions for such equilibrium points to be asymptotically
stable by using the theory of Lyapunov stability. Stability implies absence of
cascades and convergence of market values to a constant value. Furthermore, we
relax the aforementioned sufficient conditions and prove that equilibrium
points may exist under milder conditions on the network connectivity.

Fourth, we extend the stability analysis to the case of large-network and
reinterpret the previous stability result in terms of centrality and degree
connectivity of the network.

Fifth, we provide a constructive method for control design to obtain
stabilizing market investments in the form of feedback-feedforward control
inputs. The feedforward control blocks the propagation of the failures while
the feedback control makes the market value dynamics stable and converging to
a stationary value.

###  1.2 Related Literature.

We build on the model in [11], and we extend it to a dynamic setting in the
same spirit as in [4]. Similar models as the one in [11] can be found in the
previous literature, see for instance [10]. In [10] the authors define the
equity value in terms of debts and payments between institutions, and the
cascade of failures is propagated when a failed institution reduces its
payment to its creditors, which can generate the failure of at least one of
these creditors. On the other hand, in [11] the equity value is defined in
terms of shares held by the organizations as well as by external shareholders,
and the propagation of failure is greatly affected by a failure cost.

Dynamic threshold models are a classic of the literature on cascading failures
and for further details we refer the reader to the seminal works [1], [14],
[16], [18], [24], [21] and [29]. Dynamical models to study and mitigate the
systemic risk in financial networks are developed in [5] and [7]. A dynamic
setting allows us to use tools from systems and control to examine the
stability of a system with multiple equilibrium points. In particular, our
model uses quaternion transformation matrices and this allows to view it as a
positive system [13], and analyse the transient of the system when there
exists any perturbation that drives it away from an equilibrium point.
Furthermore, thanks to this dynamic framework, in our work we design a
feedback-feedforward control to determine the market investments that
stabilize the system and stop the cascade of failures.

The literature on cascading failures is vast and grows in several directions
some of which depart from the scope of the current paper. For instance,
trading uncertainty and market frictions are the main focus in [9]. The
vulnerability to failures of interconnected economic networks, also referred
to as network of networks, is studied in [17]. Contagious links and their
long-term impact on the resilience of random networks are studied in [2].
Stochastic models of cascading failures to determine the moment of the first
failure in financial networks are developed in [27] and [19]. The underlying
perspective in this paper follows the idea of viewing “economists as
engineers” in the context of market design [25] and [26].

Previous attempts of assimilating cascading failures to epidemic models seem
to fall short of satisfactory answers when it comes to understanding pairwise
interactions, which is the main focus in [28]. The use of contagion models to
understand the dynamics of financial and biological networks is analyzed in
[23]. In [3] the authors adapt the framework of epidemiology models to default
propagation in financial networks. In this sense, neighbour correlation may
have an impact on preventing exponential growth of contagion [22]. Here the
focus is on the different behaviour observed in strong and weak social links.
More specifically, a main factor of limitation of contagion is provided by
strong links for which the neighbor relations are in general transitive
(friends of friends are usually friends) in accordance to the definition in
[16]. Another factor limiting contagion are random links between otherwise
distant nodes. This is observed for instance whenever the contagion requires
simultaneous exposure to multiple sources of activation [6]. Influence
maximization is also a challenging factor in networks with threshold models.
The idea here is to find the smallest set of nodes with maximal aggregated
influence [8], [15]. A game-theoretic perspective of contagion assimilates it
to a coordination game [12]. Cascades of failures follow the rules of
contagion in pairwise interactions between nodes in a network. This places
emphasis on certain network properties such as connectivity or centrality.
While connectivity allows the failure to propagate, a high connectivity is
proven to limit the diffusion due to the impact of high-degree nodes which are
very stable [20].

###  1.3 Organization of the manuscript.

This manuscript is organized as follows. In Section 2, we develop the model.
In Section 3, we cast the model within the theory of positive systems and
provide the equilibrium point analysis including stability. In Section 5, we
perform numerical simulations. Finally in Section 6 we provide conclusions.

##  2 Model.

The financial network involves a set of companies $N=\\{1,\ldots,n\\}$ and a
set of primitive assets $M=\\{1,\ldots,m\\}$. Asset $h\in M$ is characterized
by the market price $p_{h}\in\mathbb{R}$. Company $i\in N$ is characterized by
an equity value $V_{i}\in\mathbb{R}$ and a market value $v_{i}\in\mathbb{R}$.
Company $i$ owns a fraction $c_{ij}\in[0,1]$ of company $j$ for each pair
$i,j\in N$, $i\not=j$. No company has shares in itself, i.e., $c_{ii}=0$
$\forall i\in N$. The fraction of company $i$ which is not owned by other
companies is $\hat{c}_{ii}:=1-\sum_{j\in N,\,i\not=j}{c_{ji}}$,
$\hat{c}_{ii}>0$. Likewise, company $i$ has a share $D_{ih}$ of asset $h$, for
each $i\in N$ and $h\in M$. To introduce a compact notation, denote by
$V=[V_{i}]_{i\in N}\in\mathbb{R}^{n}$ and $v=[v_{i}]_{i\in
N}\in\mathbb{R}^{n}$ the vector of equity values and market values,
respectively. Also let $C:=[c_{il}]_{i,l\in N}\in\mathbb{R}^{n\times n}$ be
the matrix of cross-holding, $D:=[D_{ih}]_{i\in N,\,h\in M}$ be the matrix of
shares on primitive assets, and $p=[p_{h}]_{h\in M}$ the vector of market
prices for the primitive assets. Let also denote by
$\hat{C}:=diag(\hat{c}_{ii})\in\mathbb{R}^{n\times n}$ the diagonal matrix
containing all the fractions of companies not owned by other companies. When
the market value $v_{i}$ of company $i$ goes below a given threshold
$\underline{v}_{i}$ the company ‘fails’ and incurs a fixed nonnegative failure
cost $\beta_{i}\in\mathbb{R}_{+}$. The equity value is determined by the value
of the company primitive assets, the value of its shares on other companies
and the failure cost according to the following discrete-time dynamic model:

$\begin{cases}V(t+1)=CV(t)+Dp-\phi(v),\\\ v(t)=\hat{C}V(t),\end{cases}$ (1)

where $\phi(v)=[\beta_{i}I_{v_{i}<\underline{v}_{i}}(v_{i})]_{i\in
N}\in\mathbb{R}^{n}$ and $I_{v_{i}<\underline{v}_{i}}(v_{i})=1$ if
$v_{i}<\underline{v}_{i}$ and $I_{v_{i}<\underline{v}_{i}}(v_{i})=0$
otherwise.

To investigate the properties of (1) in the context of positive systems let us
rewrite it in terms of the new variables

$x_{i}(t)=v_{i}(t)-\underline{v}_{i},\quad\mbox{for all $i\in N$}.$ (2)

Company $i$ fails if $x_{i}(t)<0$, and does not fail if $x_{i}(t)\geq 0$.
Therefore the orthant that contains the current state $x(t)=[x_{i}(t)]_{i\in
N}$ determines also which companies have failed at time $t$. When an orthant
is invariant, namely when the state trajectory does not transition between
different orthants, the sign of $x_{i}(t)$ does not change over time and there
is no cascade of failure.

The term $\phi(v)$ can be modelled as a bang-bang input. Let us consider the
offset $U_{i}^{(offset)}=\frac{\beta}{2}$ and let us introduce the function

$U_{i}(x):=\begin{cases}\frac{\beta}{2},&\text{if}{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}...}x_{i}<0,\\\
-\frac{\beta}{2},&\text{if}{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>0.\end{cases}$ (3)

Let the following vectors be given $U(x):=[U_{i}(x)]_{i\in
N}\in\mathbb{R}^{n}$ and $U^{(offset)}:=[U_{i}^{(offset)}]_{i\in
N}\in\mathbb{R}^{n}$. We can rewrite $\phi(v)$ as follows:

$\phi(v)=U(x)+U^{(offset)}=U(x)+\mathbf{1}\frac{\beta}{2}.$ (4)

The term $\phi(v)$ as in (4) represents an offset of the penalty term and can
be interpreted as a static nonlinearity. From (1), (2) and (4), after denoting
$\underline{v}:=[\underline{v}_{i}]_{i\in N}$, dynamics (1) can be
reformulated as a dynamics on the error $x(t)$ as:

$\displaystyle\begin{array}[]{lll}x(t+1)&=v(t+1)-\underline{v}=\hat{C}V(t+1)-\underline{v}=\hat{C}\Big{[}CV(t)+Dp-\phi(v)\Big{]}-\underline{v}\\\
&=\hat{C}C\hat{C}^{-1}x(t)+(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}(Dp-U(x(t))-\mathbf{1}\frac{\beta}{2}).\end{array}$
(7)

In the last equality we have used the fact that

$V(t)=\hat{C}^{-1}v(t)=\hat{C}^{-1}(x(t)+\underline{v}).$

Let us enumerate the $2^{n}$ orthants according to a specific but random order
by using the index $r=1,\ldots,\hat{r}$ where $\hat{r}:=2^{n}$. Let us denote
by $\mathcal{O}^{(r)}$ the $r$th orthant. For each orthant $\mathcal{O}^{(r)}$
$r=1,\ldots,\hat{r}$ consider the following quaternion matrix
$J^{(r)}=[j_{il}^{(r)}]_{i,l\in N}$ where $j^{(r)}_{il}:=0$, for all $i,l\in
N,\,i\not=l,$ and for any generic $x$ in $\mathcal{O}^{(r)}$

$j^{(r)}_{ii}:=\begin{cases}-1,&{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}....}\text{if}{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}.}x_{i}<0,\\\
1,&{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}....}\text{if}{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>0.\end{cases}$ (8)

This transformation allows to turn any orthant into the positive orthant. In
other words, for any initial state $x(0)$ in the $r$th orthant
$\mathcal{O}^{(r)}$ we have a corresponding initial state in the positive
orthant given by

$X(0)=J^{(r)}x(0)\in\mathbb{R}_{+}^{n}.$ (9)

The transformation matrix can be used to define any generic orthant as follows
$\mathcal{O}^{(r)}:=\\{x\in\mathbb{R}^{n}|\,J^{(r)}x\geq 0\\}$ for all
$r=1,\ldots,\hat{r}$. Using the quaternion matrix $J^{(r)}$, dynamics (7) can
be written in terms of the new variable $X(t)$ as:

$\begin{array}[]{ll}X(t+1)=J^{(r)}\hat{C}C\hat{C}^{-1}J^{(r)}X(t)+J^{(r)}\Big{[}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}\Big{(}Dp-U(J^{(r)}X(t))-\mathbf{1}\frac{\beta}{2}\Big{)}\Big{]}.\end{array}$
(10)

After introducing the new matrix

$\begin{array}[]{lll}A^{(r)}:=[a_{ij}]_{i,j=1,\ldots,n}:=J^{(r)}\hat{C}C\hat{C}^{-1}J^{(r)}=\begin{bmatrix}0&\dots&j_{nn}\frac{1}{\hat{c}_{nn}}c_{1n}j_{11}\hat{c}_{11}\\\
j_{11}\frac{1}{\hat{c}_{11}}c_{21}j_{22}\hat{c}_{22}&\dots&j_{nn}\frac{1}{\hat{c}_{nn}}c_{2n}j_{22}\hat{c}_{22}\\\
\vdots&\ddots&\vdots\\\
j_{11}\frac{1}{\hat{c}_{11}}c_{n1}j_{nn}\hat{c}_{nn}&\dots&0\\\
\end{bmatrix},\end{array}$ (11)

and vector

$b^{(r)}(X):=J^{(r)}\Big{[}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}\Big{(}Dp-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2}\Big{)}\Big{]}$
(12)

the model can be written in compact form as:

$X(t+1)=A^{(r)}X(t)+b^{(r)}(X(t)).$ (13)

Equivalently we can write the above in input-output form as

$\begin{cases}X(t+1)=A^{(r)}X(t)+b^{(r)}(X(t)),\\\
v(t)=J^{(r)}X(t)+\underline{v}.\end{cases}$ (14)

Figure 1 depicts a block diagram representation of (14) in feedback-form.

$\displaystyle\begin{array}[]{ll}X(t+1)=A^{(r)}X(t)+b^{(r)}(X(t))\\\
v(t)=J^{(r)}X(t)+\underline{v}\end{array}$
$v(t)$$v(t)$$X(t)=J^{(r)}x(t)$$+$$-$$\underline{v}$$x(t)$$U(.)$$-$$-$$\frac{\beta}{2}$$-\phi(v)$$Dp$
Figure 1: Block diagram representation of (14) in feedback-form.

##  3 Analysis.

We recall that a square matrix $A\in\mathbb{R}^{n\times n}$ is nonnegative if
the entries $a_{ij}\geq 0$ for all $i,j=1,\ldots n$. We denote $A\geq 0$ and
$v\geq 0$ for a nonnegative matrix $A$ and vector $v$. If the inequalities are
strict, namely $a_{ij}>0$ then the matrix is said positive. We denote $A>0$
and $v>0$ for a positive matrix A and vector $v$.

An affine system of type $x(k+1)=Ax(k)+b$ is positive if $A\geq 0$ and $b\geq
0$. For a positive system it holds that $x(0)\geq 0$ implies $x(k)\geq 0$ for
all $k=1,2,\ldots$. Also, if $A$ is convergent then
$\lim_{k\rightarrow\infty}x(k)=(I-A)^{-1}b\geq 0$. We recall that $A$ is
convergent if all eigenvalues are in the unit disc, namely $|\lambda|<1$ where
$\lambda\in spec(A)$ and $spec(A)$ denotes the spectrum of $A$ (set of all
eigenvalues).

Consider the generic orthant $\mathcal{O}^{(r)}$ and let
$\mathcal{P}^{(r)}:=\\{i\in N:\,j_{ii}=1\\}$ and $\mathcal{N}^{(r)}:=\\{i\in
N:\,j_{ii}=-1\\}$. Evidently it holds
$N=\mathcal{P}^{(r)}\cup\mathcal{N}^{(r)}$ and $\mathcal{P}^{(r)}$,
$\mathcal{N}^{(r)}$ are a partition of $N$. Let us define
$b^{(r)}_{+}(X):=[b^{(r)}_{i}(X)]_{i\in\mathcal{P}^{(r)}}\in\mathbb{R}^{|\mathcal{P}^{(r)}|}$.
In other words, $b^{(r)}_{+}(X)$ is the vector obtained by extracting only the
components $i\in\mathcal{P}^{(r)}$ from $b^{(r)}(X)$.

The first result is about the sign of the vector function
$b^{(r)}(X)\in\mathbb{R}^{n}$ as defined in (12).

###### Lemma 1

Let $x\in\mathcal{O}^{(r)}$ be given, and consider $X=J^{(r)}x$. We have
$b^{(r)}(X)\geq 0$ if for all $i=1,\ldots,n$:

$\displaystyle(Dp)_{i}\leq\Big{(}(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}\Big{)}_{i}+\beta\qquad\text{
if }i\in\mathcal{N}^{(r)},$ (15)
$\displaystyle(Dp)_{i}\geq\Big{(}(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}\Big{)}_{i}\qquad\text{if
}i\in\mathcal{P}^{(r)}.$ (16)

###### Proof 1

Let us start by considering that for the $i$th component we have
$b^{(r)}_{i}(X)\geq 0$ if and only if

$\displaystyle\begin{array}[]{lll}b^{(r)}_{i}(X)&=\Big{(}J^{(r)}\Big{[}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}\Big{(}Dp-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2}\Big{)}\Big{]}\Big{)}_{i}\\\
&=j^{(r)}_{ii}\Big{(}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}\Big{(}Dp-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2}\Big{)}\Big{)}_{i}\geq
0.\end{array}$ (19)

For the first part of the proof let $j^{(r)}_{ii}=-1$. From $j^{(r)}_{ii}=-1$,
inequality (19) is true if and only if

$\Big{(}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}\Big{(}Dp-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2}\Big{)}\Big{)}_{i}\leq
0.$

From $\hat{C}$ being positive and diagonal the pre-multiplication of the
vector in the left-hand-side by $\hat{C}^{-1}$ keeps the sign of the elements
of the vector. Hence, the above inequality is true if and only if

$\begin{array}[]{cc}\Big{(}\hat{C}^{-1}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}^{-1}\hat{C}(Dp-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2})\Big{)}_{i}\\\
=\Big{(}(C\hat{C}^{-1}-\hat{C}^{-1})\underline{v}+Dp-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2}\Big{)}_{i}\leq
0.\end{array}$ (20)

From the above we obtain

$(Dp)_{i}\leq\Big{(}-(C\hat{C}^{-1}-\hat{C}^{-1})\underline{v}+U(J^{(r)}X)+\mathbf{1}\frac{\beta}{2}\Big{)}_{i}.$
(21)

From $j^{(r)}_{ii}=-1$ we have
$\Big{(}U(J^{(r)}X)+\mathbf{1}\frac{\beta}{2}\Big{)}_{i}=\beta$ and (15) is
proven.

For the second inequality (16), let $j^{(r)}_{ii}=1$. Then inequality (19) is
true if and only if

$\Big{(}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}\Big{(}Dp-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2}\Big{)}\Big{)}_{i}\geq
0.$

The above inequality is true if and only if

$\begin{array}[]{cc}\Big{(}\hat{C}^{-1}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}^{-1}\hat{C}\Big{(}Dp-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2}\Big{)}\Big{)}_{i}\\\
=\Big{(}(C\hat{C}^{-1}-\hat{C}^{-1})\underline{v}+Dp-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2}\Big{)}_{i}\geq
0.\end{array}$ (22)

As $j_{ii}^{(r)}=1$ then we have $U(J^{(r)}X)+\mathbf{1}\frac{\beta}{2}=0$ and
for the above we obtain

$(Dp)_{i}\geq\Big{(}-(C\hat{C}^{-1}-\hat{C}^{-1})\underline{v}\Big{)}_{i},$
(23)

and (16) is proven. $\square$

The next result addresses the case where not all companies are inter-
dependent. In terms of the interdependent graph, this means that such graph is
not completely connected.

In preparation to the next result, let
$\hat{A}^{(r)}:=[A_{i\bullet}^{(r)}]_{i\in\mathcal{P}^{(r)}}$. Furthermore,
let $\hat{A}^{(r)}_{+}:=[a_{ij}^{(r)}]_{i,j\in\mathcal{P}^{(r)}}$ and
$\hat{A}^{(r)}_{-}:=[a_{ij}^{(r)}]_{i\in\mathcal{P}^{(r)},\,j\in\mathcal{N}^{(r)}}$.
In other words we can partition
$\hat{A}^{(r)}=[\hat{A}^{(r)}_{+}|\hat{A}^{(r)}_{-}]$. Also, let
$X_{+}(t):=[X_{i}(t)]_{i\in\mathcal{P}^{(r)}}$ and
$X_{-}(t):=[X_{i}(t)]_{i\in\mathcal{N}^{(r)}}$. Evidently we can partition
$X(t):=[X_{+}(t)^{T}\,|X_{-}(t)^{T}]^{T}$.

###### Lemma 2

Matrix $A^{(r)}$ is a positive matrix, i.e.

$c_{il}=0,\mbox{for all $i,l\in N$: $j^{(r)}_{ii}\neq j^{(r)}_{ll}$ }.$ (24)

Furthermore the submatrix
$\hat{A}^{(r)}:=[A_{i\bullet}^{(r)}]_{i\in\mathcal{P}^{(r)}}$ is a positive
matrix if

$c_{il}=0,\mbox{for all $i,l\in N$: $i\in\mathcal{P}^{(r)}$,
$l\in\mathcal{N}^{(r)}$}.$ (25)

###### Proof 2

The entries of matrix $A^{(r)}$ are defined by
$a^{(r)}_{il}=j^{(r)}_{ii}\hat{c}_{ii}c_{il}\frac{1}{\hat{c}_{ll}}j^{(r)}_{ll}{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>i\neq l$. If for every $j^{(r)}_{ii},j^{(r)}_{ll}\in\\{-1,1\\}$ where
$j^{(r)}_{ii}\neq j^{(r)}_{ll}$ the entry of the cross-holdings matrix is
$c_{il}=0$, then
$a^{(r)}_{il}=j^{(r)}_{ii}\hat{c}_{ii}c_{il}\frac{1}{\hat{c}_{ll}}j^{(r)}_{ll}=-1\times\hat{c}_{ii}\times
0\times\frac{1}{\hat{c}_{ll}}\times 1=0\geq 0$. In conclusion, matrix
$A^{(r)}$ is a positive matrix for as long as $c_{il}=0$ for every
$j^{(r)}_{ii},j^{(r)}_{ll}\in\\{-1,1\\}$ where $j^{(r)}_{ii}\neq
j^{(r)}_{ll}$. To prove (25) we simply apply the same reasoning only to the
generic row $i$ where $j_{ii}=1$. $\square$

For any given orthant $\mathcal{O}^{(r)}$, define the subspace
$\mathcal{S}^{(r)}:=\\{\xi\in\mathbb{R}^{n}|\,\xi_{i}\geq 0,\forall
i\in\mathcal{P}^{(r)}\\}$. In simple terms, the subspace $\mathcal{S}^{(r)}$
is the set of points where all positive components of a vector in
$\mathcal{O}^{(r)}$ remain positive (but negative components can become
positive). In other words, for a given $\mathcal{O}^{(r)}$,
$\mathcal{S}^{(r)}=\cup_{k,\mathcal{P}^{(k)}\supset\mathcal{P}^{(r)}}\mathcal{O}^{(k)}$.

We are ready to establish conditions under which $\mathcal{S}^{(r)}$ is
positively invariant.

###### Theorem 1

Let dynamics (13) be given where $X(0)\in\mathbb{R}_{+}^{n}$ and
$x=J^{(r)}X(0)$ is in $\mathcal{O}^{(r)}$. If (15)-(16) and (24) hold true,
then we have

$X(0)\geq 0\Longrightarrow X(t)\geq 0,\quad\mbox{for all $t=0,1,2,\ldots$}$
(26)

Similarly, if (16) and (25) hold true, then we have for all
$i\in\mathcal{P}^{(r)}$

$X_{i}(0)\geq 0\Longrightarrow X_{i}(t)\geq 0,\quad\mbox{for all
$t=0,1,2,\ldots$}$ (27)

###### Proof 3

From Lemmas 1 and 2, if (15)-(16) and (24) hold true, then we have that
$A^{(r)}\geq 0$ and $b^{(r)}\geq 0$. This means that (13) is a positive affine
system and this in turn implies that if the system is initialized in the
positive orthant, the trajectory will remain confined in the positive orthant.

For the second part of the proof, we can write

$X_{+}(t+1)=A^{(r)}_{+}X_{+}(t)+A^{(r)}_{-}X_{-}(t)+b^{(r)}_{+}(X(t)).$

From (25) we have that $A^{(r)}_{-}$ is a null matrix and therefore the above
is equivalent to

$X_{+}(t+1)=A^{(r)}_{+}X_{+}(t)+b^{(r)}_{+}(X(t))\geq 0.$

From (16) we have $b^{(r)}_{+}(X(t))\geq 0$. It also holds $A^{(r)}_{+}$ and
therefore we obtain (27).

This concludes our proof. $\square$

Condition (26) means that failures do not propagate and companies who fail at
the beginning will not come back to the market. Condition (27) means that
failures do not propagate but companies who fail at the beginning can come
back to the market. In both cases $\mathcal{S}^{(r)}$ is positively invariant.

In the following result we establish conditions for stability and obtain the
unique equilibrium point in the positive orthant.

###### Theorem 2

Consider dynamics (13), and let (15)-(16) and (24) hold true. Also let

$0\leq\sum_{l=1,\\\ l\neq i}^{n}\hat{c}_{ii}c_{il}\frac{1}{\hat{c}_{ll}}<1.$
(28)

Then $X(t)$ converges to the following unique equilibrium point:

$\begin{array}[]{lll}X^{*}=(I-J^{(r)}\hat{C}C\hat{C}^{-1}J^{(r)})^{-1}\cdot\Big{[}J^{(r)}\Big{(}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}(Dp-U(J^{(r)}X^{*})-\mathbf{1}\frac{\beta}{2})\Big{)}\Big{]}.\end{array}$
(29)

###### Proof 4

First let us note that since (15)-(16) and (24) hold true from Theorem 1 we
have that system (13) is positive. This means that the state trajectory never
leaves the positive orthant, namely condition (26) holds true. Within each
single orthant dynamics (14) is linear and therefore there exists one unique
equilibrium point. To obtain such an equilibrium point, note that from (10) we
have that

$\begin{array}[]{lll}X(t+1)=J^{(r)}\hat{C}C\hat{C}^{-1}J^{(r)}X(t)+J^{(r)}\Big{[}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}(Dp-U(J^{(r)}X(t))-\mathbf{1}\frac{\beta}{2})\Big{]}.\end{array}$
(30)

At the equilibrium it must hold:

$X(t+1)=X(t)=X^{*}.$ (31)

Now by substituting (31) in (30), we have

$\begin{array}[]{lll}X^{*}=J^{(r)}\hat{C}C\hat{C}^{-1}J^{(r)}X^{*}+J^{(r)}\Big{[}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}(Dp-U(J^{(r)}X^{*})-\mathbf{1}\frac{\beta}{2})\Big{]}.\end{array}$
(32)

Taking all linear terms involving $X^{*}$ to the left-hand-side we obtain

$\begin{array}[]{lll}(I-J^{(r)}\hat{C}C\hat{C}^{-1}J^{(r)})X^{*}=J^{(r)}\Big{[}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}\Big{(}Dp-U(J^{(r)}X^{*})-\mathbf{1}\frac{\beta}{2}\Big{)}\Big{]}.\end{array}$
(33)

By pre-multiplying by $(I-J^{(r)}\hat{C}C\hat{C}^{-1}J^{(r)})^{-1}$ both left
and right-hand side we obtain the equilibrium point in 29.

Now it is left to prove that such an equilibrium is asymptotically stable,
namely that it holds $\lim_{t\rightarrow\infty}X(t)=X^{*}\geq 0$. To do this,
it suffices to show that the matrix
$A^{(r)}:=J^{(r)}\hat{C}C\hat{C}^{-1}J^{(r)}$ is convergent and we know that
this is true if all eigenvalues are in the unit disc, namely $|\lambda|<1$
where $\lambda\in spec(A)$.

The above condition is true if the following holds:

$0\leq\min_{i}\sum_{l=1,\\\ l\neq
i}^{n}|a^{(r)}_{il}|\leq\max_{i}\sum_{l=1,\\\ l\neq i}^{n}|a^{(r)}_{il}|\leq
1.$ (34)

The above implies that:

$0\leq\sum_{l=1,\\\ l\neq i}^{n}|a^{(r)}_{il}|<1.$ (35)

Now we know that
$|a^{(r)}_{ij}|=|j^{(r)}_{ii}\hat{c}_{ii}c_{il}\frac{1}{\hat{c}_{ll}}j^{(r)}_{ll}|=|\hat{c}_{ii}c_{il}\frac{1}{\hat{c}_{ll}}|{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>i\neq l$. Condition (35) is then equivalent to condition (28) and this
concludes the proof. $\square$

To understand condition (28) let us assume that company $i$ owns $90\%$ of
company $j$ and has no share in other companies. Assume also that company $j$
in turns does not invest in other companies, then $C_{ii}=0.1$ and
$\sum_{j=1}^{n}\frac{C_{ij}}{\hat{C}_{jj}}=0.9$. Then
$\hat{C}_{ii}\Big{(}\sum_{j=1}^{n}\frac{C_{ij}}{\hat{C}_{jj}}\Big{)}=0.09<1$.
From this example we can see that when there is a company that owns a large
share of another company which in turn does not invest much in other
companies, then the former company represents a stable node in the network.

Condition (25) implies that $A^{(r)}_{-}$ is a null matrix. We can relax this
condition and assume $A^{(r)}_{-}$ being nonnull. We can still prove positive
invariance in the case where $X_{-}(t)$ is bounded.

Consider the hyperbox
$H:=\\{X_{-}\in\mathbb{R}^{|\mathcal{N}^{(r)}|}|\,X_{l}\leq\bar{X}_{l},\,\forall\,l\in\mathcal{N}^{(r)}\\}\subset\mathbb{R}^{|\mathcal{N}^{(r)}|}$
where $\bar{X}_{l}$ are given upper bounds. Assume that $X_{-}(t)\in H$ for
all $t\geq 0$. Also denote
$\bar{X}_{-}:=[\bar{X}_{l}]_{l\in\mathcal{N}^{(r)}}$.

###### Theorem 3

Consider dynamics (13), and assume $X_{-}(t)\in H$ for all $t\geq 0$. If

$(Dp)_{i}\geq\Big{(}(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}\Big{)}_{i}-\Big{(}A^{(r)}_{-}\bar{X}_{-}\Big{)}_{i}\forall
i\in\mathcal{P}^{(r)}$ (36)

then for all $i\in\mathcal{P}^{(r)}$

$X_{i}(0)\geq 0\Longrightarrow X_{i}(t)\geq 0,\quad\mbox{for all
$t=0,1,2,\ldots$}.$ (37)

###### Proof 5

Let us consider dynamics

$X_{+}(t+1)=A^{(r)}_{+}X_{+}(t)+A^{(r)}_{-}X_{-}(t)+b^{(r)}_{+}(X(t)).$ (38)

We know that $A_{+}^{(r)}\geq 0$.

Our aim is now to show that

$A^{(r)}_{-}X_{-}(t)+b^{(r)}_{+}(X(t))\geq 0.$

To do this, we recall that for all $i\in\mathcal{P}^{(r)}$,

$b^{(r)}_{i}(X)\geq-\Big{(}A^{(r)}_{-}\bar{X}_{-}\Big{)}_{i}$

is true if and only if

$\begin{array}[]{r}\Big{(}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}\Big{(}Dp-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2}\Big{)}\Big{)}_{i}\geq-\Big{(}A^{(r)}_{-}\bar{X}_{-}\Big{)}_{i}.\end{array}$
(39)

The above inequality is true if and only if

$\displaystyle\begin{array}[]{cc}&\Big{(}\hat{C}^{-1}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}^{-1}\hat{C}\Big{(}Dp-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2}\Big{)}\Big{)}_{i}\\\
&=\Big{(}(C\hat{C}^{-1}-\hat{C}^{-1})\underline{v}+Dp-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2}\Big{)}_{i}\geq-\hat{C}^{-1}_{ii}\Big{(}A^{(r)}_{-}\bar{X}_{-}\Big{)}_{i}.\end{array}$
(42)

Now, for all $i\in\mathcal{P}^{(r)}$ then we have
$U(J^{(r)}X)+\mathbf{1}\frac{\beta}{2}=0$ and for the above we obtain

$(Dp)_{i}\geq\Big{(}-(C\hat{C}^{-1}-\hat{C}^{-1})\underline{v}\Big{)}_{i}-\hat{C}^{-1}_{ii}\Big{(}A^{(r)}_{-}\bar{X}_{-}\Big{)}_{i},$
(43)

which corresponds to (36).

To complete the proof it is left to notice that condition (36) also implies

$b^{(r)}_{+}(X(t))\geq-A^{(r)}_{-}\bar{X}_{-}.$ (44)

From (44) we also have

$A^{(r)}_{-}X_{-}(t)+b^{(r)}_{+}(X(t))\geq 0,$

and therefore dynamics (38) satisfies the positive invariance condition

$X_{+}(0)\geq 0\Longrightarrow X_{+}(t)\geq 0,\quad\mbox{for all
$t=0,1,2,\ldots$}.$ (45)

The above condition coincides with (37). This concludes our proof. $\square$

###  3.1 Large-scale network.

Theorem 3 has a nice counterpart result in the case of a large-scale network.
For the sake of simplicity, assume homogeneity in the shares, namely for a
given $0<\alpha<1$,

$a_{ij}=\left\\{\begin{array}[]{lll}\alpha&\mbox{with prob.
$\tilde{\Theta}_{i}\in[0,1]$}\\\ 0&\mbox{with prob.
$1-\tilde{\Theta}_{i}$}\end{array}\right.,\quad\forall
i\in\mathcal{P}^{(r)},\,j\in\mathcal{N}^{(r)}.$

Then we have

$A^{(r)}_{-}:=[a_{ij}^{(r)}\in\\{0,\alpha\\}]_{i\in\mathcal{P}^{(r)},\,j\in\mathcal{N}^{(r)}}\in\\{0,\alpha\\}^{|\mathcal{P}^{(r)}|\times|\mathcal{N}^{(r)}|}.$

Note that $\tilde{\Theta}_{i}n$ is the average connectivity (degree) of the
network.

Also assume that for all $j\in\mathcal{N}^{(r)}$ $\bar{X}_{j}=\bar{\xi}$. Then
we can obtain the maximum number of links (cross-sharing) to failed company
for which company $i$ is still safe, and this value is given by

$\hat{k}_{i}=\Big{\lfloor}\frac{b_{i}^{(r)}(X)}{\alpha\bar{\xi}}\Big{\rfloor}.$

It is worth mentioning that the above condition is obtained by taking into
account a worst case analysis, since it derives from (36) once we replace
$\Big{(}A^{(r)}_{-}\bar{X}_{-}\Big{)}_{i}=\tilde{n}\alpha\bar{\xi}$ and check
for the maximum $\tilde{n}$ for which (36) is satisfied.

The probability of a link to a failed node for each company
$i\in\mathcal{P}^{(r)}$ is

$\Theta=\frac{|\mathcal{N}^{(r)}|}{|\mathcal{P}^{(r)}|+|\mathcal{N}^{(r)}|}=\frac{|\mathcal{N}^{(r)}|}{n}.$

Let $\Delta_{i}$ be the degree of company $i$, which is given by

$\Delta_{i}=\tilde{\Theta}_{i}n.$

In other words, for large $n$ the probability of creating a link is equal to
the number of links divided by the number of companies, i.e.,
$\tilde{\Theta}_{i}=\frac{\Delta_{i}}{n}$. Let
$\binom{n}{k}=\frac{n!}{k!(n-k)!}$. Then we have

$\mbox{Pr
\Big{\\{}$\Big{(}A^{(r)}_{-}\bar{X}_{-}\Big{)}_{i}\leq\hat{k}_{i}\alpha\bar{\xi}$\Big{\\}}
}=\sum_{j=0}^{\hat{k}_{i}}\binom{\Delta_{i}}{j}\Theta^{j}_{i}(1-\Theta_{i})^{\Delta_{i}-j}.$

The above is obtained from the fact that the probability of having $j$ links
to failed companies follows the binomial distribution

$\mbox{Pr
\Big{\\{}$\Big{(}A^{(r)}_{-}\bar{X}_{-}\Big{)}_{i}=j\alpha\bar{\xi}$\Big{\\}}
}=\binom{\Delta_{i}}{j}\Theta^{j}(1-\Theta)^{\Delta_{i}-j}.$

The probability of having at most $\hat{k}$ links to failed companies is then
obtained by calculating the cumulative distribution function of the latter.

Define function

$\hat{F}_{n}(\tau)=\sum_{i=1}^{n}\mathbf{1}_{\\{\hat{k}_{i}<\tilde{\Theta}_{i}\tau\\}},\quad\forall\tau=0,\ldots,n$
(46)

where $\mathbf{1}_{\\{\hat{k}_{i}<\tilde{\Theta}_{i}\tau\\}}$ is the indicator
of event $\hat{k}_{i}<\tilde{\Theta}_{i}\tau$, namely

$\mathbf{1}_{\\{\hat{k}_{i}<\tilde{\Theta}_{i}\tau\\}}=\left\\{\begin{array}[]{ll}1&\mbox{if
$\hat{k}_{i}<\tilde{\Theta}_{i}\tau$},\\\
0&\mbox{otherwise.}\end{array}\right.$ (47)

For a population of failed companies of size $\tau$, $\tilde{\Theta}_{i}\tau$
is the expected number of links that a single company $i$ has to failed
neighbour companies. Then $\hat{F}_{n}(\tau)$ is the number of companies
characterized by a safety threshold $\hat{k}$ which is smaller than the
expected number of links to failed companies $\tilde{\Theta}_{i}\tau$. We can
obtain the estimate of an upper bound of the number of failed companies as
follows:

$\widehat{|\mathcal{N}^{(r)}|}:=\max_{\tau=0,\ldots,n}\\{\hat{F}_{n}(\tau)\geq\tau\\}.$
(48)

It is apparent that $\widehat{|\mathcal{N}^{(r)}|}$ grows with increasing
average connectivity $\tilde{\Theta}n$.

##  4 Control Design.

In this section, we build on previous results on conditions for system (13) to
be positive and stable, to design an external input to guarantee such
conditions. The external input can be viewed as a discrete-time control signal
and can be interpreted as an investment strategy on external assets. Namely,
we design a feasible investment on external assets that guarantee the system
to be positive and stable, which means that cascading failures are prevented.

###  4.1 Feedback-feedforward structure.

We recall that the physical interpretation of the control is in terms of the
investments in external assets $Dp$. The correct investments in external
assets ensures that the market values of the non-failing companies remain
above the threshold market values. We partition control $Dp$ into two parts:

  * i)

A constant feedforward part $u^{(1)}$ to ensure that the input vector
$b^{(r)}$ is non-negative,

  * ii)

A dynamical feedback part $u^{(2)}$ to ensure that state space matrix
$A^{(r)}$ is positive and satisfies the sufficient condition for stability
$0\leq\sum_{l=1}^{n}a_{il}^{(r)}<1$.

Then, we have

$\displaystyle Dp=u^{(1)}$ $\displaystyle+u^{(2)},$ (49)
$\displaystyle\text{Where:}{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}...}u^{(1)}$
$\displaystyle=\underline{u}^{(r)},$ (50) $\displaystyle u^{(2)}$
$\displaystyle=K^{(r)}X(t).$ (51)

It is important to mention that the design of the control signal in terms of
investments in external assets is performed at each time period $t$. From
(13), replacing the term on external investment $Dp$ by the feedback-
feedforward input as in (49)-(51) yields the following closed-loop dynamical
system:

$\left\\{\begin{array}[]{ll}X(t+1)=(A^{(r)}+J^{(r)}\hat{C}K^{(r)})X(t)+J^{(r)}((\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}(\underline{u}^{(r)}-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2}),\\\
v(t)=J^{(r)}X(t)+\underline{v}.\end{array}\right.$ (52)

Let $K^{(r)}:=J^{(r)}\tilde{K}^{(r)}$, then the feedback matrix is
$J^{(r)}\hat{C}J^{(r)}\tilde{K}^{(r)}=J^{(r)}\hat{C}K^{(r)}$ and the closed-
loop matrix $\tilde{A}^{(r)}$ is:

$\begin{array}[]{lll}\tilde{A}^{(r)}:=[\tilde{a}_{ij}]_{i,j=1,\ldots,n}=(A^{(r)}+\hat{C}\tilde{K}^{(r)})=\begin{bmatrix}\hat{c}_{11}\tilde{k}_{11}&\ldots&a_{1n}+\hat{c}_{11}\tilde{k}_{1n}\\\
\vdots&\ddots&\vdots\\\
a_{n1}+\hat{c}_{nn}\tilde{k}_{n1}&\ldots&\hat{c}_{nn}\tilde{k}_{nn}\\\
\end{bmatrix}.\end{array}$ (53)

Let us also introduce the vector

$\tilde{b}^{(r)}:=J^{(r)}\Big{(}(\hat{C}C\hat{C}^{-1}-I)\underline{v}+\hat{C}(\underline{u}^{(r)}-U(J^{(r)}X)-\mathbf{1}\frac{\beta}{2})\Big{)}.$
(54)

Note that $\tilde{b}^{(r)}$ in (54) differs from ${b}^{(r)}$ in (12) as the
external assets’ investment term $Dp$ is now replaced by the yet to be
designed control input $\underline{u}^{(r)}$. System (52) can be rewritten in
compact way as follows:

$\left\\{\begin{array}[]{ll}X(t+1)=\tilde{A}^{(r)}X(t)+\tilde{b}^{(r)},\\\
v(t)=J^{(r)}X(t)+\underline{v}.\end{array}\right.$ (55)

In the following of this manuscript, system (55) is referred to as the closed-
loop dynamical system and all results will relate to it. To visualize the
dynamics of the closed-loop dynamical system, a block diagram of the model is
represented in 2.

$\displaystyle\begin{array}[]{ll}X(t+1)=A^{(r)}X(t)+b^{(r)}(X(t))\\\
v(t)=J^{(r)}X(t)+\underline{v}\end{array}$
$v(t)$$v(t)$$X(t)=J^{(r)}x(t)$$+$$-$$\underline{v}$$x(t)$$U(.)$$-$$-$$\frac{\beta}{2}$$Dp$$K$$+$$+$$u^{(2)}$$u^{(1)}$
Figure 2: Feedback-feedforward control on assets’ investments.

###  4.2 Ensuring $\tilde{b}^{(r)}$ is non-negative.

In this section, the constant feedforward control
$u^{(1)}=\underline{u}^{(r)}$ is designed to ensure that $\tilde{b}^{(r)}$ is
non-negative.

###### Theorem 4

Consider system (55), and let us design

$\underline{u}^{(r)}=(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}+U(J^{(r)}X)+\mathbf{1}\frac{\beta}{2}+J^{(r)}\xi,\quad\xi>0,$
(56)

then vector $\tilde{b}^{(r)}$ is non-negative.

###### Proof 6

From Lemma 1 it is known that the following constraints should hold for
$\tilde{b}^{(r)}$ to be non-negative:

$\displaystyle\underline{u}_{i}^{(r)}\leq\Big{(}(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}\Big{)}_{i}+\beta\qquad\text{
if }i\in\mathcal{N}^{(r)},$ (57)
$\displaystyle\underline{u}_{i}^{(r)}\geq\Big{(}(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}\Big{)}_{i}\qquad\text{if
}i\in\mathcal{P}^{(r)}.$ (58)

Substituting (56) in (57) we have that for all $i\in\mathcal{N}^{(r)}$:

$\Big{(}(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}\Big{)}_{i}+\beta+j_{ii}^{(r)}\xi_{i}\leq\Big{(}(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}\Big{)}_{i}+\beta,$

which is satisfied as $j_{ii}^{(r)}\xi_{i}<0$.

Analogously, substituting (56) in (58) we have that for all
$i\in\mathcal{P}^{(r)}$:

$\Big{(}(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}\Big{)}_{i}+j_{ii}^{(r)}\xi_{i}\geq\Big{(}(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}\Big{)}_{i},$

which is satisfied as $j_{ii}^{(r)}\xi_{i}>0$. $\square$

The above designed control $u^{(1)}$ once substituted in system (55) yields
the following dynamics:

$\left\\{\begin{array}[]{ll}X(t+1)=(A^{(r)}+J^{(r)}\hat{C}K^{(r)})X(t)+J^{(r)}\hat{C}J^{(r)}\xi,\\\
v(t)=J^{(r)}X(t)+\underline{v},\end{array}\right.$ (59)

where $\tilde{b}^{(r)}=\hat{C}J^{(r)}x^{(r)}\geq 0$.

###### Theorem 5

Consider system (55), and assume $X_{-}(t)\in H$ for all $t\geq 0$. Let us
design $\forall i\in\mathcal{P}^{(r)}$

$\underline{u}^{(r)}_{i}=\Big{(}(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}\Big{)}_{i}-\hat{C}^{-1}_{ii}\Big{(}A^{(r)}_{-}\bar{X}_{-}\Big{)}_{i}+\xi_{i},\quad\xi_{i}>0,$
(60)

then for all $i\in\mathcal{P}^{(r)}$

$X_{i}(0)\geq 0\Longrightarrow X_{i}(t)\geq 0,\quad\mbox{for all
$t=0,1,2,\ldots$}.$ (61)

###### Proof 7

From Theorem 3, $\forall i\in\mathcal{P}^{(r)}$ we have that
$\tilde{b}^{(r)}_{i}$ is non-negative if:

$\displaystyle\underline{u}_{i}^{(r)}\geq\Big{(}(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}\Big{)}_{i}-\hat{C}^{-1}_{ii}\Big{(}A^{(r)}_{-}\bar{X}_{-}\Big{)}_{i}.$
(62)

Substituting (60) in (62) we have that for all $i\in\mathcal{P}^{(r)}$:

$\begin{array}[]{l}\Big{(}(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}\Big{)}_{i}-\hat{C}^{-1}_{ii}\Big{(}A^{(r)}_{-}\bar{X}_{-}\Big{)}_{i}+\xi_{i}\geq\Big{(}(\hat{C}^{-1}-C\hat{C}^{-1})\underline{v}\Big{)}_{i}-\hat{C}^{-1}_{ii}\Big{(}A^{(r)}_{-}\bar{X}_{-}\Big{)}_{i},\end{array}$
(63)

which is satisfied as $\xi_{i}>0$. $\square$

###  4.3 Ensuring $\tilde{A}^{(r)}$ positive and stable.

In this section, the dynamical feedback control $u^{(2)}$ is designed, to
ensure that $\tilde{A}^{(r)}$ is positive and guarantees that the sufficient
condition for stability in Theorem 2 holds. As previously stated, we assume
that $K^{(r)}:=J^{(r)}\tilde{K}^{(r)}$. Matrix $\tilde{K}^{(r)}$ is hereafter
referred to as the control input matrix. To find the entries of the control
input matrix, the following linear program is developed.

Let us introduce the following new notation

$\tilde{k}:=vec(\tilde{K})=[\tilde{K}_{\bullet 1}^{T}\ldots\tilde{K}_{\bullet
n}^{T}]^{T},\quad\delta:=\Big{[}\frac{a_{\bullet
1}}{\hat{c}_{11}}^{T}\ldots\frac{a_{\bullet n}}{\hat{c}_{nn}}^{T}\Big{]}^{T}.$

Note that $\tilde{k}$ is the vectorization of $\tilde{K}$ obtained by stacking
the columns of the matrix $\tilde{K}$ on top of one another. Also, let $e_{i}$
be the $i$-th canonical basis vector for the $n$-dimensional space, i.e.,
$e_{i}:=[0\ldots 0\,1\,0\ldots 0]^{T}$ and let us denote

$\mathbf{e}_{i}:=e_{i}\otimes\mathbf{1}_{n}=[\mathbf{0}^{T}\ldots\mathbf{0}^{T}\,\mathbf{1}^{T}_{n}\,\mathbf{0}^{T}\ldots\mathbf{0}^{T}]^{T}.$

Furthermore, for $\epsilon>0$, let us denote

$\gamma_{i}:=\frac{1}{\hat{c}_{ii}}(1-\displaystyle\sum\limits_{l=1}^{n}a_{il})-\epsilon,\quad\mu_{i}:=\frac{1}{\hat{c}_{ii}}\displaystyle\sum\limits_{l=1}^{n}a_{il}.$

The linear program is given by

(LP1)min $\displaystyle-\displaystyle\mathbf{1}^{T}_{n^{2}}\tilde{k}$ (64)
s.t. $\displaystyle-\tilde{k}\leq\delta$ (65)
$\displaystyle\displaystyle\mathbf{e}_{i}^{T}\tilde{k}\leq\gamma_{i}$
$\displaystyle\forall i=1,\ldots,n$ (66)
$\displaystyle-\mathbf{e}_{i}^{T}\tilde{k}\leq\mu_{i}$ $\displaystyle\forall
i=1,\ldots,n$ (67) $\displaystyle\tilde{k}_{il}=0$ $\displaystyle\forall i=l.$
(68)

###### Theorem 6

Consider system (55), and let $\tilde{K}^{(r)}$ be designed according to
(64)-(68), then matrix $\tilde{A}^{(r)}$ is positive, and satisfies the
sufficient conditions for stability in (28).

###### Proof 8

For the first part, we wish to show that matrix $\tilde{A}^{(r)}$ is a
positive matrix, i.e. $\tilde{A}^{(r)}>0$. This is true if all entries are
non-negative, i.e. $\tilde{a}_{il}\geq
0{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>i,l=1,...,n$. From (53) we know that the generic off-diagonal entry is
$\tilde{a}_{il}=a_{il}+\hat{c}_{ii}\tilde{k}_{il}$. Constraint (65) then
implies

$\tilde{a}_{il}=a_{il}+\hat{c}_{ii}\tilde{k}_{il}\geq 0.$

The above can be rewritten as

$\tilde{k}_{il}\leq\frac{a_{il}}{\hat{c}_{ii}}\quad\forall i,l=1,\ldots,n.$

This is the component-wise version of (65).

For the second part, we aim at showing that all eigenvalues of
$\tilde{A}^{(r)}$ are within the unit disc. It is well-known that all
eigenvalues of a matrix $\tilde{A}^{(r)}$ lay in the unit disc if all diagonal
values are zero and the inequalities
$0\leq\sum_{l=1}^{n}\tilde{a}_{il}<1{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>i=1,...,n$ hold. Constraint (68) ensures that all diagonal entries are null,
i.e.
$\tilde{a}_{il}=0{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>i=l$. Furthermore, constraints (66) and (67) ensure that each row sum of
matrix $\tilde{A}^{(r)}$ is smaller than one and not smaller than 0. These
constraints are obtained by setting $\sum_{l=1}^{n}\tilde{a}_{il}<1$ and
$\sum_{l=1}^{n}\tilde{a}_{il}\geq 0$. From (53) the latter can be equivalently
written as:

$\displaystyle\sum_{l=1}^{n}a_{il}+\sum_{l=1}^{n}\hat{c}_{ii}\tilde{k}_{il}<1,\quad\sum_{l=1}^{n}a_{il}+\sum_{l=1}^{n}\hat{c}_{ii}\tilde{k}_{il}\geq
0.$ (69)

Let us rewrite (69) by isolating $\tilde{k}_{il}$ to the left-hand-side. This
yields the following two sets of inequalities:

$\displaystyle\sum_{l=1}^{n}\tilde{k}_{il}<\frac{1-\sum_{l=1}^{n}a_{il}}{\hat{c}_{ii}},\quad\sum_{l=1}^{n}\tilde{k}_{il}\geq\frac{-\sum_{l=1}^{n}a_{il}}{\hat{c}_{ii}},$
(70)

which correspond to (66) and (67). This concludes the proof. $\square$

$k_{is}$$k_{ir}$$k_{ir}\geq-\frac{a_{ir}}{\hat{c}_{ii}}$$k_{il}\geq-\frac{a_{il}}{\hat{c}_{ii}}$$\displaystyle\sum\limits_{l=1}^{n}k_{il}\leq\gamma_{i}$$\displaystyle\sum\limits_{l=1}^{n}k_{il}\geq\mu_{i}$
Figure 3: Graphical interpretation of the linear program (65)-(68).

Figure 3 provides a graphical interpretation of the linear constraints
(65)-(68) projected on a plane $k_{is}-k_{ir}$. The set of feasible solutions
is in red. Constraints (65) are horizontal and vertical lines (dashed) which
correspond to inequalities $k_{is}\geq-\frac{a_{is}}{\hat{c}_{ii}}$ and
$k_{ir}\geq-\frac{a_{ir}}{\hat{c}_{ii}}$ respectively. The economic
interpretation is as follows. The horizontal and vertical lines (dashed)
identifies a minimum cash inflow from external investment to compensate the
drop in market value $i$ due to companies $s$ and $r$ failing. Constraint (66)
is represented by a diagonal line which corresponds to hyperplane $\sum
k_{il}=\gamma_{i}$ and identifies the feasible halfspace $\sum
k_{il}\leq\gamma_{i}$ (points below the diagonal line). Similarly, Constraint
(67) is represented by a diagonal line which corresponds to hyperplane $\sum
k_{il}=\mu_{i}$ and identifies the feasible halfspace $\sum k_{il}\geq\mu_{i}$
(points above the diagonal line).

###  4.4 Finding the investment choices $D$.

Once both the constant feedforward and dynamical feedback parts of the control
are designed, we need to translate such control into investment choices
through an opportune design of matrix $D$. This can be achieved by solving the
following linear program, which has the structure of a typical network flow
problem.

Let us define the new notation

$\begin{array}[]{lll}d:=vec(D)=[D_{\bullet 1}^{T}\ldots D_{\bullet
m}^{T}]^{T},\\\
\mathbf{\hat{e}_{i}}:=e_{i}\otimes\mathbf{1}_{m}=[\mathbf{0}^{T}\ldots\mathbf{0}^{T}\,\mathbf{1}^{T}_{m}\,\mathbf{0}^{T}\ldots\mathbf{0}^{T}]^{T},\\\
\mathbf{\tilde{e}_{i}}:=\mathbf{1}_{m}\otimes e_{i}=[e_{i}^{T}\ldots
e_{i}^{T}]^{T},\\\ \mathbf{\theta_{i}}:=[p_{1}e_{i}^{T}\ldots
p_{m}e_{i}^{T}]^{T}.\end{array}$

The linear program is then given by

(LP2)min $\displaystyle-\displaystyle\mathbf{1}_{nm}^{T}d$ (71) s.t.
$\displaystyle\displaystyle\mathbf{\theta_{i}}^{T}d=u^{(1)}_{i}+u^{(2)}_{i}$
$\displaystyle\forall i=1,\ldots,n$ (72)
$\displaystyle\displaystyle\mathbf{\tilde{e}_{i}}^{T}d\leq 1$
$\displaystyle\forall i=1,\ldots,n$ (73) $\displaystyle-d\leq 0$ (74)
$\displaystyle\displaystyle\mathbf{\hat{e}_{l}}^{T}d\leq 1$
$\displaystyle\forall l=1,\ldots,m$ (75)

The linear equality constraints (72) receive $u^{(1)}$ and $u^{(2)}$ as input
and return the entries $d_{il}$, $i,l=1,\ldots,n$. Furthermore, (72) works
only if the entries of $u^{(1)}_{i}+u^{(2)}_{i}$ are non-negative. However,
there exist possibilities that $u^{(1)}_{i}+u^{(2)}_{i}$ has negative entries.
To tackle this problem, the different cases of negative entries are discussed
below:

  * i)

If a company $i$ has failed, its future investments in external assets do not
matter anymore. Therefore, for a failed company $i$, entry
$u^{(1)}_{i}+u^{(2)}_{i}$ should always be set to zero.

  * ii)

If a company $i$ has not failed but its entry $u^{(1)}_{i}+u^{(2)}_{i}$ is
negative, this entry should be set to zero as well, because this company
requires no further investments to survive.

Therefore, all negative entries $u^{(1)}_{i}+u^{(2)}_{i}<0$ should be set to
zero prior to running the above linear program.

Consider the designed control $Dp$, where $D$ denotes the investment choices.
Then the linear program (71)-(75) finds the investment choices for all
companies. To see this, note that constraint (72) ensures that the choice of
investments in external assets for company $i$ results in the desired market
value to survive, as the required investments in external assets for the non-
failing companies are defined by $u^{(1)}_{i}+u^{(2)}_{i}$. Constraint (73)
ensures that a company can invest for no more than 100% to ensure that
companies do not spend more than possible. Constraint (74) ensures that only
positive investments are made since it is not possible to negatively invest in
external assets. Constraint (75) ensures that companies can never invest in
more than what is left of a certain external asset. Therefore, they have to
strategically decide in which asset they want to invest.

$p_{1}$$\sum\limits_{i=1}^{n}p_{1}d_{i1}\leq
p_{1}$$p_{2}$$\sum\limits_{i=1}^{n}p_{2}d_{i2}\leq
p_{2}$$\vdots$$p_{m}$$\sum\limits_{i=1}^{n}p_{m}d_{im}\leq
p_{m}$$1$$w_{1}=u^{(1)}_{1}+u^{(2)}_{1}$$2$$w_{2}=u^{(1)}_{2}+u^{(2)}_{2}$$\vdots$$n$$w_{n}=u^{(1)}_{n}+u^{(2)}_{n}$
Figure 4: Network flow problem representation of linear program (71)-(75).

Figure 4 provides a network flow problem representation of linear program
(71)-(75). The network topology is the one of a bipartite graph where the
nodes on the left represents assets and the nodes on the right represent
companies. We observe capacity constraints of type
$\sum\limits_{i=1}^{n}p_{l}d_{il}\leq p_{l}$ for each asset $l=1,\ldots,m$.
Such constraints are equivalent to (75). Similarly we have demand constraints
of type $\sum\limits_{l=1}^{m}p_{l}d_{il}=w_{i}=u^{(1)}_{i}+u^{(2)}_{i}$ for
each company $i=1,\ldots,n$. Such constraints are (72).

##  5 Numerical analysis.

Let us simulate for a set of one-hundred companies. All parameters are listed
in Table 1. In particular, we set the initial states within a range between
$1$ and $X(0)^{+}:=\max_{i=1,\ldots,n}X_{i}(0)=5000$. We set the threshold
$\underline{v}=100$ namely $2\%$ of $X(0)^{+}$. We set the incurred failure
cost $\beta=X(0)^{+}:=5000$. We set the asset values within the range from $1$
to about $0.1$ of $X(0)^{+}$. Each company invests on a different external
asset (namely matrix $D=\mathbb{I}_{(100\times 100)}$) and the values of the
assets are increasing with step size equal to $6$. By doing this, we obtain a
first set of (seventeen) companies for which condition (16) is violated and a
second set of (eighty-three) companies for which the same condition holds
true. All simulations are carried out with MATLAB on an Intel(R) Core(TM)2
Duo, CPU P8400 at 2.27 GHz and a 3GB of RAM. The horizon window consists of
$T=300$ iterations.

Parameters | Value
---|---
Number of companies | $n=100$
Initial states | $X(0)=\begin{bmatrix}1&5000\end{bmatrix}^{50}$
Threshold market value | $\underline{v}=100$
Incurred failure costs | $\beta=5000$
Assets values | $p=\begin{bmatrix}1&7&13&\ldots&595\end{bmatrix}^{\prime}$
Investments on assets | $D=\mathbb{I}_{(100\times 100)}$
Table 1: Parameters’ values

We build a cross-holding matrix as follows: $C_{ij}=\frac{1}{10n}$ with
probability 20% and zero otherwise. In other words, each company has a share
of 0.001% in two out of ten other companies. Note that this is the adjacency
matrix of a random network with average degree equal to 20. In other words,
all companies have approximately the same degree equal to 20. For each
company, we simulate the sample error according to dynamics (13). The pseudo-
code of the simulation algorithm is in Table 2.

Input: Set of parameters as in Table 1
---
Output: Time-evolution of error $X(t)$
$\quad 1:$ Initialize. Generate random matrix $C$,
initial values $X(0)$, assets’ values $p$
$\quad 2:$ for time $t=1,\ldots,T$ do
$\quad 3:\quad$ Calculate stabilizing control $u^{(1)}$ as in (60)
$\quad 4:\quad$ Calculate control $u^{(2)}$ solving LP (64)-(68)
$\quad 5:\quad$ Calculate error according to dynamics (13);
$\quad 6:$ end for
$\quad 7:$ Display cluster of nodes violating (16) to the left
and all other nodes to the right
$\quad 8:$ Color red failed companies and blue all others
$\quad 9:$ STOP
Table 2: Simulation algorithm

The outputs of the simulations are visualized in Figs. 5-7. In particular,
Fig. 5 displays a slowed down (step size $dt=0.1$) time plot of the evolution
of the error $X(t)$ for $t=1,\ldots,T$. Seventeen companies fail in the first
time instances and eventually four more companies fail. So in total we observe
twenty-one failed companies at the end of the horizon (stead-state). In Fig.
6, companies (nodes) for which condition (16) is violated are clusterized to
the left and all other nodes to the right. The size of the node increases with
the investment on external assets (bigger nodes invest more on external
assets).

Figure 5: Time plot of the evolution of the error $X(t)$ obtained from (13).
Threshold is indicated with a dashed line.

Simulations are initialized in the orthant where all companies are initially
healthy (no failures), see plot at time $t=0$ (top-left). However, as
condition (16) does not hold true for a subset of companies (vector
$b^{(r)}(X)$ is not positive), then we have that (26) does not hold true as
well, and the trajectory of the error $X(t)$ leaves the orthant. Some
companies start failing and failures propagate throughout the network quickly,
see plots at time $t=30$ (top-right) and $t=60$ (bottom-left). Failure
involves all companies for which (16) does not hold, but eventually also
companies for which (16) holds true. This occurs as condition (25) does not
hold, namely we have healthy companies from the right cluster which are still
connected with failed companies from the left cluster. In other words, it is
not true that at time $t=60$, $c_{il}=0$, for all $i,l\in N$,
$i\in\mathcal{P}^{(r)}$, $l\in\mathcal{N}^{(r)}$. Therefore also (27) does not
hold and three healthy companies eventually fail as well. You can think of
$\mathcal{P}^{(r)}$ at time $t=60$ as the cluster to the right and the fact
that three companies in this cluster eventually fail implies that (27) does
not hold. In other words, the set $\mathcal{S}^{(r)}$ at time $t=60$ is not
positively invariant. Note that eventually $20\%$ of the companies fail, which
is in agreement with the fact that each company has on average 20 neighbors
and $\frac{1}{5}$ of such neighbors is on average failed, namely each company
has on average 4 failed neighbor companies. If we look at the values of
$\hat{k}$ we notice that the first 20 companies have $\hat{k}\leq 4$. In Fig.
7 we plot function $\hat{F}_{n}(\tau)$ as in (46) (thick line). Estimate (48)
is the intersection point with bisect (thin line) which amounts to a bit less
the 20% of the total.

Figure 6: Cascading failure for a random network with average connectivity
$\langle\kappa\rangle=20$: failed companies in red at $t=0,30,60,T$. Figure 7:
Function $\hat{F}_{n}(\tau)$ as in (46) (thick line). Estimate (48) is the
intersection point with bisect (thin line).

We repeat the experiment with a dense graph and the result is displayed in
Figs. 8-10. Here we increase the probability of cross-holding from 20% to 80%.
This means that we have a random network with an average degree equal to 80.
Each company has cross-shares to about 80% of the other companies. In
particular, Fig. 8 displays a slowed down (step size $dt=0.1$) time plot of
the evolution of the error $X(t)$ for $t=1,\ldots,T$. Eventually forty-eight
companies fail at the end of the horizon (stead-state). This amounts to about
$50\%$ of the companies, which is in agreement with the fact that each company
has on average 80 neighbors and $\frac{1}{2}$ of such neighbors is on average
failed, namely each company has on average 40 failed neighbor companies. If we
look at the values of $\hat{k}$ we notice that the first 50 companies have
$\hat{k}\leq 40$. In Fig. 10 we plot function $\hat{F}_{n}(\tau)$ as in (46)
(thick line). Estimate (48) is the intersection point with bisect (thin line)
which amounts to about 50% of the total.

Figure 8: Time plot of the evolution of the error $X(t)$ obtained from (13).
Figure 9: Cascading failure for a random network with average connectivity
$\langle\kappa\rangle=80$: failed companies in red at $t=0,30,60,T$. Figure
10: Function $\hat{F}_{n}(\tau)$ as in (46) (thick line). Estimate (48) is the
intersection point with bisect (thin line).

We repeat the experiment a third time with a sparse graph and the result is
displayed in Figs. 11-13. Here the probability of cross-holding follows a
power law $P(\kappa)\sim\kappa^{-\rho}$, where $P(\kappa)$ is the fraction of
nodes in the network having $\kappa$ out-degree connections to other nodes and
$\rho=2.1$. This distribution generates a scale-free networks with very many
nodes poorly connected and very few hubs. The average connectivity is
$\langle\kappa\rangle=20$. This time the set $\mathcal{S}^{(r)}$ at time
$t=60$ is positively invariant (no healthy companies from the right cluster
fail). In particular, Fig. 11 displays a slowed down (step size $dt=0.1$) time
plot of the evolution of the error $X(t)$ for $t=1,\ldots,T$. Eventually
seventeen companies fail at the end of the horizon (stead-state). Note that
eventually about $20\%$ of the companies fail, which is in agreement with the
fact that each company has on average 20 neighbors and $\frac{1}{5}$ of such
neighbors is on average failed, namely each company has on average 4 failed
neighbor companies. In Fig. 13 we plot function $\hat{F}_{n}(\tau)$ as in (46)
(thick line). Estimate (48) is the intersection point with bisect (thin line)
which amounts to a bit less the 20% of the total.

Figure 11: Time plot of the evolution of the error $X(t)$ obtained from (13).
Figure 12: Cascading failure for a scale-free network with average
connectivity $\langle\kappa\rangle=2$: failed companies in red at
$t=0,30,60,T$. Figure 13: Function $\hat{F}_{n}(\tau)$ as in (46) (thick
line). Estimate (48) is the intersection point with bisect (thin line).

In a fourth set of simulations we consider again a random network with an
average degree $\langle\kappa\rangle=80$. Figure 14 displays the time plot of
the error (left). Failures propagate through the whole network. Figure 14
(right) shows the effect of the stabilizing control $u^{(1)}$ as in (60). The
control becomes active at time $t=60$ thus preventing further cascade of
failures. The y-axis on the right shows the magnitude of $u^{(1)}$, which can
be viewed as cash flow deriving from external assets.

---
Figure 14: Left: Failures propagate through the whole network. Right:
Stabilizing control $u^{(1)}$ (60) activates at time $t=60$.

In a fifth set of simulations we consider again a random network with an
average degree $\langle\kappa\rangle\approx 20$. Figure 15 displays the time
plot of the error (left). Dynamics is very slow and close to instability and
indeed we obtain eigenvalues very close to 1 in absolute value. Figure 15
(right) shows the effect of the stabilizing control $u^{(2)}$ where
$\hat{K}^{(r)}$ is designed according to the linear program (64)-(68).

---
Figure 15: Left: Dynamics is very slow and close to instability. Right: effect
of the stabilizing control $u^{(2)}$ where $\hat{K}^{(r)}$ is from (65)-(68).

##  6 Conclusions.

For a financial network we have developed a dynamic model of cascading
failures. Theoretical analysis has provided a deeper understanding of the
conditions that prevent cascades to occur and lead to stability of the market
values around certain equilibrium points. Such conditions involve the
structure of the cross-holding matrix as well as the magnitude of the external
assets. Results are illustrated on a large network with different
connectivity. Finally, we have also provided a control design method to
optimize the market investments which are now viewed as feedback-feedforward
control inputs.

## Acknowledgments.

We thank Prof. Dr. Rakesh Vohra (University of Pennsylvania), for the useful
discussions. In particular we thank him for his insights to improve the
explanation of the motivation and novelty of our work, for the inputs to
clarify the assumptions and conditions of our model, and for his comments
about the control design.

## References

  * [1] Adam, E. M. and Dahleh, M. A. and Ozdaglar, A. (2012). On the behavior of threshold models over finite networks. Proceedings of the IEEE 51st Annual Conference on Decision and Control (CDC), 2672–2677.
  * [2] Amini, H. and Cont, R. and Minca, A. (2016). Resilience to contagion in financial networks. Mathematical Finance, 26(2):329–365.
  * [3] Barja, A. and Martínez, A. and Arenas, A. and Fleurquin, P. and Nin, J. and Ramasco, J. and Tomás E. (2019). Assesing the risk of default propagation in interconnected sectoreal financial networks. EPJ Data Science, 8–32.
  * [4] Calafiore, G. C. and Fracastoro, G. and Proskurnikov, A. V. (2022). Control of Dynamic Financial Networks. IEEE Control Systems Letters, 3206–3211.
  * [5] Capponi, A. and Chen, P. C. (2015). Systemic risk mitigation in financial networks. Journal of Economic Dynamics and Control, 58:3206–3211.
  * [6] Centola, D. and Eguíluz, V. M. and Macy, M. W. (2007). Cascade dynamics of complex propagation. Phisica, A(374):449–456.
  * [7] Chen, H. and Wang, T. and Yao, D. (2021). Financial network and systemic risk-a dynamic model. Production and Operations Management, 30(8):2441–2466.
  * [8] Chen, W. and Yuan, Y. and Zhang, L. (2010). Scalable influence maximization in social networks under the linear threshold model. Proceedings of the 10th IEEE International Conference on Data Mining (ICDM), 88–97.
  * [9] Douglas, M. G. and Shachar, K. (2007). Financial Networks. The American Economic Review, 97(2):99–103.
  * [10] Eisenberg, L. and Noe, T. H. (2001) Systemic risk in financial systems. Management Science, 47(2):236–249.
  * [11] Elliott, M. and Golub, B. and Jackson, M. O. (2014) Financial networks and contagion. American Economic Review, 104(10):3115–53.
  * [12] Ellison, G. (1993) Learning, local interaction, and coordination. Econometrica, 61(5):1047–1071.
  * [13] Farina, L. and Rinaldi, S. (2000) Positive linear systems: theory and applications, volume 50 (John Wiley & Sons).
  * [14] Glasserman, P. and Young, H. P. (2016) Contagion in financial networks. Journal of Economic Literature 54(3):779–831.
  * [15] Goyal, A. and Lu, W. and Lakshmanan, L. V. S. (2011) Simpath: An efficient algorithm for influence maximization under the linear threshold model. Proceedings of the 11th IEEE International Conference on Data Mining (ICDM) 211–220.
  * [16] Granovetter, M. (1978) Threshold models of collective behavior. American Journal of Sociology 1420–1443.
  * [17] Havlin, S. and Kenett, D. Y. (2015) Cascading failures in interdependent economic networks. Proceedings of the International Conference on Social Modeling and Simulation, plus Econophysics Colloquium 2014 87–97 (Springer, Cham).
  * [18] Jackson, M. O. (2008) Social and economic networks. Princeton University Press, NJ, USA.
  * [19] Lee, H. (2020) On a new stochastic model for cascading failures. Journal of Applied Probability 57(4):1150–1161.
  * [20] Lelarge, M. (2012) Diffusion and cascading behavior in random networks. Games and Economic Behavior 75(2):752–775.
  * [21] Massai, L. and Como, G. and Fagnani, F. (2021) Equilibria and systemic risk in saturated networks. Mathematics of Operation Research online, as https://doi.org/10.1287/moor.2021.1188.
  * [22] Morris, S. (2000) Contagion. The Review of Economic Studies 67(1):57–78.
  * [23] Peckham, R. (2013) Economies of contagion: financial crisis and pandemic. Economy and Society 42(2):226–248.
  * [24] Rossi, W. S. and Como, G. and Fagnani, F. (2017) Threshold models of cascades in large-scale networks. IEEE Transactions on Network Science and Engineeringy 6(2):158–172.
  * [25] Roth, A. (2002) The Economist as Engineer: Game Theory, Experimentation, and Computation as Tools for Design Economics. Econometrica 70(4):1341–1378.
  * [26] Smith, N. (2020) Economists are starting to think like engineers. experimentation, and computation as tools for design economics. https://www.bloomberg.com/opinion/articles/2020-10-12/economists-are-starting-to-think-like-engineers?sref=R8NfLgwS.
  * [27] Swift, A. W. (2008) Stochastic models of cascading failures, Experimentation, and Computation as Tools for Design Economics. Journal of Applied Probability 45:907–921.
  * [28] Vega-Redondo, F. (2007) Complex social networks, Ch. 4. (Econometric Society Monographs, Cambridge University Press).
  * [29] Victor, A. and Santosh, V. and Rakesh, V. (2021) Contagion and Equilibria in Diversified Financial Networks. PIER Working Paper 21(29).
  * [30] Victor, A. and Santosh, V. and Rakesh, V. (1998) Individual strategy and social structure: an evolutionary theory of institutions, Ch. 3. (Princeton University Press).