# Causal Rate Distortion Function on Abstract Alphabets and Optimal
Reconstruction Kernel

Charalambos D. Charalambous, Photios A. Stavrou and Christos K. Kourtellaris
C. D. Charalambous (chadcha@ucy.ac.cy). P. A. Stavrou
(stavrou.fotios@ucy.ac.cy). C. K. Kourtellaris
(kourtellaris.christos@ucy.ac.cy). The authors are with the Department of
Electrical and Computer Engineering, University of Cyprus, Nicosia, CYPRUS

###### Abstract

A Causal rate distortion function with a general fidelity criterion is
formulated on abstract alphabets and the optimal reconstruction kernel is
derived, which consists of a product of causal kernels. In the process,
general abstract spaces are introduced to show existence of the minimizing
kernel using weak∗-convergence. Certain properties of the causal rate
distortion function are presented.

##  I INTRODUCTION

This paper is concerned with lossy data compression subject to distortion or
fidelity criterion and causal decoding on abstract alphabets. Its information
theoretic interpretation is the causal rate distortion function formulated via
the directed information between the source sequence
$X^{n}\stackrel{{\scriptstyle\triangle}}{{=}}\\{X_{0},X_{1},\ldots,X_{n}\\}$
and its reproduction sequence
$Y^{n}\stackrel{{\scriptstyle\triangle}}{{=}}\\{Y_{0},Y_{1},\ldots,Y_{n}\\}$
defined by

$\displaystyle I(X^{n}{\rightarrow}Y^{n})$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\sum_{i=0}^{n}I(X^{i};Y_{i}|Y^{i-1})$
(1)

The average distortion constraint is

$\displaystyle
E\\{d_{0,n}(X^{n},Y^{n})\\}\leq{D},\>d_{0,n}(x^{n},y^{n})\stackrel{{\scriptstyle\triangle}}{{=}}\sum^{n}_{i=0}\rho_{0,i}(x^{i},y^{i})$
(2)

where $D\geq 0$, $d_{0,n}(\cdot,\cdot)$ a non-negative distortion function.
Define the causal product of conditional distributions by

$\displaystyle{\overrightarrow{P}}_{Y^{n}|X^{n}}(dy^{n}|x^{n})\stackrel{{\scriptstyle\triangle}}{{=}}\otimes^{n}_{i=0}P_{Y_{i}|Y^{i-1},X^{i}}(dy_{i}|y^{i-1},x^{i})$
(3)

where $P_{Y_{i}|Y^{i-1},X^{i}}(dy_{i}|y^{i-1},x^{i})$ denotes the conditional
distribution of $Y_{i}$ given $(Y^{i-1},X^{i}),~{}i=0,1,\ldots,n.$
Since causal codes as defined in [4] satisfy
$P_{X_{i}|X^{i-1},Y^{i-1}}(dx_{i}|x^{i-1},y^{i-1})=P_{X_{i}|X^{i-1}}(dx_{i}|x^{i-1})$.
$P-a.s$ (see also Lemma II.4), in the analysis it is convenient to express
$I(X^{n}\rightarrow{Y^{n}})$ as a functional of
${\overrightarrow{P}}_{Y^{n}|X^{n}}(dy^{n}|x^{n})$ as follows.

$\displaystyle I(X^{n}\rightarrow{Y^{n}})$
$\displaystyle=\int\log(\frac{{\overrightarrow{P}}_{Y^{n}|X^{n}}(dy^{n}|x^{n})}{{P}_{Y^{n}}(dy^{n})})$
$\displaystyle\times{\overrightarrow{P}}_{Y^{n}|X^{n}}(dy^{n}|x^{n})P_{X^{n}}(dx^{n})$
(4) $\displaystyle={\mathbb{I}}(P_{X^{n}},{\overrightarrow{P}}_{Y^{n}|X^{n}})$
(5)

where ${\mathbb{I}}(P_{X^{n}},{\overrightarrow{P}}_{Y^{n}|X^{n}})$ indicates
the functional dependence of $I(X^{n}\rightarrow{Y^{n}})$ on
$\\{P_{X^{n}},{\overrightarrow{P}}_{Y^{n}|X^{n}}\\}$.
The causal information rate distortion function investigated is

$\displaystyle\inf_{{\overrightarrow{P}}_{Y^{n}|X^{n}}(dy^{n}|x^{n}):E\big{\\{}d_{0,n}(X^{n},Y^{n})\big{\\}}\leq
D}I(X^{n}\rightarrow{Y^{n}})$ (6)

Under appropriate assumptions on $d_{0,n}(\cdot,\cdot)$ it is shown that the
optimal causal product (reproduction channel)
${\overrightarrow{P}}^{*}_{Y^{n}|X^{n}}$ which achieves the infimum in (6) is
given by

$\displaystyle{\overrightarrow{P}}^{*}_{Y^{n}|X^{n}}(dy^{n}|x^{n})=\otimes^{n}_{i=0}\frac{e^{s\rho_{i}(x^{i},y^{i})}P^{*}_{Y_{i}|Y_{i-1}}(dy_{i}|y^{i-1})}{\int_{{\cal
Y}_{i}}e^{\rho_{i}(x^{i},y^{i})}P^{*}_{Y_{i}|Y_{i-1}}(dy_{i}|y^{i-1})}$ (7)

where $s\leq 0$ is the Lagrange multiplier associated with the fidelity
constraint. The operational meaning of (6) is shown in [5] via coding theorems
(called sequential code), hence this aspect will not be discussed. Rather, the
main emphasis of the paper is the mathematical formulation, the prove of
existence of solution to (6), the derivation of (7), the derivation of a
closed form expression for the causal rate distortion function, and some of
its properties.
The Shannon source code consists of an encoder-decoder pair. The encoder
observes a source sequence
$X^{\infty}\stackrel{{\scriptstyle\triangle}}{{=}}\\{X_{0},X_{1},\ldots\\}$
and generates a compressed representation $\\{Z_{0},Z_{1},\ldots\\}$. The
decoder upon observing the representation sequence $\\{Z_{0},Z_{1},\ldots\\}$
generates a reproduction sequence $Y_{i}=f_{i}(X^{\infty})$ of $X_{i}$, for
every time step $i$. The dependence of the reproduction sequence on the future
source symbols, in addition to its past and present symbols makes such a
decoder non-causal. In Neuhoff and Gilbert [4], a source code is defined as
causal if the reproduction sequence is such that
$f_{i}(X^{\infty})=f_{i}(\tilde{X}^{\infty})$ whenever
$X^{i}={\tilde{X}}^{i},~{}\forall i=0,1,\ldots$. The definition of a causal
code necessitates that any information theoretic causal rate distortion
function should lead to an optimal reconstruction conditional distribution
which is causally dependent on the source symbols, and (7) has this property.
The classical rate distortion function is defined via the mutual information
between $X^{n}$ and $Y^{n}$, namely, $I(X^{n};Y^{n})$ with average distortion
(2), and the code is assumed non-causal, leading to the well known optimal
reconstruction [1, 3]

$\displaystyle
P_{Y^{n}|X^{n}}^{*}(dy^{n}|x^{n})=\frac{e^{s\sum_{i=0}^{n}\rho_{0,i}(x^{i},y^{i})}P_{Y^{n}}^{*}(dy^{n})}{\int_{{\cal
Y}_{0,n}}e^{s\sum_{i=0}^{n}\rho_{0,i}(x^{i},y^{i})}P_{Y^{n}}^{*}(dy^{n})}$ (8)

Since by chain rule
$P_{Y^{n}|X^{n}}(dy^{n}|X^{n}=x^{n})=\otimes_{i=0}^{n}P_{Y_{i}|Y^{i-1},X^{n}=x^{n}}(dy_{i}|y^{i-1}=y^{i-1},X^{n}=x^{n})$,
the classical rate distortion theory gives a reconstruction $Y_{i}=y_{i}$
which depends on future values of the source symbols,
$(X_{i+1}=x_{i+1},\ldots,X_{n}=x_{n})$ in addition to its past reconstructions
$Y^{i-1}=y^{i-1}$, and past and present source symbols $X^{i}=x^{i}$. The
point to be made here is that, in general, aside from some special examples,
such as the i.i.d source and single letter distortion
$d_{0,n}=\sum^{n}_{i=0}\rho_{i}(x_{i},y_{i})$ [2] the reconstruction
conditional distribution and hence the decoder of the classical rate
distortion function is non-causal. On the other hand, a code is causal if the
reconstruction distribution is causal.

##  II PROBLEM FORMULATION

In this section, we introduce the set up of the problem on discrete time sets
$\mathbb{N}^{n}\stackrel{{\scriptstyle\triangle}}{{=}}\\{0,1,\ldots,n\\}$,
$n\in\mathbb{N}\stackrel{{\scriptstyle\triangle}}{{=}}\\{0,1,2,\ldots\\}$.
Assume all processes are defined on a complete probability space
$(\Omega,{\cal F}(\Omega),\mathbb{P})$ with filtration $\\{{\cal
F}_{t}\\}_{t\geq 0}$. The source and reconstruction alphabets are sequences of
Polish spaces [11] $\\{{\cal X}_{t}:t\in\mathbb{N}\\}$ and $\\{{\cal
Y}_{t}:t\in\mathbb{N}\\}$, respectively, (e.g., ${\cal Y}_{t},{\cal X}_{t}$
are complete separable metric spaces), associated with their corresponding
measurable spaces $({\cal X}_{t},{\cal B}({\cal X}_{t}))$ and $({\cal
Y}_{t},{\cal B}({\cal Y}_{t}))$ (e.g., ${\cal B}({\cal X}_{t})$ is a Borel
$\sigma-$algebra of subsets of the set ${\cal X}_{t}$ generated by closed
sets), $t\in\mathbb{N}$. Sequences of alphabets are identified with the
product spaces $({\cal X}_{0,n},{\cal B}({\cal
X}_{0,n}))\stackrel{{\scriptstyle\triangle}}{{=}}\times_{k=0}^{n}({\cal
X}_{k},{\cal B}({\cal X}_{k}))$, and $({\cal Y}_{0,n},{\cal B}({\cal
Y}_{0,n}))\stackrel{{\scriptstyle\triangle}}{{=}}\times_{k=0}^{n}({\cal
Y}_{k},{\cal B}({\cal Y}_{k}))$. The source and reconstruction are processes
denoted by
$X^{n}\stackrel{{\scriptstyle\triangle}}{{=}}\\{X_{t}:t\in\mathbb{N}^{n}\\}$,
$X:\mathbb{N}^{n}\times\Omega\mapsto{\cal X}_{t}$, and by
$Y^{n}\stackrel{{\scriptstyle\triangle}}{{=}}\\{Y_{t}:t\in\mathbb{N}^{n}\\}$,
$Y:\mathbb{N}^{n}\times\Omega\mapsto{\cal Y}_{t}$, respectively. Probability
measures on any measurable space $({\cal Z},{\cal B}({\cal Z}))$ are denoted
by ${\cal M}_{1}({\cal Z})$. It is assumed that the $\sigma$-algebras
$\sigma\\{X^{-1}\\}=\sigma\\{Y^{-1}\\}=\\{\emptyset,\Omega\\}$.

###### Definition II.1

Let $({\cal X},{\cal B}({\cal X})),({\cal Y},{\cal B}({\cal Y}))$ be
measurable spaces in which $\cal Y$ is a Polish Space.
A stochastic Kernel on $\cal Y$ given $\cal X$ is a mapping $q:{\cal B}({\cal
Y})\times{\cal X}\rightarrow[0,1]$ satisfying the following two properties:

1) For every $x\in{\cal X}$, the set function $q(\cdot;x)$ is a probability
measure (possibly finitely additive) on ${\cal B}({\cal Y}).$

2) For every $F\in{\cal B}({\cal Y})$, the function $q(F;\cdot)$ is ${\cal
B}({\cal X})$-measurable.
The set of all such stochastic Kernels is denoted by ${\cal Q}({\cal Y};{\cal
X})$.

An important notion is conditional independence. The Random Variable (R.V.)
${Z}$ is called conditional independent of R.V. $X$ given the R.V. $Y$ if and
only if $X\leftrightarrow Y\leftrightarrow Z$ forms a Markov chain in both
directions.

Stochastic kernels can be used to define non-causal and causal product
reconstruction kernels and associated rate distortion functions.

###### Definition II.2

Given measurable spaces $({\cal X}_{0,n},{\cal B}({\cal X}_{0,n}))$, $({\cal
Y}_{0,n},{\cal B}({\cal Y}_{0,n}))$, and their product spaces, data
compression channels are defined as follows.

  1. 1.

A Non-Causal Data Compression Channel is a stochastic kernel
$q_{0,n}(dy^{n};x^{n})\in{\cal Q}({\cal Y}_{0,n};{\cal
X}_{0,n}),n\in\mathbb{N}$.

  2. 2.

A Causal Product Data Compression Channel is a product of a sequence of causal
stochastic kernels defined by

$\displaystyle{\overrightarrow{q}}_{0,n}(dy^{n};x^{n})$
$\displaystyle=\otimes_{i=0}^{n}q_{i}(dy_{i};y^{i-1},x^{i})$

where $q_{i}\in{\cal Q}({\cal Y}_{i};{\cal Y}_{0,i-1}\times{\cal
X}_{0,i}),i=0,\ldots,n,~{}n\in\mathbb{N}$.

Note that classical rate distortion theory is concerned with finding the
optimal $P_{Y^{n}|X^{n}}(dy^{n}|X^{n}=x^{n})$, which is generally non-causal,
while in this paper the interest is to find the optimal causal product kernel.

###  II-A Causal and Classical Rate Distortion Functions

In this section the classical rate distortion function which has a non-causal
structure is reviewed, and then the causal rate distortion function is
defined.
Given a source probability measure ${\cal\mu}_{0,n}\in{\cal M}_{1}({\cal
X}_{0,n})$ (possibly finite additive) and a reconstruction Kernel
$q_{0,n}\in{\cal Q}({\cal Y}_{0,n};{\cal X}_{0,n})$, one can define three
probability measures as follows.

(P1): The joint measure $P_{0,n}\in{\cal M}_{1}({\cal Y}_{0,n}\times{\cal
X}_{0,n})$:

$\displaystyle P_{0,n}(G_{0,n})$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}(\mu_{0,n}\otimes
q_{0,n})(G_{0,n}),\>G_{0,n}\in{\cal B}({\cal X}_{0,n})\times{\cal B}({\cal
Y}_{0,n})$ $\displaystyle=\int_{{\cal
X}_{0,n}}q_{0,n}(G_{0,n,x^{n}};x^{n})\mu_{0,n}(d{x^{n}})$

where $G_{0,n,x^{n}}$ is the $x^{n}-$section of $G_{0,n}$ at point ${x^{n}}$
defined by
$G_{0,n,x^{n}}\stackrel{{\scriptstyle\triangle}}{{=}}\\{y^{n}\in{\cal
Y}_{0,n}:(x^{n},y^{n})\in G_{0,n}\\}$ and $\otimes$ denotes the convolution.

(P2): The marginal measure $\nu_{0,n}\in{\cal M}_{1}({\cal Y}_{0,n})$:

$\displaystyle\nu_{0,n}(F_{0,n})$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}P_{0,n}({\cal
X}_{0,n}\times F_{0,n}),~{}F_{0,n}\in{\cal B}({\cal Y}_{0,n})$
$\displaystyle=\int_{{\cal X}_{0,n}}q_{0,n}(({\cal X}_{0,n}\times
F_{0,n})_{{x}^{n}};{x}^{n})\mu_{0,n}(d{x^{n}})$ $\displaystyle=\int_{{\cal
X}_{0,n}}q_{0,n}(F_{0,n};x^{n})\mu_{0,n}(dx^{n})$

(P3): The product measure $\pi_{0,n}:{\cal B}({\cal X}_{0,n})\times{\cal
B}({\cal Y}_{0,n})\mapsto[0,1]$ of $\mu_{0,n}\in{\cal M}_{1}({\cal X}_{0,n})$
and $\nu_{0,n}\in{\cal M}_{1}({\cal Y}_{0,n})$:

$\displaystyle\pi_{0,n}(G_{0,n})$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}(\mu_{0,n}\times\nu_{0,n})(G_{0,n}),~{}G_{0,n}\in{\cal
B}({\cal X}_{0,n})\times{\cal B}({\cal Y}_{0,n})$ $\displaystyle=\int_{{\cal
X}_{0,n}}\nu_{0,n}(G_{0,n,x^{n}})\mu_{0,n}(dx^{n})$

The precise definition of mutual information between two sequences of Random
Variables $X^{n}$ and $Y^{n}$, denoted $I(X^{n};Y^{n})$ is defined via the
Kullback-Leibler distance (or relative entropy) between the joint probability
distribution of $(X^{n},Y^{n})$ and the product of its marginal probability
distributions of $X^{n}$ and $Y^{n}$, using the Radon-Nikodym derivative.
Hence, by the construction of probability measures (P1)-(P3), and the chain
rule of relative entropy [11]:

$\displaystyle
I(X^{n};Y^{n})\stackrel{{\scriptstyle\triangle}}{{=}}\mathbb{D}(P_{0,n}||\pi_{0,n})$
(9) $\displaystyle=\int_{{\cal X}_{0,n}\times{\cal
Y}_{0,n}}\log\Big{(}\frac{d(\mu_{0,n}\otimes
q_{0,n})}{d(\mu_{0,n}\times\nu_{0,n})}\Big{)}d(\mu_{0,n}\otimes q_{0,n})$
$\displaystyle=\int_{{\cal X}_{0,n}\times{\cal
Y}_{0,n}}\log\Big{(}\frac{q_{0,n}(dy^{n};x^{n})}{\nu_{0,n}(dy^{n})}\Big{)}$
$\displaystyle q_{0,n}(dy^{n};dx^{n})\mu_{0,n}(dx^{n})$
$\displaystyle=\int_{{\cal
X}_{0,n}}\mathbb{D}(q_{0,n}(\cdot;x^{n})||\nu_{0,n}(\cdot))\mu_{0,n}(dx^{n})$
$\displaystyle\equiv\mathbb{I}(\mu_{0,n};q_{0,n})$ (10)

Note that $(\ref{re3})$ states that mutual information is expressed as a
functional of $\\{\mu_{0,n},q_{0,n}\\}$ and it is denoted by
$\mathbb{I}(\mu_{0,n};q_{0,n})$. Note that necessary and sufficient conditions
for existence of a Radon-Nikodym derivative for finitely additive measures can
be found in [13]. Moreover, $I(X^{n};Y^{n})$ is also expressed by the sum of
two directed information as follows

$\displaystyle I(X^{n};Y^{n})$
$\displaystyle=I(X^{n}{\rightarrow}Y^{n})+I(X^{n}{\leftarrow}Y^{n})$ (11)

where

$\displaystyle I(X^{n}{\rightarrow}Y^{n})$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\sum_{i=0}^{n}I(X^{i};Y_{i}|Y^{i-1})$
(12) $\displaystyle I(X^{n}{\leftarrow}Y^{n})$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\sum_{i=0}^{n}I(Y^{i-1};X_{i}|X^{i-1})$
(13)

###### Definition II.3

(Classical Rate Distortion Function) Let $d_{0,n}:{\cal X}_{0,n}\times{\cal
Y}_{0,n}\rightarrow[0,\infty)$, be an ${\cal B}({\cal X}_{0,n})\times{\cal
B}({\cal Y}_{0,n})$-measurable distortion function, and let
$Q_{0,n}(D)\subset{\cal Q}({\cal Y}_{0,n};{\cal X}_{0,n})$ (assuming is non-
empty) denotes the average distortion or fidelity constraint defined by

$\displaystyle
Q_{0,n}(D)\stackrel{{\scriptstyle\triangle}}{{=}}\Big{\\{}q_{0,n}\in{\cal
Q}({\cal Y}_{0,n};{\cal X}_{0,n}):$ $\displaystyle\frac{1}{n+1}\int_{{\cal
X}_{0,n}}\int_{{\cal Y}_{0,n}}d_{0,n}(x^{n},y^{n})q_{0,n}(dy^{n};x^{n})$
$\displaystyle\mu_{0,n}(dx^{n})\leq D\Big{\\}},~{}D\geq 0$ (14)

The classical rate distortion function associated with the non-causal kernel
$q_{0,n}\in{\cal Q}({\cal Y}_{0,n};{\cal X}_{0,n})$ is defined by

$\displaystyle
R_{0,n}(D)\stackrel{{\scriptstyle\triangle}}{{=}}\inf_{q_{0,n}\in
Q_{0,n}(D)}\frac{1}{n+1}\mathbb{I}(\mu_{0,n};q_{0,n})$ (15)

while its operational meaning can be established via
${\lim}{\sup_{n\rightarrow\infty}}{R_{0,n}}$.

Existence in (15) is shown assuming $d_{0,n}(x^{n};\cdot)$ is bounded
continuous on ${\cal Y}_{0,n}$ and ${\cal Y}_{0,n}$ is compact, using weak-
convergence of probability measures in [3], and for more general
$d_{0,n}(x^{n};\cdot)$ which is only continuous in ${\cal Y}_{0,n}$ using
weak*-convergence of measures [14] on Polish spaces.
A version of the optimal reconstruction kernel which attains the infimum in
(15), [3] is

$\displaystyle q_{0,n}^{*}(dy^{n};x^{n})=\frac{\
e^{sd_{0,n}(x^{n},y^{n})}\nu_{0,n}^{*}(dy^{n})}{\int_{{\cal
Y}_{0,n}}e^{sd_{0,n}(x^{n},y^{n})}\nu_{0,n}^{*}(dy^{n})},\quad s\leq 0$ (16)

where $\nu_{0,n}^{*}\in{\cal M}_{1}({\cal Y}_{0,n})$ is the marginal of
$P_{0,n}^{*}=\mu_{0,n}\otimes q_{0,n}^{*}\in{\cal M}_{1}({\cal
X}_{0,n}\times{\cal Y}_{0,n})$ and $s\leq 0$ is the Lagrange multiplier
associated with the fidelity constraint $(\ref{dc1})$. Unfortunately, for
general sources and distortion function $d_{0,n}$, the optimal reconstruction
$q^{*}_{0,n}(dy^{n};x^{n})=\otimes^{n}_{i=0}q^{*}_{i}(dy_{i};y^{i-1},x^{n})$
is non-causal and introduces delay in the reconstruction processes. On the
other hand, if the solution (16) gives a reconstruction such that
$q^{*}_{0,n}(dy^{n};x^{n})={\overrightarrow{q}}^{*}_{0,n}(dy^{n};x^{n})=\otimes^{n}_{i=0}q^{*}_{i}(dy_{i};y^{i-1},x^{i})$
it will be causal. However, there are only limited examples in which
$(\ref{f6a})$ is causal on the source sequence. For single letter distortion
function
$d_{0,n}(x^{n},y^{n})=\frac{1}{n+1}\sum^{n}_{i=0}\rho_{i}(x_{i},y_{i})$ and
independent sources $\mu_{0,n}(dx^{n})=\otimes^{n}_{i=0}\mu_{i}(dx^{i})$
(e.g., $\\{X_{i}:i\in\mathbb{N}\\}$ are independent) the optimal
reconstruction $q^{*}_{0,n}(dy^{n};x^{n})$ factors into a product of causal
kernels $q^{*}_{0,n}(dy^{n};x^{n})=\otimes^{n}_{i=0}q_{i}(dy_{i},x_{i})$ [2].
This raises the question whether the classical rate distortion function can be
reformulated using the causal product
${\overrightarrow{q}}_{0,n}(dy^{n};x^{n})$.

The next lemma relates causal product reconstruction kernels, mutual
information, directed information, and conditional independence.

###### Lemma II.4

The following are equivalent for each $n\in\mathbb{N}$.

  1. 1.

$q_{0,n}(dy^{n};x^{n})={\overrightarrow{q}}_{0,n}(dy^{n};x^{n})$, as defined
in Definition II.2-2)

  2. 2.

For each $i=0,1,\ldots,n-1$,
$Y_{i}\leftrightarrow(X^{i},Y^{i-1})\leftrightarrow(X_{i+1},X_{i+2},\ldots,X_{n})$,
forms a Markov chain

  3. 3.

$I(X^{n};Y^{n})=I(X^{n}\rightarrow Y^{n})$

  4. 4.

$I(X^{n}\leftarrow Y^{n})=0$

  5. 5.

For each $i=0,1,\ldots,n-1$, $Y^{i}\leftrightarrow X^{i}\leftrightarrow
X_{i+1}$ forms a Markov chain

Proof. Omitted due to space limitation.
According to Lemma II.4 any source with a satisfying conditional distribution
$P_{X_{i}|X^{i-1},Y^{i-1}}(dx_{i}|X^{i-1}=x^{i-1},Y^{i-1}=y^{i-1})=P_{X_{i}|X^{i-1}}(dx_{i}|X^{i-1}=x^{i-1}),~{}P-a.s.,$
$\forall{i}\in\mathbb{N}$ is equivalent to any of the equivalent statements of
Lemma II.4. Therefore, for such a source the mutual information becomes

$\displaystyle I(X^{n};Y^{n})=I(X^{n}{\rightarrow}Y^{n})$
$\displaystyle=\int_{{\cal X}_{0,n}\times{\cal
Y}_{0,n}}\log\Big{(}\frac{\overrightarrow{q}_{0,n}(dy^{n};x^{n})}{\nu_{0,n}(dy^{n})}\Big{)}$
$\displaystyle\overrightarrow{q}_{0,n}(dy^{n};dx^{n})\mu_{0,n}(dx^{n})$ (17)
$\displaystyle\equiv{\mathbb{I}}(\mu_{0,n};\overrightarrow{q}_{0,n})$ (18)

where (18) states that $I(X^{n};Y^{n})$ is a functional of
$\\{\mu_{0,n},{\overrightarrow{q}}_{0,n}\\}$. Hence, causal rate distortion is
defined by optimizing ${\mathbb{I}}(\mu_{0,n};\overrightarrow{q}_{0,n})$ over
${\overrightarrow{q}}_{0,n}$ which satisfies a distortion constraint.

###### Definition II.5

(Causal Rate Distortion Function) Suppose
$d_{0,n}\stackrel{{\scriptstyle\triangle}}{{=}}\sum^{n}_{i=0}\rho_{0,i}(x^{i},y^{i})$,
where $\rho_{0,i}:{\cal X}_{0,i}\times{\cal Y}_{0,i}\rightarrow[0,\infty)$, is
a sequence of ${\cal B}({\cal X}_{0,i})\times{\cal B}({\cal
Y}_{0,i})$-measurable distortion functions, and let
$\overrightarrow{Q}_{0,n}(D)$ (assuming is non-empty) denotes the average
distortion or fidelity constraint defined by

$\displaystyle\overrightarrow{Q}_{0,n}(D)\stackrel{{\scriptstyle\triangle}}{{=}}\Big{\\{}\overrightarrow{q}_{0,i}\in{\cal
M}_{1}({\cal Y}_{0,i}),0\leq{i}\leq{n}:$
$\displaystyle\frac{1}{n+1}\sum_{i=0}^{n}\int_{{\cal X}_{0,i}}\int_{{{\cal
Y}}_{0,i}}\rho_{0,i}({x^{i}},{y^{i}})\overrightarrow{q}_{0,i}(d{y}^{i};{x}^{i})$
$\displaystyle\mu_{0,i}(d{x}^{i})\leq D\Big{\\}},~{}D\geq 0$ (19)

The causal rate distortion function associated with the causal product kernel
${\overrightarrow{q}}_{0,n}\in{\overrightarrow{Q}}_{0,n}(D)$ is defined by

$\displaystyle{\overrightarrow{R}}_{0,n}(D)\stackrel{{\scriptstyle\triangle}}{{=}}\inf_{{\overrightarrow{q}_{0,n}\in\overrightarrow{Q}_{0,n}(D)}}\frac{1}{n+1}{\mathbb{I}}(\mu_{0,n};\overrightarrow{q}_{0,n})$
(20)

while its operational meaning can be established via
$\lim\sup_{n\rightarrow{\infty}}{\overrightarrow{R}}_{0,n}$.

Clearly, ${\overrightarrow{R}}_{0,n}(D)$ is characterized by minimizing
directed information or equivalently
$\mathbb{I}(\mu_{0,n};\overrightarrow{q}_{0,n})$ over the causal product
measure ${\overrightarrow{q}}_{0,n}\in{\overrightarrow{Q}}_{0,n}(D)$.

###### Lemma II.6

$\overrightarrow{q}_{0,n}\in{\cal M}_{1}({\cal Y}_{0,n})$ is uniquely
determined by $\\{q_{i}\in{\cal Q}_{i}({\cal Y}_{i};{\cal
Y}_{0,i-1}\times{\cal X}_{0,i})\\}_{i=0}^{n}$ and vice-versa, $P-a.s$.

Proof. For densities this result is derived in [15].

##  III EXISTENCE OF OPTIMAL CAUSAL PRODUCT RECONSTRUCTION KERNEL

In this section, appropriate topologies and function spaces are employed to
show existence of the minimizing causal product kernel in $(\ref{ex12})$. In
the process we also show existence for $R_{0,n}(D)$.

###  III-A Abstract Spaces

Let $BC({\cal Y}_{0,n})$ denote the vector space of bounded continuous real
valued functions defined on the Polish space ${\cal Y}_{0,n}$. Furnished with
the sup norm topology, this is a Banach space. The topological dual of
$BC({\cal Y}_{0,n})$ denoted by $\Big{(}BC({\cal Y}_{0,n})\Big{)}^{*}$ is
isometrically isomorphic to the Banach space of finitely additive regular
bounded signed measures on ${\cal Y}_{0,n}$ [7], denoted by $M_{rba}({\cal
Y}_{0,n})$. Let $\Pi_{rba}({\cal Y}_{0,n})\subset M_{rba}({\cal Y}_{0,n})$
denote the set of regular bounded finitely additive probability measures on
${\cal Y}_{0,n}$. Clearly if ${\cal Y}_{0,n}$ is compact, then
$\Big{(}BC({\cal Y}_{0,n})\Big{)}^{*}$ will be isometrically isomorphic to the
space of countably additive signed measures, as in [3]. Denote by
$L_{1}(\mu_{0,n},BC({\cal Y}_{0,n}))$ the space of all $\mu_{0,n}$-integrable
functions defined on ${\cal X}_{0,n}$ with values in $BC({\cal Y}_{0,n}),$ so
that for each $\phi\in L_{1}(\mu_{0,n},BC({\cal Y}_{0,n}))$ its norm is
defined by

$\displaystyle\parallel\phi\parallel_{\mu_{0,n}}\stackrel{{\scriptstyle\triangle}}{{=}}\int_{{\cal
X}_{0,n}}||\phi(x^{n})(\cdot)||_{BC({\cal Y}_{0,n})}\mu_{0,n}(dx^{n})<\infty$

The norm topology $\parallel{\phi}\parallel_{\mu_{0,n}}$, makes
$L_{1}(\mu_{0,n},BC({\cal Y}_{0,n}))$ a Banach space, and it follows from the
theory of “lifting” [10] that the dual of this space is
$L_{\infty}^{w}(\mu_{0,n},M_{rba}({\cal Y}_{0,n}))$, denoting the space of all
$M_{rba}({\cal Y}_{0,n})$ valued functions $\\{q\\}$ which are
weak∗-measurable in the sense that for each $\phi\in BC({\cal Y}_{0,n}),$
$x^{n}\longrightarrow
q_{x^{n}}(\phi)\stackrel{{\scriptstyle\triangle}}{{=}}\int_{{\cal
Y}_{0,n}}\phi(y^{n})q(dy^{n};x^{n})$ is $\mu_{0,n}$-measurable and
$\mu_{0,n}$-essentially bounded.

###  III-B Weak∗-Compactness and Existence

Define an admissible set of stochastic kernels associated with classical rate
distortion function by

$\displaystyle
Q_{ad}\stackrel{{\scriptstyle\triangle}}{{=}}L_{\infty}^{w}(\mu_{0,n},\Pi_{rba}({\cal
Y}_{0,n}))\subset L_{\infty}^{w}(\mu_{0,n},M_{rba}({\cal Y}_{0,n}))$

Clearly, $Q_{ad}$ is a unit sphere in $L_{\infty}^{w}(\mu_{0,n},M_{rba}({\cal
Y}_{0,n}))$. For each $\phi{\in}L_{1}(\mu_{0,n},BC({\cal Y}_{0,n}))$ we can
define a linear functional on $L_{\infty}^{w}(\mu_{0,n},M_{rba}({\cal
Y}_{0,n}))$ by

$\displaystyle\ell_{\phi}(q_{0,n})\stackrel{{\scriptstyle\triangle}}{{=}}\frac{1}{n+1}\int_{{\cal
X}_{0,n}}\Big{(}\int_{{\cal Y}_{0,n}}\phi(x^{n},y^{n})$ $\displaystyle
q_{0,n}(dy^{n};x^{n})\Big{)}\mu_{0,n}(dx^{n})$

This is a bounded, linear and weak∗-continuous functional on
$L_{\infty}^{w}(\mu_{0,n},M_{rba}({\cal Y}_{0,n}))$. For $d_{0,n}:{\cal
X}_{0,n}\times{\cal Y}_{0,n}\rightarrow[0,\infty)$ measurable and
$d_{0,n}{\in}L_{1}(\mu_{0,n},BC({\cal Y}_{0,n}))$ the distortion constraint
set of the classical rate distortion function is

$\displaystyle
Q_{0,n}(D)\stackrel{{\scriptstyle\triangle}}{{=}}\\{q{\in}Q_{ad}:\frac{1}{n+1}\ell_{d_{0,n}}(q_{0,n}){\leq}D\\}$

It can be shown that $Q_{0,n}(D)$ is bounded and weak∗-closed subset of
$Q_{ad}$ and hence weak∗-compact (Compactness of $Q_{ad}$ follows from
Alaoglu’s Theorem [7],[12]).
Next, we define the set of causal product kernels as follows.

$\displaystyle{\overrightarrow{\Pi}}_{rba}({\cal Y}_{0,n})$
$\displaystyle=\Big{\\{}{\overrightarrow{q}}_{0,n}(dy^{n};x^{n})\stackrel{{\scriptstyle\triangle}}{{=}}\otimes_{i=1}^{n}{q}_{i}(dy_{i};y^{i-1},x^{i}):$
$\displaystyle{q}_{i}(dy_{i};y^{i-1},x^{i})\in{\Pi}_{rba}({\cal
Y}_{i}),\>i\in\mathbb{N}^{n}\Big{\\}}$

where $L_{\infty}^{w}(\mu_{0,n},{\overrightarrow{\Pi}}_{rba}({\cal Y}_{0,n}))$
denotes the space of all ${\overrightarrow{\Pi}}_{rba}({\cal Y}_{0,n})$ valued
functions $\\{\overrightarrow{q}\\}$ which are weak∗-measurable in the sense
that for each $\phi\in BC({\cal Y}_{0,n}),$
$x^{n}\rightarrow{\overrightarrow{q}}_{x^{n}}(\phi)\stackrel{{\scriptstyle\triangle}}{{=}}\int_{{\cal
Y}_{0,n}}\phi(y^{n}){\overrightarrow{q}}(dy^{n};x^{n})$ is
$\mu_{0,n}$-measurable and $\mu_{0,n}$-essentially bounded.
Define the admissible set of causal product stochastic kernels associated with
the causal rate distortion function by

$\displaystyle{\overrightarrow{Q}}_{ad}$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}L_{\infty}^{w}(\mu_{0,n},{\overrightarrow{\Pi}}_{rba}({\cal
Y}_{0,n}))$

Clearly, ${\overrightarrow{Q}}_{ad}=\\{q_{0,n}\in
Q_{ad}:q_{0,n}(dy^{n};x^{n})=\overrightarrow{q}_{0,n}(dy^{n};x^{n})\\}$. For
$d_{0,n}:{\cal X}_{0,n}\times{\cal Y}_{0,n}\rightarrow[0,\infty)$ which is
measurable and $d_{0,n}{\in}L_{1}(\mu_{0,n},BC({\cal Y}_{0,n}))$ the
distortion constraint of causal rate distortion function is

$\displaystyle{\overrightarrow{Q}_{{0,n}}(D)}\stackrel{{\scriptstyle\triangle}}{{=}}\Big{\\{}{\overrightarrow{q}}_{0,n}\in{\overrightarrow{Q}}_{ad}:$
$\displaystyle\frac{1}{n+1}\ell_{d_{0,n}}({\overrightarrow{q}}_{0,n})\stackrel{{\scriptstyle\triangle}}{{=}}\int_{{\cal
X}_{0,n}}\biggr{(}\int_{{\cal Y}_{0,n}}d_{0,n}(x^{n},y^{n})$
$\displaystyle{\overrightarrow{q}}_{0,n}(dy^{n};x^{n})\biggr{)}\mu_{0,n}(dx^{n})\leq
D\Big{\\}}$

###### Assumptions III.1

We make the following assumptions.

  1. 1.

The set $\overrightarrow{Q}_{ad}$ is weak∗-closed.

  2. 2.

The set ${\overrightarrow{Q}_{0,n}(D)}$ is non-empty.

###### Lemma III.2

Suppose Assumptions III.1 hold. Let ${\cal X}_{0,n},{\cal Y}_{0,n}$ be two
Polish spaces and $d_{0,n}:{\cal X}_{0,n}\times{\cal
Y}_{0,n}\rightarrow[0,\infty],$ a measurable, non-negative, extended real
valued function, such that $d_{0,n}{\in}L_{1}(\mu_{0,n},BC({\cal Y}_{0,n}))$.
For any $D\in[0,\infty)$, the set ${\overrightarrow{Q}_{0,n}(D)}$ is
weak∗-compact.

Proof. By Assumptions III.1, $\overrightarrow{Q}_{ad}$ is a weak∗-closed,
hence as a subset of a weak∗-compact set $Q_{ad}$ it is weak∗-compact. Also,
under assumptions III.1, ${\overrightarrow{Q}_{0,n}(D)}$ is bounded and
weak∗-closed and hence it is weak∗-compact (as a weak∗-closed subset of the
weak∗-compact set ${\overrightarrow{Q}}_{ad}$) $\bullet$

###### Theorem III.3

Under Assumptions III.1, ${\overrightarrow{R}}_{0,n}(D)$ has a minimum.

Proof. Follows from Lemma III.2 and the lower semi-continuity of
${\mathbb{I}}(\mu_{0,n};\cdot)$ on ${\overrightarrow{Q}}_{ad}$ $\bullet$

##  IV NECESSARY CONDITIONS OF OPTIMALITY OF CAUSAL PRODUCT RATE DISTORTION
FUNCTION

In this section the form of the optimal causal product reconstruction kernels
is derived. The method is based on calculus of variations on the space of
measures [9].

###### Theorem IV.1

Suppose
${\mathbb{I}}_{\mu_{0,n}}({\overrightarrow{q}}_{0,n})\stackrel{{\scriptstyle\triangle}}{{=}}{\mathbb{I}}(\mu_{0,n};\overrightarrow{q}_{0,n})$
is well defined for every ${\overrightarrow{q}}_{0,n}\in
L_{\infty}^{w}({\mu_{0,n},\overrightarrow{\Pi}}_{rba}({\cal Y}_{0,n}))$
possibly taking values from the set $[0,\infty].$ Then
${\overrightarrow{q}}_{0,n}\rightarrow{\mathbb{I}}_{\mu_{0,n}}({\overrightarrow{q}}_{0,n})$
is Gateaux differentiable at every point in
$L_{\infty}^{w}({\mu_{0,n},\overrightarrow{\Pi}}_{rba}({\cal Y}_{0,n})),$ and
the Gateaux derivative at the point ${\overrightarrow{q}}_{0,n}^{0}$ in the
direction ${\overrightarrow{q}}_{0,n}-{\overrightarrow{q}}_{0,n}^{0}$ is given
by

$\displaystyle\delta{\mathbb{I}}_{\mu_{0,n}}({\overrightarrow{q}}_{0,n}^{0};{\overrightarrow{q}}_{0,n}-{\overrightarrow{q}}_{0,n}^{0})$
$\displaystyle=\int_{{\cal X}_{0,n}}\int_{{\cal
Y}_{0,n}}\log\Bigg{(}\frac{{\overrightarrow{q}}_{0,n}^{0}(dy^{n};x^{n})}{\nu_{0,n}^{0}(dy^{n})}\Bigg{)}$
$\displaystyle({\overrightarrow{q}}_{0,n}-{\overrightarrow{q}}_{0,n}^{0})(dy^{n};x^{n})\mu_{0,n}(dx^{n})$

where $\nu_{0,n}^{0}\in{\cal M}_{1}({\cal Y}_{0,n})$ is the marginal measure
corresponding to
${\overrightarrow{q}}_{0,n}^{0}\otimes\mu_{0,n}(dx^{n})\in{\cal M}_{1}({\cal
Y}_{0,n}\times{\cal X}_{0,n})$.

Proof. The proof is based on the fact that the causal product stochastic
kernel ${\overrightarrow{q}}_{0,n}$ is used to show the existence of Gateaux
Differential [9] rather than for individual causal stochastic kernel
$q_{i}(dy_{i};y^{i-1},x^{i})$, $i\in\mathbb{N}^{n}$ $\bullet$

The constrained problem defined by (20) can be reformulated using Lagrange
multipliers as follows (equivalence of constrained and unconstrained problems
follows from [9]).

$\displaystyle{\overrightarrow{R}}_{0,n}(D)=\inf_{{\overrightarrow{q}}_{0,n}\in{\overrightarrow{Q}}_{ad}}\Big{\\{}\frac{1}{n+1}{{\mathbb{I}}}(\mu_{0,n};{\overrightarrow{q}}_{0,n})$
$\displaystyle-s(\ell_{{d}_{0,n}}({\overrightarrow{q}}_{0,n})-D)\Big{\\}}$
(21)

and $s\in(-\infty,0]$ is the Lagrange multiplier.

###### Theorem IV.2

Suppose $d_{0,n}(x^{n},y^{n})=\sum_{i=0}^{n}\rho_{0,i}(x^{i},y^{i})$ and the
assumptions of Lemma III.2 hold. The infimum in $(\ref{ex13})$ is attained at
$\overrightarrow{q}^{*}_{0,n}\in
L_{\infty}^{w}(\mu_{0,n},{\overrightarrow{\Pi}}_{rba}({\cal Y}_{0,n}))$ given
by

$\displaystyle\overrightarrow{q}^{*}_{0,n}(dy^{n};x^{n})=\otimes_{i=0}^{n}\frac{e^{s\rho_{i}(x^{i},y^{i})}\nu^{*}_{i}(dy^{i};y^{i-1})}{\int_{{\cal
Y}_{i}}e^{s\rho_{i}(x^{i},y^{i})}\nu^{*}_{i}(dy_{i};y^{i-1})}$ (22)

and $\nu^{*}_{i}(dy_{i};y^{i-1})\in{\cal Q}({\cal Y}_{i};{\cal Y}_{0,{i-1}})$.
The causal rate distortion function is given by

$\displaystyle{\overrightarrow{R}}_{0,n}(D)=sD-\frac{1}{n+1}\sum_{i=0}^{n}\int_{{{\cal
X}_{0,i}}\times{{\cal Y}_{0,i-1}}}$ $\displaystyle\log\Big{(}\int_{{\cal
Y}_{i}}e^{s\rho_{i}(x^{i},y^{i})}\nu^{*}_{i}(dy_{i};y^{i-1})\Big{)}$
$\displaystyle{{\overrightarrow{q}}^{*}_{0,i-1}}(dy^{i-1};x^{i-1})\otimes\mu_{0,i}(dx^{i})$
(23)

If ${\overrightarrow{R}}_{0,n}(D)>0$ then $s<0$ and

$\displaystyle\frac{1}{n+1}\sum_{i=0}^{n}\int_{{\cal X}_{0,i}}\int_{{\cal
Y}_{0,i}}\rho_{0,i}(x^{i},y^{i}){\overrightarrow{q}}^{*}_{0,i}(dy^{i};x^{i})\mu_{0,i}(dx^{i})=D$

Proof. The fully unconstraint problem of (21) is obtained by introducing
another Lagrange multiplier. Using this and Theorem IV.1 we obtain (22) and
(23) $\bullet$

##  V PROPERTIES OF CAUSAL RATE DISTORTION FUNCTION

In this section, we present some important properties of the causal rate
distortion function as it is defined in (20).

###### Theorem V.1

  1. 1.

${\overrightarrow{R}}_{0,n}(D)$ is a convex, non-increasing function of $D$

  2. 2.

If $\rho_{0,i}\in L^{1}(\pi_{0,i})$ then
a)
${\overrightarrow{R}}_{0,n}(\frac{1}{n+1}\sum_{i=0}^{n}E_{\pi_{0,i}}(\rho_{0,i}))=0$;
b) ${\overrightarrow{R}}_{0,n}(D)$ is non-increasing for $D\in[0,D_{max}]$
where $D_{max}=\frac{1}{n+1}\sum_{i=0}^{n}E_{\pi_{0,i}}(\rho_{0,i})$ and
${\overrightarrow{R}}_{0,n}(D)=0$ for any $D\geq D_{max}$

  3. 3.

${\overrightarrow{R}}_{0,n}(D)>0$ for all $D<D_{max}$ and
${\overrightarrow{R}}_{0,n}(D)=0$ for all $D\geq D_{max}$, where

$\displaystyle D_{max}=\min_{\\{y^{n}\\}\in{\cal
Y}_{0,n}}\frac{1}{n+1}\sum_{i=0}^{n}\int_{{\cal
X}_{0,i}}\rho_{0,i}(x^{i},y^{i})\mu_{0,i}(dx^{i})$

if such a minimum exists.

Proof. Omitted due to space limitation.

##  VI CONCLUSION AND FUTURE WORK

###  VI-A Conclusion

The solution of the causal rate distortion function subject to a reproduction
kernel which is a product of causal kernels is presented, on abstract
alphabets. Some of its properties are also presented. It is believed that the
optimal reconstruction kernel as a product of causal kernels has several
implications in applications where causality of the decoder as a function of
the source is of concern.

###  VI-B Future Work

Examples are currently under investigation, and will be presented at the final
version of the paper.

##  VII APPENDIX

## References

  * [1] T. Berger, Rate Distortion Theory: A Mathematical Basis for Data Compression.  Prentice Hall, Englewood Cliffs, NJ, 1971.
  * [2] T. Cover and J. Thomas, Elements of Information Theory.  John Wiley & Sons, 1991.
  * [3] I. Csiszár, “On an extremum problem of information theory”, Studia Scientiarum Mathematicarum Hungarica, vol. 9, pp. 57–71, 1974.
  * [4] D. L. Neuhoff and R. Kent Gilbert, “Causal Source Codes”, IEEE Transactions on Information Theory, vol. IT-28, No.5, pp. 701–713, 1982.
  * [5] S. Tatikonda, “Control Over Communication Constraints”, PhD Dissertation, M.I.T., Cambridge, MA, 2000.
  * [6] J. Massey, “Causality, Feedback and Directed Information”, in the IEEE International Symposium on Information Theory and its Applications, pp. 303–305, Nov. 27–30, Hawaii, U.S.A.1990.
  * [7] N. Dunford and J. T. Schwartz, Linear Operators, Part I: General Theory.  Interscience Publishers, Inc., New York, 1958.
  * [8] R. M. Gray, Entropy and Information Theory.  Springer-Verlag, 1990.
  * [9] D. G. Luenberger, Optimization by Vector Space Methods.  John Wiley & Sons, 1969.
  * [10] A. Ionescu Tulcea & C. Ionescu Tulcea, Topics in the Theory of Lifting,  Springer Verlag, Berlin, Heidelberg, New York, 1969.
  * [11] P. Dupuis and R. S. Ellis, A Weak Convergence Approach to the theory of Large Deviations.  John Wiley & Sons, 1997.
  * [12] W. Rudin, Functional analysis.  McGraw-Hill, 1991.
  * [13] H.B.Maynard, A Radon-Nikodym Theorem for Finitely Additive Bounded Measures,  Pacific Journal of Mathematics, 83(2), 1979, pp. 401-413.
  * [14] F. Rezaei, N. U. Ahmed and C. D. Charalambous, Rate Distortion Theory for General Sources With Potential Application to Image Compression,  International Journal of Applied Mathematical Sciences, vol. 3 No. 2, 2006, pp. 141-165.
  * [15] H. H. Permuter, T. Weissman, A. Goldsmith, “Finite State Channels with Time-Invariant Deterministic Feedback”, IEEE Transactions on Information Theory, vol.IT-55, No. 2, pp. 644-662, February 2009.


11institutetext: Vladislav Taynitskiy 22institutetext: St. Petersburg State
University,
Faculty of Applied Mathematics and Control Processes,
Universitetskii prospekt 35, Petergof, Saint-Petersburg, Russia, 198504;
22email<EMAIL_ADDRESS>33institutetext: Elena Gubar 44institutetext: St.
Petersburg State University, St. Petersburg State University,
Faculty of Applied Mathematics and Control Processes,
Universitetskii prospekt 35, Petergof, Saint-Petersburg, Russia, 198504;
44email<EMAIL_ADDRESS>55institutetext: Quanyan Zhu 66institutetext:
Department of Electrical and Computer Engineering, Tandon School of
Engineering,
New York University, Brooklyn, NY, USA, 11201. 66email<EMAIL_ADDRESS>
# Optimal Impulse Control of SIR Epidemics over Scale-Free Networks

Vladislav Taynitskiy1 Elena Gubar1 Quanyan Zhu2

###### Abstract

Recent wide spreading of Ransomware has created new challenges for
cybersecurity over large-scale networks. The densely connected networks can
exacerbate the spreading and makes the containment and control of the malware
more challenging. In this work, we propose an impulse optimal control
framework for epidemics over networks. The hybrid nature of discrete-time
control policy of continuous-time epidemic dynamics together with the network
structure poses a challenging optimal control problem. We leverage the
Pontryagin’s minimum principle for impulsive systems to obtain an optimal
structure of the controller and use numerical experiments to corroborate our
results.

##  1 Introduction

Malware spreading becomes a more prevalent issue recently as the number of
devices and their connections grow exponentially. Many devices that are
connected to the Internet do not have strong protections, and they contain
cyber vulnerabilities that create a fast spreading of malware over large
networks. A higher level of connectivity of the network is often desired for
information spreading. However, in the context of malware, the high
connectivity can exacerbate the spreading and makes the containment and
control of the malware more challenging. One example is the recent Ransomware
Ransomware ; Luo  that spreads over the Internet with the objective to lock
the files of a victim using cryptographic techniques and demand a ransom
payment to decrypt them. The worldwide spread of WannaCry ransomware has
affected more than 200,000 computers across 150 countries and caused billions
of dollars of damages. Hence it is critical to take into account the network
structure when developing control policies to control the infection dynamics.

In this paper, we investigate a continuous-time Susceptible-Infected-Recovered
(SIR) epidemic model over large-scale networks. The malware control mechanism
is to patch an optimal fraction of the infected nodes at discrete points in
time. Such mechanism is also known as an optimal impulse controller. The
hybrid nature of discrete-time control policy of continuous-time epidemic
dynamics together with the network structure poses a challenging optimal
control problem. We leverage the Pontryagin’s minimum principle for impulsive
systems to obtain an optimal structure of the controller and use numerical
experiments to demonstrate the computation of the optimal control and the
controlled dynamics. This work extends the investigation of previous related
works Zhu ; Gubar ; Taynitskiy  to a new paradigm of coupled epidemic models
and the regime of optimal impulsive control.

The rest of the paper is organized as follows. Section II presents the
controlled SIR mathematical model. Section III shows the structure of optimal
control policies. Section IV presents numerical examples. Section V concludes
the paper.

##  2 The model

In this section, we formulate a model of spreading of malware in the network
of $N$ nodes use the modification of classical SIR model. As in previous works
Gubar , Taynitskiy:2016  two different forms of malware with different
strengths spread over the network simultaneously, we denote them as $M_{1}$
and $M_{2}$. We also assume that a structure of population is described by the
scale free network Vespignani , Porokhnyavaya . Normally, as SIR model points,
all nodes in the population are divided into three groups: Susceptible $(S)$,
Infected $(I)$ and Recovered $(R)$, Altman. Susceptible is a group of nodes
which are not infected by any malware, but may be invaded by any forms of
virus. The Infected nodes are those that have been attacked by the virus and
the Recovered is a group of recovered nodes. In modified model subgroup of
Infected nodes also is brunched into two subgroups $I_{1}$ and $I_{2}$, where
nodes in $I_{i}$ are infected by malware $M_{i},i=1,2$ respectively. We
formulate the epidemic process as a system of nonlinear differential
equations, where $n_{S}$, $n_{M_{1}}$, $n_{M_{2}}$ and $n_{R}$ correspond to
the number of susceptible, infected and recovered nodes, respectively. In
current model the connections between nodes are described by the scale-free
network, then we will use the following notation: $S_{k}(t)$ and $R_{k}(t)$
are fractions of Susceptible and Recovered nodes with degree $k$ at time
moment $t$, $I^{1}_{k}(t),I^{2}_{k}(t)$ are fractions of Infected nodes with
degree $k$. At each time moment $t\in[0,T]$ the number of nodes is constant
and equal $N$, and the following condition
$S_{k}(t)+I^{1}_{k}(t)+I^{2}_{k}(t)+R_{k}(t)=1$ is satisfied. The process of
spreading is defined by the system of ordinary differential equations:

$\begin{array}[]{l}\displaystyle\frac{dS_{k}}{dt}=-\delta_{1k}S_{k}I^{1}_{k}\Theta_{1}-\delta_{2k}S_{k}I^{2}_{k}\Theta_{2};\\\
\displaystyle\frac{dI^{1}_{k}}{dt}=\delta_{1k}S_{k}I^{1}_{k}\Theta_{1}-\sigma^{1}_{k}I^{1}_{k};\\\
\displaystyle\frac{dI^{2}_{k}}{dt}=\delta_{2k}S_{k}I^{2}_{k}\Theta_{2}-\sigma^{2}_{k}I^{2}_{k},\\\
\displaystyle\frac{dR_{k}}{dt}=\sigma^{1}_{k}I^{1}_{k}+\sigma^{2}_{k}I^{2}_{k},\end{array}$
(1)

where $\delta_{ik}(k)$ is the infections rate for the first type of malware
$i$ if a susceptible node has a contact with infected node with the degree
$k$, $\sigma^{i}_{k}$ is recovery rate.

We consider the graph generated by using the algorithm devised in barabasi .
We start from a small number $m_{0}$ of disconnected nodes; every time step a
new node is added, with $m$ links that are connected to an old node $i$ with
$k_{i}$ links according to the probability $k_{i}/\sum_{j}k_{j}$. After
iterating this scheme a sufficient number of times, we obtain a network
composed by $N$ nodes with connectivity distribution $P(k)\approx k^{-3}$ and
average connectivity $\langle k\rangle=2m$. In this work we take $m=4$.

At the initial time moment $t=0$, the most number of nodes belong to
Susceptible group and only a small fraction of Infected by malwares $M_{1}$ or
$M_{2}$. Initial state for system (1) is $0<S_{k}(0)<1$, $0<I^{1}_{k}(0)<1$,
$0<I^{2}_{k}(0)<1$, $R(0)=1-S_{k}(0)-I^{1}_{k}(0)-I^{2}_{k}(0)$. Analogously
with Fu , Vespignani  we define parameter $\Theta_{i}(t)$ as

$\Theta_{i}(\lambda_{i})=\sum\limits_{k^{\prime}}\frac{\delta_{ik}P(k^{\prime}|k)I^{i}_{k^{\prime}}}{k^{\prime}},\
i=1,2,$ (2)

where $\delta_{ik}$ denotes the infectivity of a node with degree $k$ and
$\lambda_{i}=\delta_{ik}/\sigma^{i}_{k}$ an effective spreading rate.
$P(k^{\prime}|k)$ describes the probability of a node with degree $k$ pointing
to a node with degree $k^{\prime}$, and
$P(k^{\prime}|k)=\frac{k^{\prime}P(k^{\prime})}{k^{\prime}}$, where $\langle
k\rangle=\sum\limits_{k^{\prime}}kP(k)$. For scale-free node distribution
$P(k)=C^{-1}k^{-2-\gamma},\ 0<\gamma\leq 1$, where $C=\zeta(2+\gamma)$ is
Riemann’s zeta function, which provides an appropriate normalization constant
for sufficiently large networks. In the SF model considered here, we have a
connectivity distribution $P(k)=2m^{2}/k^{-3}$, where $k$ is approximated as a
continuous variable. According to Vespignani  we can rewrite (2) as

$\Theta_{i}(\lambda_{i})=\frac{e^{-1/m\lambda_{i}}}{m\lambda_{i}},\ i=1,2.$
(3)

##  3 Impulse control problem

Previously it was shown in Vespignani  a small fraction of the infected nodes
might be survived on small segments of the network and can provoke new waves
of epidemics. This cycled process recalls the behavior the virus of influenza
which causes a seasonally periodic epidemic, Agur . Hence the control of the
epidemic process can be formulated as an impulse control problem in which a
series of impulses of antivirus patches are designed to reduce the
periodically incipient zones of infected nodes. We extend the model (1) to
present an impulse control problem for episodic attacks of the malware and
obtain the optimal strategy of application of antivirus software to damp the
spreading of malware at discrete time moments.

We suppose that impulses occur at time
$\tau^{i}_{k,1},\ldots,\tau^{i}_{k,q_{i}(k)}$, where $q_{i}(k)$ describes the
number of launching of impulse controls for nodes with $k$ degrees, index $i$
indicates the type of malware. We also assume that on the time intervals
$(\tau^{i}_{k,j},\tau^{i}_{k,j+1}]$ system (1) describes the behavior of
malware in the network. We have reformulated epidemic model to describe the
situation with two types of malware for all time periods except the sequence
of times $\tau_{k,j}^{i+}$, $j=1,\ldots,q_{i}(k)$, $i=1,2$. Additionally, we
set $S(\tau_{k,j}^{i})=S(\tau_{k,j}^{i-})$,
$I_{1}(\tau_{k,j}^{i})=I_{1}(\tau_{k,j}^{i-})$,
$I_{2}(\tau_{k,j}^{i})=I_{2}(\tau_{k,j}^{i-})$,
$R(\tau_{k,j}^{i})=R(\tau_{k,j}^{i-})$.

The system after activation of impulses at time moment $\tau_{k,j}^{i+}$ is:

$\begin{array}[]{l}S_{k}(\tau_{k,j}^{i+})=S_{k}(\tau_{k,j}^{i}),\\\
{I^{1}_{k}}({\tau_{k,j}^{i+}})=I^{1}_{k}(\tau_{k,j}^{i})-\nu^{1}_{k}(\tau_{k,j}^{i}),\\\
{I^{2}_{k}}({\tau_{k,j}^{i+}})=I^{2}_{k}(\tau_{k,j}^{i})-\nu^{2}_{k}(\tau_{k,j}^{i}),\\\
{R_{k}}({\tau_{k,j}^{i+}})={R_{k}}(\tau_{k,j}^{i})+\nu^{1}_{k}(\tau_{k,j}^{i})+\nu^{2}_{k}(\tau_{k,j}^{i}).\end{array}$
(4)

Variables $\nu^{i}_{k}=({\nu^{i}_{k,1}},\ldots,{\nu^{i}_{k,q_{i}(k)}})$,
$i=1,2,$ correspond to control impulses applied at the discrete time moments
$\tau_{k,1},\ldots,\tau_{k,q_{i}(k)}$ and represent the fraction of recovered
nodes. Let be $\nu^{i}_{k,j}={c^{i}_{k,j}}\delta(t-{\tau^{i}_{k,j}})$, where
$\delta(t-{\tau^{i}_{k,j}})$ is Dirac function,
${c^{i}_{k,j}}\in[0,{\overline{u}^{i}_{k,j}}]$ is the value of impulse, leads
to changes of the dynamical system, ${\overline{u}^{i}_{k,j}}$ is the maximum
value for control Agur .

Functional: the objective function of the combined system (4) is represented
by the aggregated costs on the time interval $[0,T]$ including the costs of
control impulses. The aggregated costs for continuous system (1) are defined
as follows: at time moment $t\neq{\tau^{i}_{k,j}},\ j=1,\ldots,q_{i}(k)$,
$i=1,2$, we have the costs from infected nodes $f^{1}_{k}(I^{1}_{k}(t))$ and
$f^{2}_{k}(I^{2}_{k}(t))$. Functions $f^{i}_{k}(\cdot)$ are non-decreasing and
twice-differentiable, such that $f^{i}_{k}(0)=0$, $f^{i}_{k}(I^{i}_{k}(t))>0$
for $I^{i}_{k}(t)>0$ with $t\in({\tau^{i}_{k,j-1}},{\tau^{i}_{k,j}}]$. For
system (4), we define the treatment costs as functions
$h^{i}_{k}(\nu^{i}_{k,j}({\tau^{i+}_{k,j}}))$, $j=1,\ldots,{q_{i}(k)}$, where
$h^{i}_{k}(\nu^{i}_{k,j}({\tau^{i+}_{k,j}}))>0$,
$\nu^{i}_{k,j}({\tau^{i+}_{k,j}})>0$ for $i=1,2$. Functions $g(R_{k}(t))$ are
non-decreasing and capture the benefit rates from Recovered nodes. The
aggregated system costs are defined by the functional:

$\begin{array}[]{l}J=\sum\limits_{k\in\mathbb{N}}[\int_{0}^{T}f^{1}_{k}(I^{1}_{k}(t))+f^{2}_{k}(I^{2}_{k}(t))-g(R(t))dt+\\\
\hskip
40.0pt\sum\limits_{j=1}^{q_{1}(k)}h^{1}_{k}(\nu^{1}_{k,j}({\tau^{1}_{k,j}}))+\sum\limits_{j=1}^{q_{2}(k)}h^{2}_{k}(\nu^{2}_{k,j}({\tau^{2}_{k,j}}))].\end{array}$
(5)

##  4 The structure of impulse control

According to principle maximum in impulse form Blaquiere , Chahim , Dykhta ,
Taynitskiy  we write Hamiltonian for dynamic system (1)

$\begin{array}[]{l}H^{0}_{k}(t)=-(f^{1}_{k}(I^{1}_{k}(t))+f^{2}_{k}(I^{2}_{k}(t))-g(R_{k}(t)))+(\lambda_{I^{1}_{k}}(t)-\lambda_{S_{k}}(t))\delta_{1k}S_{k}(t)I^{1}_{k}(t)\Theta_{1}(t)+\\\
\hskip
35.0pt(\lambda_{I^{2}_{k}}(t)-\lambda_{S_{k}}(t))\delta_{2k}S_{k}(t)I^{2}_{k}(t)\Theta_{2}(t)+(\lambda_{R_{k}}-\lambda_{I^{1}_{k}})\sigma^{1}_{k}I^{1}_{k}+(\lambda_{R_{k}}-\lambda_{I^{2}_{k}})\sigma^{2}_{k}I^{2}_{k};\end{array}$
(6)

and construct adjoint system as follows:

$\begin{array}[]{l}\displaystyle\dot{\lambda}_{S_{k}}(t)=(\lambda_{S_{k}}(t)-\lambda_{I^{1}_{k}}(t))\delta_{1k}I^{1}_{k}(t)\Theta_{1}(t)+(\lambda_{S_{k}}(t)-\lambda_{I^{2}_{k}}(t))\delta_{2k}I^{2}_{k}(t)\Theta_{2}(t);\\\
\displaystyle\dot{\lambda}_{I^{1}_{k}}(t)=\frac{df^{1}_{k}(I^{1}_{k}(t))}{dI^{1}_{k}}+(\lambda_{S_{k}}(t)-\lambda_{I^{1}_{k}}(t))\left(\delta_{1k}S_{k}(t)\Theta_{1}(t)+\frac{(\delta_{1k})^{2}S_{k}(t)I^{1}_{k}(t)P(k)}{\langle
k\rangle}\right)+\\\ \hskip
35.0pt(\lambda_{I^{1}_{k}}-\lambda_{R_{k}})\sigma^{1}_{k};\\\
\displaystyle\dot{\lambda}_{I^{2}_{k}}(t)=\frac{df^{2}_{k}(I^{2}_{k}(t))}{dI^{2}_{k}}+(\lambda_{S_{k}}(t)-\lambda_{I^{2}_{k}}(t))\left(\delta_{2k}S_{k}(t)\Theta_{2}(t)+\frac{(\delta_{2k})^{2}S_{k}(t)I^{2}_{k}(t)P(k)}{\langle
k\rangle}\right)+\\\ \hskip
35.0pt(\lambda_{I^{2}_{k}}-\lambda_{R_{k}})\sigma^{2}_{k};\\\
\displaystyle\dot{\lambda}_{R_{k}}(t)=-\frac{dg(R_{k}(t)}{dR_{k}},\end{array}$
(7)

with transversality conditions
$\lambda_{S_{k}}(T)=\lambda_{I^{1}_{k}}(T)=\lambda_{I^{2}_{k}}(T)=\lambda_{R_{k}}(T)=0$.

Following the maximum principle for impulse control (see Blaquiere , Sethi ,
Chahim ), we formulate necessary optimality conditions as in Theorem 4.1.

###### Theorem 4.1

Let $(x^{*},N,\tau_{i}^{j*},\nu_{i}^{*})$, $i=1,2,$ be an optimal solution for
the impulse control problem. Then, there exists an adjoint vector function
$\lambda(t)=(\lambda_{S}(t),\lambda_{I_{1}}(t),\lambda_{I_{2}}(t)$,
$\lambda_{R}(t))$ such that the following conditions hold:

$\begin{array}[]{l}\dot{\lambda}_{x}(t)=-\frac{\partial H_{0}}{\partial
x}(x^{*}(t),\lambda(t),t),\end{array}$ (8)

where $x(t)={S(t),I_{1}(t),I_{2}(t),R(t)}$.

At the impulse or jump points, it holds that

$\begin{array}[]{l}\frac{\partial
H_{i}^{c}}{\partial\nu_{i}}(x^{*}(\tau_{i}^{j*-}),\nu_{i},\lambda(\tau_{i}^{j*+}),\tau_{i}^{j*})(\nu_{i}^{j}-\nu_{i}^{j*})\geq
0,\\\ \\\ \end{array}$ (9)
$\begin{array}[]{l}\lambda_{x}(\tau_{i}^{j*+})-\lambda_{x}(\tau_{i}^{j*-})=\frac{\partial
H_{i}^{c}}{\partial
x}(x^{*}(\tau_{i}^{j*-}),\nu_{i}^{j*},\lambda(\tau_{i}^{j*+}),\tau_{i}^{j*}),\\\
\\\ \end{array}$ (10)
$\begin{array}[]{l}H_{0}(x^{*}(\tau_{i}^{j*+}),\lambda(\tau_{i}^{j*+}),\tau_{i}^{j*})-H_{0}(x^{*}(\tau_{i}^{j*-}),\lambda(\tau_{i}^{j*-}),\tau_{i}^{j*})-\vskip
12.0pt plus 4.0pt minus 4.0pt\\\ \hskip 1.0pt-\frac{\partial
H_{i}^{c}}{\partial\tau_{i}^{j}}(x^{*}(\tau_{i}^{j*-}),\nu_{i}^{j*},\lambda(\tau_{i}^{j*+}),\tau_{i}^{j*})\par\begin{cases}>0&\text{for
$\tau_{i}^{j*}=0$},\\\ =0&\text{for $\tau_{i}^{j*}\in(0,T)$,}\\\ <0&\text{for
$\tau_{i}^{j*}=T$.}\\\ \end{cases}\end{array}$ (11)

For all points in time at which there is no jump, i.e. $t\neq\tau_{j}\
(j=1,\ldots,k_{i})$, it holds that

$\begin{array}[]{l}\frac{\partial
H_{j}^{c}}{\partial\nu_{j}}(x^{*}(t),0,\lambda(t),t)\nu_{j}\leq 0,\end{array}$
(12)

with the transversality condition $\lambda(T)=0.$

Hamiltonian in impulsive form is

$\begin{array}[]{l}H^{c}_{k}(\tau_{k,j}^{1+})=-h^{1}_{k}(\nu^{1}_{k,j}(\tau_{k,j}^{1+}))+(\lambda_{R_{k}}(\tau_{k,j}^{1+})-\lambda_{I^{1}_{k}}(\tau_{k,j}^{1+}))\nu^{1}_{k,j}(\tau_{k,j}^{1+});\\\
H^{c}_{k}(\tau_{k,j}^{2+})=-h^{2}_{k}(\nu^{2}_{k,j}(\tau_{k,j}^{2+}))+(\lambda_{R_{k}}(\tau_{k,j}^{2+})-\lambda_{I^{2}_{k}}(\tau_{k,j}^{2+}))\nu^{2}_{k,j}(\tau_{k,j}^{2+}).\end{array}$
(13)

Here we assume that for each type of malwares $M_{1}$ and $M_{2}$ and for each
$k$ we have own set of control impulses
$\nu^{1}_{k}=({\nu^{1}_{k,1}},\ldots,{\nu^{1}_{k,q_{1}(k)}})$ and
$\nu^{2}_{k}=({\nu^{2}_{k,1}},\ldots,{\nu^{2}_{k,q_{2}(k)}})$.

Adjoin system for system(4) is ($i=1,2$):

$\begin{array}[]{l}\displaystyle\lambda_{S_{k}}(\tau_{k,j}^{i+})=\lambda_{S_{k}}(\tau_{k,j}^{i});\\\
\displaystyle\lambda_{I^{1}_{k}}(\tau_{k,j}^{i+})=\lambda_{I^{1}_{k}}(\tau_{k,j}^{i});\\\
\displaystyle\lambda_{I^{2}_{k}}(\tau_{k,j}^{i+})=\lambda_{I^{2}_{k}}(\tau_{k,j}^{i});\\\
\displaystyle\lambda_{R_{k}}(\tau_{k,j}^{i+})=\lambda_{R_{k}}(\tau_{k,j}^{i}).\end{array}$
(14)

Here is the conditions for $\Delta_{i}$ for each $I^{i}_{k}$ from the theorem
4.1:

$\begin{array}[]{l}\displaystyle\Delta_{1}=f^{1}_{k}(I^{1}_{k}(\tau_{k,j}^{1}))-f^{1}_{k}(I^{1}_{k}(\tau_{k,j}^{1+}))-g(R_{k}(\tau_{k,j}^{1}))+g(R_{k}(\tau_{k,j}^{1+}))+\\\
\hskip 23.0pt\displaystyle
c^{1}_{k,j}[\frac{dg(R_{k}(\tau_{k,j}^{1+}))}{dR_{k}(\tau_{k,j}^{1+})}+\frac{df^{1}_{k}(I^{1}_{k}(\tau_{k,j}^{1+}))}{dI^{1}_{k}(\tau_{k,j}^{1+})}]+2\sigma^{1}_{k}c^{1}_{k,j}(\lambda_{R_{k}}(\tau_{k,j}^{1+})-\lambda_{I^{1}_{k}}(\tau_{k,j}^{1+}))\\\
\hskip
23.0pt\displaystyle\delta_{1k}S_{k}(\tau_{k,j}^{1})c^{1}_{k,j}(\lambda_{S_{k}}(\tau_{k,j}^{1})-\lambda_{I^{1}_{k}}(\tau_{k,j}^{1}))[2\Theta_{1}(\tau_{k,j}^{1+})+\frac{\delta_{1k}P(k)}{\langle
k\rangle}(1+I^{1}_{k}(\tau_{k,j}^{1})-c^{1}_{k,j})].\end{array}$ (15)

Here is the conditions for $\Delta$ for each $\Delta$ for each $I^{2}_{k}$
from theorem 4.1:

$\begin{array}[]{l}\displaystyle\Delta_{2}=f^{2}_{k}(I^{2}_{k}(\tau_{k,j}^{2}))-f^{2}_{k}(I^{2}_{k}(\tau_{k,j}^{2+}))-g(R_{k}(\tau_{k,j}^{2}))+g(R_{k}(\tau_{k,j}^{2+}))+\\\
\hskip 23.0pt\displaystyle
c^{2}_{k,j}[\frac{dg(R_{k}(\tau_{k,j}^{2+}))}{dR_{k}(\tau_{k,j}^{2+})}+\frac{df^{2}_{k}(I^{2}_{k}(\tau_{k,j}^{2+}))}{dI^{2}_{k}(\tau_{k,j}^{2+})}]+2\sigma^{2}_{k}c^{2}_{k,j}(\lambda_{R_{k}}(\tau_{k,j}^{2+})-\lambda_{I^{2}_{k}}(\tau_{k,j}^{2+}));\\\
\hskip
23.0pt\displaystyle\delta_{2k}S_{k}(\tau_{k,j}^{2})c^{2}_{k,j}(\lambda_{S_{k}}(\tau_{k,j}^{2})-\lambda_{I^{2}_{k}}(\tau_{k,j}^{2}))[2\Theta_{2}(\tau_{k,j}^{2+})+\frac{\delta_{2k}P(k)}{\langle
k\rangle}(1+I^{2}_{k}(\tau_{k,j}^{2})-c^{2}_{k,j})].\end{array}$ (16)

According to Theorem 4.1 at time $\tau^{i}_{k,j}\in(0,T)$ $\Delta_{i}$ should
be equal to zero. Therefore, we deal with two different problems: firstly, if
the intensity of impulses $c^{i}_{k,j}$ are fixed, then from (15) and (16), we
can find the optimal time $\tau^{i*}_{k,j}$ of using impulses; secondly, if
the sequence of time $\tau^{i}_{k,j}$ are fixed, then we obtain the optimal
level of the intensity of impulses $c^{i*}_{k,j}$, $j=1,\ldots,q_{i}$,
$i=1,2$.

##  5 Numerical simulations

In this paragraph we present numerical experiments to depict theoretical
results and to study the behavior of malwares and show the application of
control impulses. Here we use the following set of the initial states and
values of parameters of the system (1): initial system states and parameters
are $S_{k}(0)=0.4,$ $I^{1}_{k}(0)=0.3,$ $I^{2}_{k}(0)=0.2$ and $R_{k}(0)=0$,
spreading rates are $\delta_{1k}=0.075k$, $\delta_{2k}=0.1k$, self-recovery
rates are $\sigma^{1}_{k}=0.0005k$ and $\sigma^{1}_{k}=0.0003k$. We set costs
functions for infectious subgroups as
$f^{i}_{k}(I^{i}_{k}(t))=A^{i}_{k}I^{i}_{k}(t)$ with coefficients
$A^{1}_{k}=2k$, $A^{2}_{k}=3k$ and treatment costs functions as
$h^{i}_{k}(\nu^{i}_{k,j}(\tau^{i+}_{k,j})=B^{i}_{k}c^{i}_{k,j}I^{i}_{k}(\tau^{i+}_{k,j})$,
where coefficients are equal to $B^{1}_{k}=3k$, $B^{2}_{k}=4k$,
$c^{1}_{k,j}=0.1$, $c^{2}_{k,j}=0.08$ for $i=1,2$, utility function is
$g(R_{k}(t))=0.1R_{k}(t)$.

Case 1. In case 1 we present the initial example of the behavior of the system
and aggregated system costs if an average number of links between $i$-th node
and its neighbors is $k=4$. Figs. 1-2 show the spreading on two modification
of malwares and corresponding total system costs.

Figure 1: Evolution of the system in Case 1. Number of links: $k=4$, spreading
rates: $\delta_{1k}=0.075k$ and $\delta_{2k}=0.1k$.  Figure 2: Aggregated
system costs are equal to $J=37.65$.

Aggregated system costs in this case are equal to $J=37.65$. By applying the
control impulses at discrete time moments we received that an amount of
impulses are equal to $p_{1}(4)=37$ and $p_{2}(4)=49$.

Case 2. In this experiment we use the same parameters for initial data, but in
contrast tot xase 1 an average number on neighbors is equal $k=7$. In this
case, we obtain that the aggregated costs are $J=73.93$, and an amount of
impulses are equal to $p_{1}(7)=29$ and $p_{2}(7)=44$. We may notice that
increasing the number of neighbor links increases the costs of the system.
Since there are less nodes with connectivity $k=7$ which is more than average
connectivity $\langle k\rangle=4$, we need less impulse treatment to vaccinate
the network, thereby if we apply control to more connected nodes we reduce the
costs of treatment.

Figure 3: 3 Evolution of the system in Case 2. Number of links: $k=7$,
spreading rates: $\delta_{1k}=0.075k$ and $\delta_{2k}=0.1k$. 3 Aggregated
system costs are equal to $J=73.93$.

Case 3. In case 3, by using the same initial set of data we variate the
spreading rate for malwares and consider $\delta_{1k}=0.075k$ and
$\delta_{2k}=0.1k$. Here we receive that the aggregated costs are $J=122.27$
and a number of impulses are $p_{1}(4)=43$ and $p_{2}(4)=55$, then increasing
the spreading rates are leading to increasing aggregated costs and number of
impulses which are needed to heal the network.

Figure 4: 4 Evolution of the system in Case 3. Number of links: $k=4$,
spreading rates: $\delta_{1k}=0.1k$ and $\delta_{2k}=0.2k$. 4 Aggregated
system costs are equal to $J=122.27$.

##  6 conclusion

This work addresses the spreading of cyber threats over large-scale networks
by investigating the optimal control policies in the impulsive form for SIR-
type of epidemics over scale-free networks. We have applied the impulse
optimal control framework to the epidemics over networks to devise impulsive
protection policies to mitigate the impact of the epidemics and contain the
spreading of the malware. With the application of the maximum principle, we
have obtained the structure of the optimal control impulses where actions are
taken at discrete-time moments. We have corroborated the obtained results
using numerical examples.

###### Acknowledgements.

The work of the second author was supported by the research grant “Optimal
Behavior in Conflict-Controlled Systems” (17-11-01079) from Russian Science
Foundation.

## References

  * (1) Kharraz, A., Robertson, W., Balzarotti, D., Bilge, L., Kirda, E,  Cutting the gordian knot: A look under the hood of ransomware attacks.  In _International Conference on Detection of Intrusions and Malware, and Vulnerability Assessment_ , pp. 3-24. Springer, 2015.
  * (2) Blaqui$\grave{e}$re A.  Impulsive Optimal Control with Finite or Infinite Time Horizon.  _Journal Of Optimization Theory And Applications._ , vol. 46. 1985.
  * (3) Chahim, M., Harti, R., Kort, P.  A tutorial on the deterministic Impulse Control Maximum Principle: Necessary and sufficient optimality conditions.  _European Journal of Operational Research._ , 219, 18–26, 2012.
  * (4) Dykhta V. A., Samsonyuk O. N.  A maximum principle for smooth optimal impulsive control problems with multipoint state constraints.  _Computational Mathematics and Mathematical Physics._ , 49, 942–957, 2009.
  * (5) Taynitskiy V., Gubar E., Zhu Q. Optimal Impulse Control of Bi-Virus SIR Epidemics with Application to Heterogeneous Internet of Things. _Constructive Nonsmooth Analysis and Related Topics. Abstracts of the International Conference. Dedicated to the Memory of Professor V.F. Demyanov_. p. 113-116, 2017.
  * (6) Barab si, A. L., Albert, R. Emergence of scaling in random networks. science, 286(5439), pp. 509-512, 1999.
  * (7) Gubar E., Zhu Q., Taynitskiy V. Optimal Control of Multi-strain Epidemic Processes in Complex Networks. _Game Theory for Networks. GameNets 2017. Lecture Notes of the Institute for Computer Sciences, Social Informatics and Telecommunications Engineering, vol 212_. Springer, Cham. p. 108–117, 2017.
  * (8) Gubar E., Kumacheva S., Zhitkova E., Porokhnyavaya O. Impact of Propagation Information in the Model of Tax Audit. _Recent Advances in Game Theory and Applications. ”Static and Dynamic Game Theory: Foundations and Applications”_. Switzerland. p. 91–110, 2015.
  * (9) Fu X., Small M., Walker D. M., Zhang H. Epidemic dynamics on scale-free networks with piecewise linear infectivity and immunization. _Phys. Rev. E. 77, 3, 036113_ , 2008.
  * (10) Pastor-Satorras R., Vespignani A. Epidemic spreading in scale-free networks. _Phys. Rev. Lett. 86, 14, 3200,_ 2001\.
  * (11) Agur Z., Cojocaru L., Mazor G., Anderson R. M., Danon Y. L. Pulse mass measles vaccination across age cohorts. _Proceedings of the National Academy of Sciences of the United States of America._ 90:11698–11702, 1993.
  * (12) Luo X, Liao Q. Ransomware: A new cyber hijacking threat to enterprises. _Handbook of research on information security and assurance._ 2009:1-6.
  * (13) Masuda N., Konno N. Multi-state epidemic processes on complex networks. _J. Theor. Biol.243, 1, 64–75_ , 2006.
  * (14) Sethi, S.P., Thompson, G.L.  Optimal Control Theory: Applications to Management Science and Economics.  Springer, Berlin, 2006.
  * (15) Taynitskiy V. A., Gubar E. A., Zhitkova E. M.  Optimization of protection of computer networks against malicious software, _Proc. of International Conference ”Stability and Oscillations of Nonlinear Control Systems” (Pyatnitskiy’s Conference)_ , 2016.
  * (16) Zaccour, G., Reddy, P., Wrzaczek, S.  Quality effects in different advertising models - An impulse control approach. _European Journal of Operational Research._ , 255, 984–995, 2016.
  * (17) Leander R., Lenhart S. and Vladimir Protopopescu V.  Optimal control of continuous systems with impulse controls.  _Optimal Control Applications and Methods Optim. Control Appl. Meth._ , 36:–549, 2015.
  * (18) Gubar, E., Zhu, Q.  Optimal Control of Influenza Epidemic Model with Virus Mutations.  _Proceedings 12th Biannual European Control Conference._ IEEE Control Systems Society. 3125–3130, 2013.


STFT

    short-time Fourier transform
iSTFT

    inverse short-time Fourier transform
DNN

    deep neural network
PESQ

    Perceptual Evaluation of Speech Quality
POLQA

    perceptual objectve listening quality analysis
WPE

    weighted prediction error
PSD

    power spectral density
RIR

    room impulse response
SNR

    signal-to-noise ratio
LSTM

    long short-term memory
POLQA

    Perceptual Objectve Listening Quality Analysis
SDR

    signal-to-distortion ratio
ESTOI

    extended short-term objective intelligibility
ELR

    early-to-late reverberation ratio
TCN

    temporal convolutional network
RLS

    recursive least squares
ASR

    automatic speech recognition
HA

    hearing aid
CI

    cochlear implant
MAC

    multiply-and-accumulate
VAE

    variational auto-encoder
GAN

    generative adversarial network
T-F

    time-frequency
SDE

    stochastic differential equation
DRR

    direct to reverberant ratio
LSD

    log spectral distance
SI-SDR

    scale-invariant signal to distortion ratio
MOS

    mean opinion score
AGC

    automatic gain control
VAD

    voice activity detection

Signal Processing (SP), Universität Hamburg, Germany

# Speech Signal Improvement Using Causal Generative Diffusion Models

###### Abstract

In this paper, we present a causal speech signal improvement system that is
designed to handle different types of distortions. The method is based on a
generative diffusion model which has been shown to work well in scenarios with
missing data and non-linear corruptions. To guarantee causal processing, we
modify the network architecture of our previous work and replace global
normalization with causal adaptive gain control. We generate diverse training
data containing a broad range of distortions. This work was performed in the
context of an “ICASSP Signal Processing Grand Challenge” and submitted to the
non-real-time track of the “Speech Signal Improvement Challenge 2023”, where
it was ranked fifth.

Index Terms—  Speech signal improvement, universal speech enhancement,
diffusion models, causal processing

##  1 Introduction

High-quality voice communication requires clear audio and natural-sounding
speech. However, there are numerous factors that can degrade speech signals
including background noise, room acoustics, transmission errors, limited
bandwidth, and codec artifacts. Prior works on improving speech signals have
typically studied each type of distortion separately. However, there has been
some recent interest in developing universal approaches that address a broader
range of distortions. These approaches typically use generative modeling,
which works particularly well in scenarios with missing data and non-linear
corruptions [1].

In this paper, we present our diffusion-based speech enhancement method
submitted to the non-real-time track of the “Speech Signal Improvement
Challenge 2023“ [2] as part of the “ICASSP Signal Processing Grand
Challenges“. The proposed model extends our previous work [3], incorporating
significant modifications in the network architecture to meet the causality
requirement and to output super wideband speech. Furthermore, we devise a data
corruption approach to generate diverse training data resembling distortions
observed in the blind data. Strong variations in loudness are compensated by
using causal adaptive gain control.

In the challenge’s subjective test, our proposed method yields a final score
of 0.445 using ITU-T P.863.2 and P.804 with mean opinion scores of 2.570
(Overall), 2.998 (Signal), 3.765 (Noise), 3.330 (Coloration), 3.674
(Discontinuity), 3.435 (Loudness), and 4.241 (Reverberation). In the
subjective test based on ITU-T P.835 our method yields a final score of 0.495
with MOSs of 2.842 (Overall), 3.119 (Signal), and 3.682 (Background). It is
interesting to note that, as our approach is generative in nature, it is
capable of achieving excellent performance for moderately distorted inputs,
while it may generate phonetic confusions and insertions if the input
distortions are too strong.

##  2 Proposed System

(a) Training

(b) Inference

Fig. 1: Proposed system: (a) At training, a corruption model $D$ generates $y$
from clean speech $x$. The forward diffusion moves from the clean spectrogram
$\mathbf{x}_{0}$ to the corrupted $\mathbf{y}$, while Gaussian noise is
gradually added. The score model $\mathbf{s}_{\theta}$ is learned using score
matching. (b) At inference, the data is normalized with causal adaptive gain
control (AGC) and the reverse diffusion uses the trained score model.

Fig. 1 shows an overview of our proposed system for speech signal improvement
at training and inference time. The diffusion process (Sec. 2.1) is the core
of the method and is accompanied by other processing blocks including the
short-time Fourier transform (STFT), causal automatic gain control (AGC) (Sec.
2.3), and a data corruption model $D$ used to simulate various distortion
types (Sec. 3.1).

###  2.1 Diffusion process

The diffusion process $\\{\mathbf{x}_{t}\\}_{t=0}^{T}$ for speech enhancement
is essentially the same as in our previous works [3, 4], where
$\mathbf{x}_{t}$ is the state of the process at time step $t$. During
training, the forward process moves from clean speech
$\mathbf{x}_{0}:=\mathbf{x}$ to corrupted speech $\mathbf{y}$, while
increasing amounts of Gaussian noise are gradually added. At inference, the
corresponding reverse process [5] is used to progressively remove the
corruption and therefore generate an estimate of $\mathbf{x}$ starting from
$\mathbf{x}_{T}\sim\mathcal{N}_{\mathbb{C}}(\mathbf{x}_{T};\mathbf{y},\sigma_{t}^{2})$.
This reverse process involves the _score function_ of $\mathbf{x}_{t}$, i.e.,
the gradient of its log-probability. Functioning as a prior for clean speech,
it is unavailable at inference time and is thus approximated by a trained deep
neural network (DNN) $\mathbf{s}_{\theta}$ called the _score model_. We make
no further adaptations to this process as it is defined for every STFT bin
independently, and is thus inherently causal as long as $\mathbf{s}_{\theta}$
is implemented with a causal network architecture.

###  2.2 Network architecture

We use a modified version of NCSN++ [5] for score estimation. The network is
an encoder-decoder architecture based on 2D convolutions, taking complex
spectrograms $\mathbf{x}_{t}$ and $\mathbf{y}$ and the process time $t$ as
input. Real and imaginary parts are considered as separate channels and the
convolutions are performed over time and frequency.

We apply the following modifications to the architecture to meet the causality
constraint: (1) Padding in the 2D convolutions is modified so that the
convolution along the time-dimension is causal; (2) Batch normalization is
replaced with cumulative group normalization, aggregating statistics
recursively; (3) Downsampling in the time dimension is performed with strided
convolutions and corresponding upsampling with transposed strided
convolutions. Up- and downsampling in the frequency dimension are realized
with finite impulse response filters, as in [3]. (4) All attention layers as
well as the progressive growing path are removed.

###  2.3 Automatic gain control

To match the unit-scale training condition of the score model, i.e.,
normalization of corrupted speech $y$ (see Fig. 1a), we use a causal AGC
system before feeding the mixture to the diffusion process, and again after
enhancement to maximize loudness, as part of the signal improvement task (see
Fig. 1b). To this end, we recursively track the maximum value per magnitude
frame averaged over the frequency bins to normalize the spectrogram in a
causal manner. We start tracking when speech activity is first detected, using
a voice activity detection method based on a causal speech presence
probability estimator [6]. The speech probability is fed through an ideal low-
pass filter to avoid having high-frequency noise bursts produce false
positives. Voice activity is then assumed if the speech presence probability
is higher than a threshold $\tau=0.8$, for a duration of $100$ms. When
discovering a larger maximum than the previous one, we smooth the
normalization with an exponential ramp going from the old value to the new
one. Finally, to avoid clipping, we use a causal compressor from the
pedalboard library111https://github.com/spotify/pedalboard.

##  3 Experimental Setup

###  3.1 Dataset

We use the VCTK corpus [7] as the clean speech dataset and resample all
utterances from 48kHz to our processing sampling frequency 32kHz. Using the
audiomentations library222https://github.com/iver56/audiomentations, we
simulate several corruptions observed in the blind data, namely stationary and
non-stationary noise, reverberation, clipping, gain reduction, packet loss and
lossy speech coding (the last two implemented by us). For each clean
utterance, a random corruption chain is chosen among plausible candidates,
e.g.
$\\{\mathrm{Reverb}\rightarrow\mathrm{Noise}\rightarrow\mathrm{PacketLoss}\\}$,
with the parameters of each corruption chosen randomly as well. We use the QUT
corpus [8] as the environmental noise dataset and take room impulse responses
from the DNS challenge [9] for reverberation.

###  3.2 Hyperparameters and training configuration

All processing is performed at $f_{s}=32$kHz and we upsample the processed
files to the original $48$kHz frequency. We use an STFT with a $638$-point
Hann window and $160$-point hop, which guarantees the global latency to be
below $20$ms, as we use purely causal processing. We set the lowest two
frequency bins of the spectrogram to zero, to remove the DC offset and low
frequency noise. The diffusion process hyperparameters are identical to those
in [3]. We train the score model $\mathbf{s}_{\theta}$ with the denoising
score matching objective [5] using the Adam optimizer with mini-batch size of
8 and exponential moving average over the parameters with a factor of $0.999$.
Training takes around two days on two NVIDIA RTX A$6000$ GPUs.

##  4 Results

In addition to the challenge’s subjective evaluation [2], we use DNSMOS P.835
[10] to evaluate our method on the 500 files in the blind set. Fig. 2 depicts
the histograms of DNSMOS scores: (a), (b) speech quality (SIG) and overall
quality (OVRL) of the improved files compared to the corrupted files; (c) SIG
with and without AGC.

Fig. 2: Speech signal improvement results obtained with DNSMOS.

Computational complexity: The network has about 55.7M parameters and the time
to infer a frame on a CPU (Intel Core i7-7800X @ 3.50GHz) takes 0.89s. Please
note that the model is designed to run on a GPU, on which the inference time
is orders of magnitude faster (0.02 s / frame with an NVIDIA GeForce RTX 2080
Ti).

##  5 Conclusion

In this work, we have built upon our previous work on diffusion-based speech
enhancement, with the novel contribution in making the system causal and
training the model on different distortion types. The proposed method was
ranked fifth in the non-real-time track of the “Speech Signal Improvement
Challenge 2023”.

##  6 References

## References

  * [1] Jean-Marie Lemercier, Julius Richter, Simon Welker and Timo Gerkmann “Analysing Diffusion-based Generative Approaches versus Discriminative Approaches for Speech Restoration” In _ICASSP_ , 2023
  * [2] Ross Cutler, Ando Saabas, Babak Naderi, Nicolae-Cătălin Ristea, Sebastian Braun and Solomiya Branets “ICASSP 2023 Speech Signal Improvement Challenge”, 2023 URL: https://arxiv.org/abs/2303.06566
  * [3] Julius Richter, Simon Welker, Jean-Marie Lemercier, Bunlong Lay and Timo Gerkmann “Speech enhancement and dereverberation with diffusion-based generative models” In _arXiv preprint_ , 2022
  * [4] Simon Welker, Julius Richter and Timo Gerkmann “Speech Enhancement with Score-Based Generative Models in the Complex STFT Domain” In _Interspeech_ , 2022
  * [5] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon and Ben Poole “Score-Based Generative Modeling through Stochastic Differential Equations” In _ICLR_ , 2021
  * [6] Timo Gerkmann and Richard C. Hendriks “Noise power estimation based on the probability of speech presence” In _WASPAA_ , 2011
  * [7] Junichi Yamagishi, Christophe Veaux and Kirsten MacDonald “CSTR VCTK Corpus: English multi-speaker corpus for CSTR voice cloning toolkit”, 2019
  * [8] David Dean, Sridha Sridharan, Robert Vogt and Michael Mason “The QUT-NOISE-TIMIT corpus for evaluation of voice activity detection algorithms” In _Interspeech_ , 2010
  * [9] Harishchandra Dubey, Vishak Gopal, Ross Cutler, Sergiy Matusevych, Sebastian Braun, Emre Sefik Eskimez, Manthan Thakker, Takuya Yoshioka, Hannes Gamper and Robert Aichner “ICASSP 2022 Deep Noise Suppression Challenge” In _ICASSP_ , 2022
  * [10] Chandan KA Reddy, Vishak Gopal and Ross Cutler “DNSMOS P.835: A non-intrusive perceptual objective speech quality metric to evaluate noise suppressors” In _ICASSP_ , 2022

  *[MOSs]: mean opinion score
  *[STFT]: short-time Fourier transform
  *[AGC]: automatic gain control
  *[DNN]: deep neural network


# Skyrmion size in skyrmion crystals

H. T. Wu  Physics Department, The Hong Kong University of Science and
Technology (HKUST), Clear Water Bay, Kowloon, Hong Kong  HKUST Shenzhen
Research Institute, Shenzhen 518057, China  X. C. Hu  Physics Department, The
Hong Kong University of Science and Technology (HKUST), Clear Water Bay,
Kowloon, Hong Kong  HKUST Shenzhen Research Institute, Shenzhen 518057, China
K. Y. Jing  Physics Department, The Hong Kong University of Science and
Technology (HKUST), Clear Water Bay, Kowloon, Hong Kong  HKUST Shenzhen
Research Institute, Shenzhen 518057, China  X. R. Wang<EMAIL_ADDRESS>Physics
Department, The Hong Kong University of Science and Technology (HKUST), Clear
Water Bay, Kowloon, Hong Kong  HKUST Shenzhen Research Institute, Shenzhen
518057, China

###### Abstract

A magnetic skyrmion is a topological object that can exist as a solitary
embedded in the vast ferromagnetic phase, or coexists with a group of its
“siblings” in various stripy phases as well as skyrmion crystals (SkXs).
Isolated skyrmions and skyrmions in an SkX are circular while a skyrmion in
other phases is a stripe of various forms. Unexpectedly, the sizes of the
three different types of skyrmions depend on material parameters differently.
For chiral magnetic films with exchange stiffness constant $A$, the
Dzyaloshinskii-Moriya interaction (DMI) strength $D$, and perpendicular
magnetic anisotropy $K$, $\kappa\equiv\pi^{2}D^{2}/(16AK)=1$ separates
isolated skyrmions from condensed skyrmion states. In contrast to isolated
skyrmions whose size increases with $D/K$ and is insensitive to $\kappa\ll 1$
and stripe skyrmions whose width increases with $A/D$ and is insensitive to
$\kappa\gg 1$, the size of skyrmions in SkXs is inversely proportional to the
square root of skyrmion number density and decreases with $A/D$. This finding
has important implications in our search for stable smaller skyrmions at the
room temperature in applications.

Skyrmions, localized topological objects characterized by the skyrmion number
$Q=\frac{1}{4\pi}\int{\mathbf{m}}\cdot(\partial_{x}{\mathbf{m}}\times\partial_{y}{\mathbf{m}}){\rm\,d}x{\rm
d}y$, were originally proposed as resonance states of baryons [1], but were
first unambiguously observed in chiral magnets that involve the
Dzyaloshinskii-Moriya interaction (DMI) [2, 3, 4, 5, 6, 9, 7, 8, 10, 12, 14,
11, 13]. Here $\mathbf{m}$ is the unit vector of magnetization and $Q$
measures the number of times that spins, after moving to the origin, cover the
unit sphere of $|\mathbf{m}|=1$. Skyrmions have received intensive and
extensive studies [2, 3, 4, 5, 6, 9, 7, 8, 10, 12, 14, 15, 16, 17, 10, 18, 7,
34, 19, 21, 11, 13, 25, 26, 28, 35, 29, 8, 31, 24, 26, 22, 23, 30, 27, 33, 7,
34, 32, 20, 32, 31] in magnetics for their academic interest and promising
applications in information technology.

Our knowledge about magnetic skyrmions has been greatly advanced in recent
years [15, 16, 17]. Skyrmions in metastable states are fundamentally different
from the ones in the ground state. A skyrmion is circular and its formation
energy is positive when the skyrmion is metastable and ferromagnetic ordering
is the ground state. Skyrmion morphologies are various stripes that could be
ramified or non-ramified, long or short, and straight or curved when skyrmions
are the ground state. Thus, present knowledge of skyrmion zoo consisting of
three types of circular skyrmions [15], namely Bloch skyrmions (whirl
skyrmions), the Néel skyrmions (hedgehog skyrmions), and anti-skyrmions,
should be greatly expanded by including various stripes in disordered stripy
phases, in helical states, and in maze. Even when skyrmions are the ground
state, individual skyrmion can still be disk-like in triangular lattice of
closely packed skyrmions where skyrmion-skyrmion repulsion dominates the
morphology of individual skyrmions. It is already known [22] that the size of
isolated skyrmions increases with the DMI strength, and decreases with the
exchange stiffness coefficient and magnetic anisotropy [22, 23]. It is also
known [16] that the width of stripe skyrmions increases with exchange
stiffness coefficient and decreases with the DMI strength, but does not depend
on skyrmion number density. Since the shape of a skyrmion in a skyrmion
crystal (SkX) has different underlying physics as an isolated skyrmion and a
stripe skyrmion, the parameter dependence of the size of disk-like skyrmions
in an SkX may also be different although people do not distinguish a skyrmion
in an SkX from an isolated ones to date.

In this paper, we use a generic chiral magnetic film with exchange stiffness
constant $A$, the DMI strength $D$, and perpendicular magnetic anisotropy $K$
to show that $\kappa\equiv\pi^{2}D^{2}/(16AK)<1$ and $\kappa>1$ determines the
existence of isolated skyrmions (the former) and stripe skyrmions or SkXs (the
latter). Although skyrmions in SkXs and stripe skyrmions are “siblings”, the
dominate repulsion among skyrmions in SkX changes the parameter dependence of
skyrmion size. This finding provides a general guidance in skyrmion size
manipulation.

We consider a two-dimensional (2D) ferromagnetic film of thickness $d$ in the
$xy$-plane described by the following energy

$\displaystyle E=$ $\displaystyle
d\iint\\{A|\nabla\mathbf{m}|^{2}+K(1-m_{z}^{2})$ (1)
$\displaystyle+D[(\mathbf{m}\cdot\nabla)m_{z}-m_{z}\nabla\cdot\mathbf{m}]\\}\mathrm{\,d}S,$

where $\mu_{0}$ is the vacuum permeability. The energy is set to zero, $E=0$,
for ferromagnetic state of $|m_{z}|=1$. Films with $\kappa<1$ support isolated
circular skyrmions that are metastable state of energy $8\pi
Ad\sqrt{1-\kappa}$ [22]. The dynamics of $\mathbf{m}$ is governed by the
Landau-Lifshitz-Gilbert (LLG) equation,

$\frac{\partial\mathbf{m}}{\partial t}=-\gamma\mathbf{m}\times\mathbf{H}_{\rm
eff}+\alpha\mathbf{m}\times\frac{\partial\mathbf{m}}{\partial t},$ (2)

where $\gamma$ and $\alpha$ are respectively gyromagnetic ratio and the
Gilbert damping constant. The effective magnetic field $\mathbf{H}_{\rm
eff}=2A\nabla^{2}\mathbf{m}+2K_{1}m_{z}\hat{z}+\mathbf{H}_{\rm\,d}+\mathbf{H}_{\rm
DM}$ includes the exchange field, the magneto-crystalline anisotropy field of
coefficient $K_{1}$, the external magnetic field, the demagnetizing field
$\mathbf{H}_{\rm\,d}$, and the DMI field $\mathbf{H}_{\rm DM}$ respectively.
$\mathbf{H}_{\rm\,d}$ can be included in the effective anisotropy of
$K=K_{1}-\mu_{0}M_{\rm s}^{2}/2$. This approximation is good when $d$ is much
smaller than the exchange length [22].

In the absence of energy source such as an electric current and at zero
temperature, the LLG equation describes a dissipative system whose energy can
only decrease [37, 38]. Thus, the long-time solution of the LLG equation with
any initial configuration must be a stable spin texture of Hamiltonian Eq.
(1). In this study, we choose three groups of model parameters to simulate
those chiral magnets such as MnSi with $A=0.27\sim 0.4\,\hbox{pJ/m}$,
$D=0.15\sim 0.33\,\hbox{mJ/m}^{2}$, $K=0.02\sim 0.2\,\hbox{MJ/m}^{3}$, $M_{\rm
s}=0.15\sim 0.5\,\hbox{MA/m}$ [19]; Co-Zn-Mn with $A=5\sim 11\,\hbox{pJ/m}$,
$D=0.5\sim 2\,\hbox{mJ/m}^{2}$, $K=20\sim 80\,\hbox{kJ/m}^{3}$, $M_{\rm
s}=0.15\sim 0.35\,\hbox{MA/m}$ [27]; and PdFe/Ir and W/Co20Fe60B20/MgO [10,
18, 7, 34, 26, 21] with $A=2\sim 10\,\hbox{pJ/m}$, $D=0.68\sim
4\,\hbox{mJ/m}^{2}$, $K=0.00228\sim 5\,\hbox{MJ/m}^{3}$, $M_{\rm s}=650\sim
1160\,\hbox{kA/m}$. Periodic boundary conditions are used to eliminate
boundary effects and the MuMax3 package [36] is employed to numerically solve
the LLG equation with mesh size of $1\,{\rm nm}\times 1\,{\rm nm}\times
1\,{\rm nm}$. The number of stable states and their structures should not
depend on the Gilbert damping constant. We use a large $\alpha=0.25$ to speed
up our simulations.

Figure 1: a) Energy (in ${\rm arcsinh}(E)$ scale) as a function of $\kappa$
(in $\log(\kappa)$ scale) for various model parameters with very different
$A$, $D$ and $K$. $E$ changes sign at $\kappa=1$ that separate isolated
skyrmions from stripe skyrmions and skyrmions in SkXs. (b) Spin profile
$m_{z}(r)$ of an isolated skyrmion shown in the inset with parameters denoted
by $b$ in (a). Symbols are simulation data and curve are theoretical fit with
$R=4.513\,$nm and $w=1.356\,$nm. (c) Spin profile $m_{z}(x)$ of stripe
skyrmion with parameters denoted by $c$ in (a). Symbols are simulation data
for stripes ①-④ in the inset and the solid curve is the theoretical fit with
$L=8.95\,$nm and $w=1.82\,$nm. $x$ is measured across the stripe and $x=0$ is
the stripe centre.

Our previous studies [16, 17] demonstrate that a small nucleation domain of
$m_{z}=-1$ in the background of ferromagnetic phase of $m_{z}=1$ can develop
into a skyrmion in a chiral magnetic film. Using MuMax3 [36], we obtain the
energy of one static skyrmion. The sign of $E$ tells us whether the skyrmion
is a metastable state or a ground state because $E=0$ is set for the
ferromagnetic state of $|m_{z}|=1$. Figure 1(a) plots the energy of one
skyrmion in chiral magnets as a function of $\kappa$. The squares and
diamonds, circles and stars, up-triangles and pentagons, are respectively from
three groups of model parameters listed above. Because the huge range of $E$
is covered around zero, $E$ is plotted in ${\rm arcsinh}(E)$ scale and
$\kappa$ is in $\log(\kappa)$ scale. Clearly, $\kappa=1$ separate metastable
circular isolated skyrmions from stripe skyrmions.

The spin profile of isolated skyrmions can be approximated by
$\Theta(r)=2\arctan\left[\frac{\sinh(r/w)}{\sinh(R/w)}\right]$ [22]. $\Theta$
is the polar angle of the magnetization at position $r$ measured from the
skyrmion centre. $R$ and $w$ measure respectively the skyrmion size and
skyrmion wall thickness. Figure 1(b) shows the comparison of spin profile from
MuMax3 (symbols) and theoretical fit with $R=4.513\,$nm and $w=1.356\,$nm for
an isolated skyrmion with $A=0.4\,\hbox{pJ/m}$, $D=0.15\,\hbox{mJ/m}^{2}$,
$K=0.2\,\hbox{MJ/m}^{3}$, $M_{\rm s}=0.15\,\hbox{MA/m}$. The inset shows the
skyrmion structure. This profile leads to the skyrmion size formula of
$R=w/\sqrt{1-\kappa}$ and $w=\pi D/(4K)$ [22].

Figure 2: SkXs and skyrmion spin profiles at various skyrmion number density
$2.34\times 10^{-3}\hbox{/nm}^{2}$ (a, d), $2.89\times 10^{-3}\hbox{/nm}^{2}$
(b, e), and $3.50\times 10^{-3}\hbox{/nm}^{2}$ (c, f). (d-f) Spin profiles
$m_{z}(r)$ along the red lines in (a), (b), and (c), respectively. The symbols
are numerical data and the solid lines are the theoretical fits with
$R=7.042\,$nm (d), $6.144\,$nm (e), and $5.399\,$nm (f).

The spin profile of stripe skyrmions can be approximated by
$\Theta(x)=2\arctan\left[\frac{\sinh(L/2w)}{\sinh(|x|/w)}\right]$ for
$m_{z}<0$ and
$\Theta(x)=2\arctan\left[\frac{\sinh(|x|/w)}{\sinh(L/2w)}\right]$ for
$m_{z}>0$ [16]. $x$ is measured across the stripe and $x=0$ is at stripe
centre. $L$ and $w$ measure respectively the stripe width and skyrmion wall
thickness. Figure 1(c) is spin profile of stripe skyrmion denoted by $c$ in
Fig. 1(a) with model parameter of $A=0.3\,\hbox{pJ/m}$,
$D=0.24\,\hbox{mJ/m}^{2}$, $K=0.096\,\hbox{MJ/m}^{3}$, $M_{\rm
s}=0.15\,\hbox{MA/m}$. The symbols are the numerical data from various
locations of stripe skyrmion shown in the inset and the solid curve is the
theoretical fit with $L=8.95$ nm and $w=1.82$ nm. The excellent agreements
between the numerical data and the approximate profiles demonstrate that one
can use this theoretical profile to extract the skyrmion size accurately. All
data from different stripes fall on the same curve demonstrate that stripes,
building blocks of irregular skyrmions, are identical. The excellent stripe
skyrmion profile leads to the stripe width formula $L=g(\kappa)A/D$, where
$g(\kappa)$ is about 6.23 for $\kappa\gg 1$ [16], opposite dependence on $D$
as that for an isolated skyrmion.

In order to study skyrmions in SkXs, we consider model parameters satisfying
$\kappa>1$ below. All initial configurations have sufficient number of
nucleation domains in triangular lattices. Each domain of $m_{z}=-1$ is a disk
of 10nm in diameter. Each initial configuration leads to an SkX as shown in
Figs. 2(a)-(c) with respectively 324, 400, and 484 skyrmions in a film of
$400\,\hbox{nm}\times 346\,\hbox{nm}\times 1\,\hbox{nm}$. The model parameters
are $A=5\,\hbox{pJ/m}$, $D=3\,\hbox{mJ/m}^{2}$, $K=4.87\,\hbox{kJ/m}^{3}$,
$M_{\rm s}=0.2\,\hbox{MA/m}$. Figs. 2(a)-(c) show that the territories of
$m_{z}>0$ and $m_{z}<0$ are different, and have different spin distribution.
Lattice constant and skyrmion size are two obvious length scales. Spin
profiles of isolated skyrmions and stripe skyrmions should not work for SkXs
because of different underlying physics, a new spin polar angle distribution
of $\Theta(r)=2\arctan\left\\{\frac{\tan\left[\frac{\pi}{2}\cos(\frac{\pi
R}{L})\right]}{\tan\left[\frac{\pi}{2}\cos(\frac{\pi r}{L})\right]}\right\\}$
along the three triangular lattice directions, (1,0), (1/2, $\sqrt{3}/2$), and
(1/2, $-\sqrt{3}/2$), is proposed, where $r$ labels the points on the yellow
lines in (a), (b), and (c) with $r=0$ being the centre of one chosen skyrmion.
$R$ and $L$ are respectively the skyrmion size (the radius of $m_{z}=0$
contour) and lattice constant. $L$ relates to skyrmion number density $n$ in a
triangular SkX as $L\approx 1.07/\sqrt{n}=22.11\,$nm (a), $19.90\,$nm (b), and
$18.09\,$nm (c). Figure 2(d), (e), and (f) are $m_{z}(r)$ along the yellow
lines in (a), (b), and (c), respectively. Symbols are numerical data from
MuMax3 while the curves are the fitting of $m_{z}(r)=\cos\Theta$ with the
fitting parameter $R=7.042\,$nm (d), $6.144\,$nm (e), $5.399\,$nm (f). One of
the interesting observations is that $R$ depends only on $n$, $A/D$, and
$\kappa$. Numerical evidences on the result can be found in the Supplementary
Materials.

If we assume energy per skyrmion in an SkX is from a skyrmion whose spin
texture is described by $\Theta(r)$ for $0\leq r\leq L/2$, then energy per
skyrmion can be expressed as a function of variable $R$,

$E(R)=2\pi
d[Af_{1}(\varepsilon)+DLf_{2}(\varepsilon)+KL^{2}f_{3}(\varepsilon)],$ (3)

where $\varepsilon=\frac{2R}{L}$,
$f_{1}(\varepsilon)=4f^{2}(\varepsilon){\displaystyle\int_{0}^{1}}\frac{\left[\xi^{2}f^{\prime
2}(\varepsilon)+f^{2}(\xi)\right]}{\xi\left[f^{2}(\xi)+f^{2}(\varepsilon)\right]^{2}}{\rm\,d}\xi$,
$f_{2}(\varepsilon)=f(\varepsilon){\displaystyle\int_{0}^{1}}\left\\{\frac{\xi
f^{\prime}(\xi)}{f^{2}(\xi)+f^{2}(\varepsilon)}+\frac{f(\xi)\left[f^{2}(\varepsilon)-f^{2}(\xi)\right]}{\left[f^{2}(\xi)+f^{2}(\varepsilon)\right]^{2}}\right\\}{\rm\,d}\xi$,
$f_{3}(\varepsilon)={\displaystyle\int_{0}^{1}}\left[\frac{f(\xi)f(\epsilon)}{f^{2}(\xi)+f^{2}(\varepsilon)}\right]^{2}{\rm\,d}\xi$,
and $f(x)=\tan\left[\frac{\pi}{2}\cos(\frac{\pi}{2}x)\right]$. The skyrmion
size is then the solution of $\frac{{\rm\,d}E}{{\rm\,d}R}=0$. Under certain
reasonable assumptions (see the Supplementary Materials), one has

$R(n,A/D,\kappa)=\frac{0.54}{\sqrt{n}}-\sqrt{\frac{0.51}{\sqrt{n}\frac{A}{D}+0.56n\frac{A^{2}}{D^{2}}-0.02\frac{1}{\kappa}}}\frac{A}{D}.$
(4)

Unlike the width of stripe skyrmions that does not depend on $n$ and increase
with $A/D$, Eq. (4) says that $R$ is almost inversely proportional to
$\sqrt{n}$ and decreases with $A/D$. The interesting dependences agree with
the underlying physics of SkX that is the balance of skyrmion-skyrmion
repulsion and intrinsic stripe nature of skyrmions. $R$ should then be
proportional to the lattice constant of a SkX that is the skyrmion-skyrmion
distance and, in turn, is inversely proportional to $\sqrt{n}$. Since
skyrmions in a SkX come from the squeeze of stripe skyrmions, the amount of
width deduction from the squeeze should be proportional to the original stripe
width proportional to $A/D$. Thus, $R$ should not only be smaller than the
natural stripe width, but also decreases with $A/D$ as what Eq. (4) says. Of
course, one should not expect coefficients of $1/\sqrt{n}$ and $A/D$ to be
very accurate because the obvious approximations involved that include
impossible filling by disks of radius $L/2$ and the neglect of energy from
uncovered region in Eq. (3).

We carried out MuMax3 simulations with very different $n$, $A$, $D$, $K$ in
three different parameter regions. For each fixed parameters, skyrmion size
$R$ in the corresponding SkX can be obtained either from the size of $m_{z}=0$
contour or through fitting of numerical spin profile to $\Theta(r)$. The
difference between two extracted values is negligible (see the Supplementary
Materials), and all data shown below is from contour measurement. One can test
our theoretical predictions of Eq. (4) against numerical simulations. As
mentioned above, simulations show that $R$ depends only on $n$, $A/D$, and
$\kappa$. Fig. 3(a) plots $R$ vs $1/\sqrt{n}$ for various fixed $A/D$ and
$\kappa$ from three groups model parameters. Simulation results (symbols) are
well described by our formula (solid curves) of Eq. (4) without any fitting.
$R$ is almost linear in $1/\sqrt{n}$. The $A/D$ dependence and $\kappa$
dependence of $R$ can also be captured well by Eq. (4) as shown by good
agreement between simulation data (symbols) and theoretical formula of Eq. (4)
(the solid curve) in Fig. 3(b) ($A/D$) and Fig. 3(c) ($\kappa$), respectively.
Interestingly, Eq. (4) says that $R$ does not depend on anisotropy $K$ for
$\kappa\gg 1$. Our new spin profile captures successfully the skyrmion-
skyrmion repulsion effect that requires $R$ to be proportional to the lattice
constant, the inverse of the skyrmion number density $n$.

Figure 3: Skyrmion size in SkXs as a function of skyrmion number density $n$
for various $A/D$ and fixed $\kappa$ (a), as a function of $A/D$ for various
$n$ and fixed $\kappa$ (b), and as a function of $\kappa$ for various $n$ and
$A/D$. The symbols are numerical data from MuMax3. The solid lines are the
formula of Eq. (4).

Magnetic skyrmions can be isolated, or condensed regular or irregular stripes
distributed randomly or arranged as helical states of one-dimensional stripe
lattices, or condensed circular skyrmions in crystal forms. The irregular
stripes can even be ramified in dendrite or maze forms. The underlying physics
of skyrmion size are different for isolated skyrmions, skyrmions in an SkX and
stripe skyrmions in vastly different morphologies and structures. It is known
that skyrmion formation energy are positive for isolated skyrmions and
negative for stripe skyrmions and SkXs. Positive formation energy results in
an extra energy for skyrmion surfaces, similar to the surface tension of a
liquid droplet, that favours a circular shape. Negative skyrmion formation
energy is similar to negative surface tension of a liquid droplet that prefers
to wet its contacts and spread out. From this viewpoint, it may not be
surprising to see different parameter dependences of skyrmion size for
isolated skyrmions and skyrmions in SkXs or stripe skyrmions. It is also not
surprising that spin profile of skyrmions in SkXs has little similarity as
those of isolated and stripe skyrmions because their shape comes from the
balance of the skyrmion squeezing and their stripe nature. Two nearby
skyrmions start to repel each other when their distance is order of $2L\approx
13A/D$. Thus, the critical skyrmion number density for SkX is about
$n_{c}\approx D^{2}/(169A^{2})$. For $n>n_{c}$, skyrmions tend to form closely
packed triangular lattice with lattice constant of $1.07/\sqrt{n}$. Due to the
squeezing of skyrmions, the skyrmion size (diameter of $m_{z}=0$ contour)
should be smaller than $0.53/\sqrt{n}$, this squeezing effect leads to Eq. 4
in which skyrmion size decreases with $A/D$, opposite to the parameter
dependence of the width of stripe skyrmions.

In applications, creating stable skyrmions of smaller size at the room
temperature is a goal. According to current results, one needs to distinguish
an isolated skyrmion ($\kappa<1$) from a skyrmion in condensed skyrmion phase
($\kappa>1$). In the case of $\kappa<1$, materials with large $A$ and smaller
$D$ are preferred since larger $A$ and smaller $D$ imply a higher $T_{\rm c}$
and smaller skyrmion size. On the other hand, in the case of $\kappa>1$,
materials with larger $A$ and larger $D$ are required so that one may have
both higher $T_{\rm c}$ and smaller skyrmion size.

In conclusion, we clarify that there are three distinct skyrmions, namely
circular isolated skyrmions, stripe skyrmions, and disk-like skyrmions in
SkXs. The material parameter dependences of skyrmion size are very different.
Unlike an isolated skyrmion whose size increases with $D/K$ and $\kappa<1$, or
a stripe skyrmion whose size increases with $A/D$ and decreases with
$\kappa>1$, the size of a disk-like skyrmion in SkXs is inversely proportional
to the skyrmion number density and decreases with $A/D$. Skyrmion size is an
very important quantity, and high-density data storage require smaller
skyrmions in applications. The findings here provide a practical guidance to
search proper materials in different applications.

###### Acknowledgements.

This work is supported by National Key R&D Program of China, grant No
2020YFA0309600, the NSFC Grant (No. 11974296 and 11774296) and Hong Kong RGC
Grants (No. 16301518 and 16301619). Partial support by the National Key
Research and Development Program of China (Grant No. 2018YFB0407600) is also
acknowledged.

## References

  * [1] T. H. R. Skyrme, Nucl. Phys. 31, 556-569 (1962).
  * [2] U. K. Rößler, A. N. Bogdanov, and C. Pfleiderer, Nature 442, 797-801 (2006).
  * [3] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfeiderer, A. Rosch, A. Neubauer, G. Georgii, and P. Böni, Science 323, 915-919 (2009).
  * [4] X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, Nature 465, 901–904 (2010).
  * [5] Y. Onose, Y. Okamura, S. Seki, S. Ishiwata, and Y. Tokura, Nature 465, 901–904 (2010).
  * [6] H. S. Park, X. Z. Yu, S. Aizawa, T. Tanikaki, T. Akashi, Y. Takahashi, T. Matsuda, N. Kanazawa, Y. Onose, D. Shindo, A. Tonomura, and Y. Tokura, Nat. Nanotech. 9, 337 (2014).
  * [7] N. Romming, C. Hanneken, M. Menzel, J. E. Bickel, B. Wolter, K. von Bergmann, A. Kubetzka, and R. Wiesendanger, Science 341, 636-639 (2013).
  * [8] J. Li, A. Tan, K. W. Moon, A. Doran, M. A. Marcus, A. T. Young, E. Arenholz, S. Ma, R. F. Yang, C. Hwang, and Z. Q. Qiu, Nat. Commun. 5, 4704 (2014).
  * [9] S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, and Stefan Blügel, Nat. Phys. 7, 713-718 (2011).
  * [10] W. Jiang, P. Upadhyaya, W. Zhang, G. Yu, M. B. Jungfleisch, F. Y. Fradin, J. E. Pearson, Y. Tserkovnyak, K. L. Wang, O. Heinonen et al., Science 349, 283 (2015).
  * [11] H. Du, R. Che, L. Kong, X. Zhao, C. Jin, C. Wang, J. Yang, W. Ning, R. Li, C. Jin et al., Nat. Commun. 6, 8504 (2015).
  * [12] S. Krause and R. Wiesendanger, Nat. Mater. 15, 493-494 (2016).
  * [13] H. Y. Yuan and X. R. Wang, Sci. Rep. 6, 22638 (2016).
  * [14] S. Woo, K. Litzius, B. Krüger, M.-Y. Im, L. Caretta, K. Richter, M. Mann, A. Krone, R. M. Reeve, M. Weigand et al., Nat. Mater. 15, 501-506 (2016).
  * [15] C. Back, V. Cros, H. Ebert, K. Everschor-Sitte, A. Fert, M. Garst, Tianping Ma, S. Mankovsky, T. L. Monchesky, M. Mostovoy et al., J. Phys. D 53, 363001 (2020).
  * [16] X. R. Wang, X. C. Hu, and H. T. Wu, arXiv:2101.10561
  * [17] X. C. Hu, H. T. Wu, and X. R. Wang, arXiv:2101.10104
  * [18] S. Jaiswal, K. Litzius, I. Lemesh, F. Büttner, S. Finizio, J. Raabe, M. Weigand, K. Lee, J. Langer, B. Ocker, G. Jakob, G. S. D. Beach, and M. Kläui, Appl. Phys. Lett. 111, 022409 (2017).
  * [19] E. A. Karhu, U. K. Rößler, A. N. Bogdanov, S. Kahwaji, B. J. Kirby, H. Fritzsche, M. D. Robertson, C. F. Majkrzak, and T. L. Monchesky, Phys. Rev. B 85, 094429 (2012).
  * [20] S. Rohart, and A. Thiaville, Phys. Rev. B 88, 184422 (2013).
  * [21] E. Simon, K. Palotás, L. Rózsa, L. Udvardi, and L. Szunyogh, Phys. Rev. B 90, 094410 (2014).
  * [22] X. S. Wang, H. Y. Yuan, and X. R. Wang, Commun. Phys. 1, 31 (2018).
  * [23] D. Cortés-Ortuño, N. Romming, M. Beg, K. von Bergmann, A. Kubetzka, O. Hovorka, H. Fangohr, and R. Wiesendanger, Phys. Rev. B 99, 214408 (2019).
  * [24] J. Iwasaki, M. Mochizuki, and N. Nagaosa, Nat. Commun. 4, 1463 (2013).
  * [25] N. Nagaosa, and Y. Tokura, Nat. Nanotech. 8, 899 (2013).
  * [26] A. Fert, V. Cros, and J. Sampaio, Nat. Nanotech. 8, 152 (2013); J. Sampaio, V. Cros, S. Rohart, A. Thiaville, and A. Fert, Nat. Nanotechnol. 8, 839 (2013).
  * [27] V. Ukleev, Y. Yamasaki, D. Morikawa, K. Karube, K. Shibata, Y. Tokunaga, Y. Okamura, K. Amemiya, M. Valvidares, H. Nakao et al., Phys. Rev. B. 99, 144408 (2019).
  * [28] Y. Zhou, and M. A. Ezawa, Nat. Commun. 5, 4652 (2014).
  * [29] P. Dürrenfeld, Y. Xu, J. Åkerman, and Y. Zhou, Phys. Rev. B 96, 054430 (2017).
  * [30] L. Kong amd J. Zang, Phys. Rev. Lett. 111, 067203 (2013).
  * [31] H. Y. Yuan, X. S. Wang, Man-Hong Yung, and X. R. Wang, Phys. Rev. B 99, 014428 (2019).
  * [32] X. Gong, H. Y. Yuan, and X. R. Wang, Current-driven skyrmion motion in granular films, Phys. Rev. B 101, 064421 (2020).
  * [33] M. N. Wilson, A. B. Butenko, A. N. Bogdanov, and T. L. Monchesky, Phys. Rev. B 89, 094411 (2014).
  * [34] N. Romming, A. Kubetzka, C. Hanneken, K. von Bergmann, and R. Wiesendanger, Phys. Rev. Lett. 114, 177203 (2015).
  * [35] A. Siemens, Y. Zhang, J. Hagemeister, E. Y. Vedmedenko, and R Wiesendanger, New. J. Phys. 18, 045021 (2016).
  * [36] A. Vansteenkiste, J Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. van Waeyenberge, AIP. Adv. 4, 107133 (2014).
  * [37] X. R. Wang, P. Yan, J. Lu, and C. He, Ann. Phys. (NY) 324, 1815 (2009).
  * [38] X. R. Wang, P. Yan, and J. Lu, Europhys. Lett. 86, 67001 (2009).