# Confidence and Backaction in the Quantum Filter Equation

Wei Cui<EMAIL_ADDRESS>Advanced Science Institute, RIKEN, Wako-shi, Saitama
351-0198, Japan  Neill Lambert<EMAIL_ADDRESS>Advanced Science Institute,
RIKEN, Wako-shi, Saitama 351-0198, Japan  Yukihiro Ota  Advanced Science
Institute, RIKEN, Wako-shi, Saitama 351-0198, Japan  Xin-You Lü  Advanced
Science Institute, RIKEN, Wako-shi, Saitama 351-0198, Japan  Z. -L. Xiang
Advanced Science Institute, RIKEN, Wako-shi, Saitama 351-0198, Japan
Department of Physics, State Key Laboratory of Surface Physics, Fudan
University, Shanghai 200433, China  J. Q. You  Department of Physics, State
Key Laboratory of Surface Physics, Fudan University, Shanghai 200433, China
Beijing Computational Science Research Center, Beijing 100084, China  Advanced
Science Institute, RIKEN, Wako-shi, Saitama 351-0198, Japan  Franco Nori
Advanced Science Institute, RIKEN, Wako-shi, Saitama 351-0198, Japan  Physics
Department, The University of Michigan, Ann Arbor, Michigan 48109-1040, USA

###### Abstract

We study the confidence and backaction of state reconstruction based on a
continuous weak measurement and the quantum filter equation. As a physical
example we use the traditional model of a double quantum dot being
continuously monitored by a quantum point contact. We examine the confidence
of the estimate of a state constructed from the measurement record, and the
effect of backaction of that measurement on that state. Finally, in the case
of general measurements we show that using the relative entropy as a measure
of confidence allows us to define the lower bound on the confidence as a type
of quantum discord.

###### pacs:

03.67.-a, 03.65.Ta, 03.65.Yz, 73.63.Kv

##  I Introduction

Recently there has been a great deal of activity on the topic of “weak”
quantum measurements Wiseman2010 ; Alicki05 ; Kofman11 ; Ashhab ; Ashhab09 ;
Croke06 ; Braginsky  in both mesoscopic Korotkov99 ; Korotkov01 ; Korotkov06
and macroscopic systems Buluta11 ; Macroscopic ; You05 ; Johansson . In
contrast to projective measurement, weak measurements only perturb the system
slightly, but in turn can only provide limited information. According to the
theory of open quantum systems, both the evolution of the quantum state and
its decoherence depend on the system-apparatus coupling strength and the basis
in which the measurement system is measured. One of the advantages of such a
measurement is that, given a weak continuous measurement record one can
reconstruct the quantum system state during its evolution. One particularly
interesting approach, which we apply and investigate here, is the “quantum
filter equation” as pioneered in the quantum limit by Belavkin Belavkin1999
and others Carmichael .

The quantum filter equation has been shown to be a powerful method for state
reconstruction, and is fairly robust in terms of the resolution needed in
describing the measurement record. For example, in recent work Ralph  it has
been shown that the continuous “analogue” measurement record can be reduced to
a “one-bit record” and still the filter equation can efficiently produce an
estimate (or purification) of the system state. Similarly it has been shown
that feedback control can be used to enhance the speed of the state estimation
Filter ; Chase ; Zhang , and that it can be further optimized when combined
with a kind of process tomography Deutsch . However, successful application of
the filter equation requires accurate knowledge of the evolution, both
coherent and incoherent, that the monitored system undergoes Hill2011 .

Here we first investigate how this state estimation method can be used by
considering both how confidently Croke2006  an estimated state deduced from
the measurement record reflects the actual state of the system, and how the
effect of backaction changes the state. Here the “actual state of the system”
means the system under the influence of the back-action backaction ; Buscemi ;
Luo  of the measurement apparatus, not the initial prepared state. In other
words, we quantify the confidence of the state-estimation process
independently from the overall fidelity of the measurement. Thus, we focus on
understanding the state reconstructed using the filter equation, without
actively undoing the backaction or employing feedback. We separately define
the overall accuracy of the measurement with another quantity or distance,
which we term the “epitome”. Second, we introduce a new ensemble averaged
version of the filter equation which enables a more efficient numerical method
(in theoretical treatments) with which to consider state estimation via the
filter equation. Third, in the final section we consider the more general
situation of an asymptotic positive operator valued measure (POVM), and show
that in a limiting case the confidence has a lower bound set by the quantum
discord.

Quantities like the confidence have been commonly employed in investigating
information gain with projective and general measurements Banazek ; Fuchs .
The behavior of these quantities in the context of weak continuous measurement
has not been studied in great detail as of yet, though the concept of
information gain and measurement disturbance is well understood Jacobs ;
Jacobs2 . In addition, we point out that our approach here is different from
that used in some other works. For example, here we are concerned primarily
with the optimal measurement of an unknown state by refining our state of
knowledge, whereas in some other approaches the goal is primarily manipulating
(or purifying) one’s state of knowledge of a given quantum system, and the
initial unknown quantum state is unimportant Jacobs2 .

The model we use here is based on the continuous measurement of a double
quantum dot (DQD) charge state using a quantum point contact (QPC) Korotkov99
; DQD ; Petta ; Lambert2010 . In this situation the tunnelling barrier of the
QPC is modulated by the charge in the nearby DQD, and produces a continuous
measurement record which can be used in the filter equation. However, the
approach is quite general, and recent applications have also arisen in
circuit-QED systems Johansson ; Korotcqed .

This article is organized as follows: In Sec. II we provide a general
definition of the confidence, backaction and epitome. In Sec. III we introduce
the model for weak measurement of a DQD, and show the results given by the
quantum filter equation. In Sec. IV we speculate about a possible filter
equation based on ensemble-averaged measurements. Finally, in Sec. V we show
that using the relative entropy as a measure of confidence allows us to define
the lower bound on the confidence as the quantum discord.

Figure 1: (Color online) (a) Schematic of the components of the quantum state
estimation using the filter equation. Here, $\rho_{s}$ is the system initial
state, $\rho_{R}(t)$ is the evolution of the initial state of the system in
the measurement-induced noises, $\rho_{E}(t)$ is the estimation of the system
state, and $\mathcal{D}[\rho_{E},\rho_{R}]$ is the estimation confidence. (b)
shows a diagrammatic representation of the confidence $\mathcal{C}$ and
backaction $\mathcal{B}$ as the distance between the various states. Note that
when we use the relative entropy as a distance measure these quantities are
asymmetric (as represented by the one-way arrows).

##  II Definition of Confidence and Backaction

Consider a quantum system $\mathcal{S}$ which interacts with a measurement
apparatus $\mathcal{A}$. When implementing a general quantum measurement
$\mathcal{M}=\\{\Pi_{i}\\}$, our knowledge of the system state is based only
on the measured apparatus output $y(t)$. We denote our estimate of the system
state as $\rho_{E}(t)$. The initial state of the system is defined as
$\rho_{s}(0)$, and $\rho_{R}(t)$ is defined as the evolution of the initial
state of the system in the measurement induced noise, i.e., what we can think
of as the actual state of the system. We define the ideal state of the system,
i.e., its evolution if it were not affected by the measurement apparatus at
all, as $\rho_{I}(t)$, i.e., the entirely coherent evolution of $\rho_{s}(0)$.

Definition: The quantum state estimation confidence is defined as

$\displaystyle\mathcal{C}\equiv D[\rho_{E},\rho_{R}],$ (1)

where $D[\cdot]$ is some appropriate distance measure. Generally speaking, the
smaller $\mathcal{C}$ is the more confident we are of the estimated state, and
$\mathcal{C}=0$ if and only if $\rho_{E}=\rho_{R}$, which means that the
estimated state is totally confident. Similarly we define the backaction by

$\displaystyle\mathcal{B}\equiv D[\rho_{I},\rho_{R}],$ (2)

so that $\mathcal{B}=0$ implies no backaction. Finally, we define the overall
accuracy of the measurement with the “epitome”,

$\displaystyle\mathcal{E}\equiv D[\rho_{I},\rho_{E}],$ (3)

which naturally completes the triangle in Fig. 1. In the treatment that
follows of course we have full knowledge of all these states at all times, and
can thus in a theoretical sense identify the optimal parameters that minimize
these quantities.

There is some freedom in choosing an appropriate measure for $\mathcal{C}$,
$\mathcal{B}$, and $\mathcal{E}$. Here we explicitly consider both the
fidelity and the relative entropy Modi . The fidelity is commonly used to
measure the effectiveness of the filter equation, and is defined as

$\displaystyle F=1-\mathcal{C}=|(\sigma^{1/2}\rho\sigma^{1/2})^{1/2}|^{2}.$
(4)

We define the confidence as $\mathcal{C}=1-F$ in this case, to match our
definition of a distance measure, so that high fidelity implies
$\mathcal{C}=0$ (and the same with the other measures). However since the
fidelity is a pseudo-distance this lacks some characteristics of a true
measure. In the case of the relative entropy we define,

$\displaystyle\mathcal{C}=S(\rho_{R}||\rho_{E}),\quad\mathcal{B}=S(\rho_{I}||\rho_{R}),\quad\mathcal{E}=S(\rho_{I}||\rho_{E})$
(5)

where

$\displaystyle S(\sigma||\rho)=-\mathrm{Tr}\sigma\ln\rho-S(\sigma).$ (6)

The relative entropy can be seen as a distance measure, though as it is
asymmetric in $\sigma\leftrightarrow\rho$ it is technically not. In fact the
ordering here is important; with the inverse ordering the backaction
$B\rightarrow\infty$ as $\rho_{I}$ is a pure state in this example. Using the
relative entropy allows us to make a more direct connection to a general POVM
description of a weak continuous measurement in the final section of this
work.

Figure 2: (Color online) Schematic of a QPC used for measuring the electron
states yields backaction on the DQD.

##  III Continuous Weak measurement of a Double Quantum Dot

We now present the specific details of the DQD and QPC system. A DQD consists
of a dot $L$, connected to an emitter, and dot $R$, connected to a collector
DQD ; Ouyang10 ; Petersson10 . As is typical, we assume that the DQD is in the
strong Coulomb regime such that only one electron is allowed in the whole DQD.
Here we assume the DQD is prepared in a single electron state, then isolated
from the emitter and collector reservoirs via manipulation of gate voltages.
The two single-dot states are denoted by $|L\rangle$, and $|R\rangle$. The
Hamiltonian of the DQD can be written by

$\displaystyle H_{{\rm DQD}}=\Omega\sigma_{x}/2+\epsilon\sigma_{z},$ (7)

with $\sigma_{x}=|L\rangle\langle R|+|R\rangle\langle L|$,
$\sigma_{z}=|L\rangle\langle L|-|R\rangle\langle R|$, $\epsilon$ is the level
splitting between the two single-dot states and $\Omega$ is the coherent
tunnelling amplitude betweens the two dots. The nearby QPC has the Hamiltonian
$H_{{\rm
QPC}}=\sum_{k}\epsilon_{1k}a_{1k}^{{\dagger}}a_{1k}+\sum_{q}\epsilon_{2q}a_{2q}^{{\dagger}}a_{2q}$,
and the interaction Hamiltonian
$H_{I}=\sum_{k,q}(\chi-\chi_{L}|L\rangle\langle L|-\chi_{R}|R\rangle\langle
R|)(a_{1k}^{{\dagger}}a_{2q}+a_{2q}^{{\dagger}}a_{1k})$, which is modulated by
the electron states of the DQD. Here $\chi$ is the tunneling amplitude of an
isolated QPC, $\chi_{L}(\chi_{R})$ gives the variation in the tunneling
amplitude when the extra electron stays on the left (right) dot, and
$a_{1k}(a_{2q})$ denotes the electron annihilation operator for the source
(drain) of the QPC. Because the height of the tunneling barrier in the QPC
depends on the electron states of the DQD, it is expected that the measured
current of the QPC will vary with the DQD states.

In this simple model, the current shot noise Aguado04  of the QPC acts as a
noise source. In the low-temperature limit with $k_{B}T\ll\hbar\omega$, the
noise spectral density takes the form Aguado00 ; Aguado04 ; Gustavsson07 ;
Petersson10

$\displaystyle J_{I}(\omega)=\frac{4}{R_{K}}D(1-D)\frac{(eV_{\rm
QPC}-\hbar\omega)}{[1-\exp\\{-(eV_{\rm QPC}-\hbar\omega)/k_{B}T\\}]},$ (8)

where $R_{k}=h/e^{2}$ is the von Klitzing constant, $D$ is the transparency of
a single channel in the QPC, and $D$ is a function of $\chi$, $\chi_{L}$ and
$\chi_{R}$ Korotkov01 ; Aguado00 ; Gustavsson07 ; Lesovik . To treat the
measurement signal, or current through the QPC, as a classical variable one
must assume that the QPC evolution is much faster than the DQD, so that only
the zero-frequency component is important; this is effectively a Markovian
limit in terms of treating the QPC backaction on the DQD Korotkov99 ;
backaction .

We also treat the QPC as a perfect detector. In this limit we can define the
measurement strength of the QPC as

$\displaystyle\kappa=\frac{(\Delta I)^{2}}{16J(0)}.$ (9)

Here, $\Delta I=I_{L}-I_{R}$, where
$I_{L}=D_{L}e^{2}V/\pi\hbar$,$I_{R}=D_{R}e^{2}V/\pi\hbar$, $D_{L}=D+\Delta D$,
$D_{R}=D-\Delta D$, and $\Delta D$ is the change in the transmission of the
QPC due to the charge state of the DQD.

Figure 3: (Color online). The top row of figures shows the ensemble-averaged
behavior (over $2,000$ realizations) of (a) one minus the confidence (blue)
and one minus the backaction (red) defined in terms of the fidelity, (b) the
occupation of the left “dot” or state for the system, Eq. (III) (in red) and
the estimator, Eq. (III.1), (in blue), and (c) the confidence (blue) and
backaction (red) defined in terms of the relative entropy. In all figures
$\kappa=0.005\Omega$, and time is evolved for the equivalent of $15$ Rabi
oscillations of the bare quantum dot state. Figure (b) shows more directly how
the estimator quickly oscillates in phase with the system state, but takes
time to evolve to the same population magnitude. The small oscillations seen
in both backaction curves in (a) and (c) are typical for the definition of the
backaction, and simply represent the “nodes” in the oscillation curve of,
e.g., $P_{L}$, where the $\rho_{I}$ and $\rho_{R}$ states coincide. Figures
(d)-(e) show the same quantities as (a)–(c) except just a single realization.
The estimator state still quickly approaches the system state, as is typical
with the filter equation. The backaction changes in magnitude drastically in
both (d) and (f). There is no single-realization dephasing in (e) because we
assume the QPC to be a perfect detector. This is in contrast to recent work
Johansson  on circuit QED where some information can remain hidden due to the
finite lifetime of the measurement cavity. Note, in (a) and (d) we have
plotted $1-C$ and $1-B$, so that the fidelity result can be easily compared to
other works investigating the effectiveness of the filter equation using
fidelity.

The evolution of the real state of the DQD $\rho_{R}$ can be derived using the
Bayesian techniques of Korotkov Korotkov99 ; Korotkov01 ; Korotkov06  to give
the selective stochastic master equation (SME) in Ito form,

$\displaystyle\rm{d}{\rho_{R}}$ $\displaystyle=$
$\displaystyle-\frac{i}{\hbar}\left[H_{\rm
DQD},\rho_{R}\right]dt+\kappa\mathcal{D}[\sigma_{z}]\,\rho_{R}\,dt$
$\displaystyle+$
$\displaystyle\sqrt{2\kappa}\;\mathcal{H}[\sigma_{z}]\;\rho_{R}\;dW_{R},$

where $\kappa$ is defined above,

$\displaystyle\mathcal{H}[\sigma_{z}]\rho_{R}\equiv\sigma_{z}\rho_{R}+\rho_{R}\sigma_{z}^{\dagger}-\mathrm{Tr}(\sigma_{z}\rho_{R}+\rho_{R}\sigma_{z}^{\dagger})\rho_{R},$
(11)

and the real Wiener process $dW_{R}$ satisfies $\mathrm{E}(dW_{R})=0,$
$(dW_{R})^{2}=dt$. In Eq. (III), the super-operator $\mathcal{D}$ is defined
as

$\displaystyle\mathcal{D}[a]\rho=a\rho
a^{{\dagger}}-\frac{1}{2}(a^{{\dagger}}a\rho+\rho a^{{\dagger}}a).$ (12)

Given that, in a general sense, the measurement operator is
$y=\sqrt{2\kappa\hbar}\sigma_{z}$ the measurement record increment at a time
$t$ is,

$\displaystyle\frac{dy(t)}{\sqrt{\hbar}}$ $\displaystyle=$
$\displaystyle\sqrt{8\kappa}\langle\sigma_{z}^{R}(t)\rangle
dt+dW_{R}=\frac{\Delta
I\langle\sigma_{z}^{R}(t)\rangle}{\sqrt{2J(0)}}dt+dW_{R}.$

Here

$\displaystyle\langle\sigma_{z}^{R}(t)\rangle=\mathrm{Tr}[\sigma_{z}\rho_{R}(t)]$
(14)

is the instantaneous expectation value of $\sigma_{z}$ at time $t$ based on
the selectively evolved density matrix $\rho_{R}(t)$.

This equation of motion also relies on the assumption that

$\displaystyle|I_{L}-I_{R}|=|\Delta I|\ll I_{0}=(I_{L}+I_{R})/2,$ (15)

so that many electrons, $N\geq(I_{0}/\Delta I)^{2}\gg 1$, should pass through
the QPC before one can distinguish the quantum state. This is the weakly
responding or linear regime, and the model as we have described it is entirely
equivalent to the formulation used by Korotkov and others Korotkov99 ;
Korotkov01 ; Korotkov06 . Also, note that for consistency with other works on
the filter equation as a state estimation technique Ralph ; Hill2011  we
describe the noise as a Wiener process, so that the width of the Gaussian
distribution Korotkov01 ; Korotkov06  used to describe the weak measurement
with a QPC is absorbed into $\kappa$.

###  III.1 Quantum Filter Equation

To estimate the quantum state $\rho_{E}$ from the measurement output we employ
the quantum filter equation method Shannon ; Mitter ; Ralph ; Filter ; Chase ;
Zhang . The evolution of the estimated state $\rho_{E}$ is described by the
following stochastic master equation, identical to the “system” one, except
the measurement signal from the system evolution, described above, determines
the noise term:

$\displaystyle d\rho_{E}=-\frac{i}{\hbar}\left[H_{\rm
DQD},\rho_{E}\right]dt+\kappa\mathcal{D}[\sigma_{z}]\rho_{E}dt$
$\displaystyle+\sqrt{2\kappa}\mathcal{H}[\sigma_{z}]\rho_{E}\left[\frac{dy(t)}{\sqrt{\hbar}}-\frac{\Delta
I}{\sqrt{2J(0)}}\langle\sigma_{z}^{E}(t)\rangle\right],$ (16)

The last term is analogous to the classical _innovation_ in control theory
Shannon ; Mitter , i.e., the difference between the actually measured current
and the predicted current with the estimated state. The state estimation
process involves setting $\rho_{E}(0)=1/2$, and evolving under the noise
determined by the measurement record from the “experiment”, or in this
theoretical work Eq. (III).

As demonstrated in Ref. [Hill2011, ], even a small error in the Hamiltonian of
the above equation can induce errors in the estimate of the state provided by
the quantum filter equation. We expect the same to be true of the estimates of
the noise spectrum of the QPC. Here our goal is to study the efficiency and
robustness of quantum state estimation via the filter equation, so as in Ralph
et al. [Ralph, ] we choose the same Hamiltonian and QPC properties in both Eq.
(III.1) and Eq. (III). We leave the problem of accurately estimating the
Hamiltonian and the noise properties of the measurement in a dynamic way
Hill2011  for future work.

Figure 4: (Color online) Ensemble-averaged (a) confidence $C$ and (b)
backaction $B$ in terms of the relative entropy between system and estimator
states, as a function of the interaction time $t$ and measurement strength
$\kappa$. Both figures are derived using the ensemble-averaged equation of
motion (IV). (c) and (d) show a comparison between the probability of
occupation of the left dot $P_{L}(t)$ derived from averaging Eq. [III.1] over
many realizations (dashed blue lines) to that given by solving Eq. [IV].

Finally, we combine the time-evolution of these two equations (III) and
(III.1) to calculate the confidence of the estimated state and backaction of
the measurement using both the fidelity and the relative entropy, as described
earlier. In Fig.3, we show numerical results for these quantities. We will
explicitly discuss the epitome in the next section. We set
$\rho_{R}(0)=|L\rangle\langle L|$, $\rho_{E}(0)=1/2$, and evolve using the
standard techniques for $150000$ time steps per Rabi oscillation. The top row
of figures shows the results averaged over $2000$ realizations, while the
bottom row shows just a single realization. In these results we set $\hbar=1$,
$\epsilon=0$ , $\Omega=1$ and $\kappa=0.005\Omega$. Comparing to real
parameters from Gustavsson07 ; Petersson10 , one could choose a strong inter-
dot coupling of the order of $\Omega=32$ $\mu e\mathrm{V}$, giving a timescale
of $130$ ps for the Rabi oscillations we show in Fig. 3. This should be chosen
carefully however, to match the desired properties of the QPC timescales (or
whatever the measurement device happens to be).

When we acquire information from the measurement, it of course induces
significant backaction on the system itself. Figure 3.(a) shows the confidence
and backaction in terms of the fidelity, while Fig. 3.(c) shows the same in
terms of the relative entropy. Both give reasonable descriptions of the
distance between the estimated state and system state, and for weak
measurement strength the confidence saturates before the backaction does.
Obviously then the trade off is to measure on a time scale where both the
confidence is relatively high and the backaction is low.

To gain more insight on what is actually happening during the evolution of the
filter equation, Fig. 3(b) and (e) show the population of the left state of
the dot for both the system and estimator. In the ensemble averaged case Fig.
3(b) backaction from the measurement dephases the state, but the estimator
matches the system state before coherent information is totally lost. In a
single realization, 3(e), the system state does not dephase because the QPC
acts as a perfect detector. Compare this to the case of circuit-QED where the
lifetime of the cavity can induce dephasing on certain timescales Johansson ;
Korotcqed . The off-diagonal elements behave in a similar fashion.

##  IV A filter equation for ensemble expectation values

Solving for the ensemble-averaged results by collating many single
realizations is sometimes an arduous numerical task, though can be useful in
the stochastic Schrödinger form if one is modelling a system with a large
Hilbert space. In practise the system state density matrix ensemble averaged
over all measurement trajectories can be trivially calculated by averaging Eq.
(III), and noting $E[dW_{R}]=0$. This gives the expected Lindblad equation of
motion which induces the behavior we observe in the backaction. How about the
estimator state? Averaging Eq. (III.1) is non-trivial as the individual
trajectories determined by $\langle\sigma_{z}^{R}(t)\rangle$ are not
statistically independent of $\rho_{E}(t)$.

Rather than attempt to do so we simply write down an analogy to the quantum
filter equation which depends on ensemble-averaged quantities rather than
individual trajectories. We now define $\rho_{E}=\mathrm{E}[\rho_{E}]$,
$\langle\sigma_{z}^{R}(t)\rangle=\mathrm{E}[\langle\sigma_{z}^{R}(t)\rangle]$,
and
$\langle\sigma_{z}^{E}(t)\rangle=\mathrm{E}[\langle\sigma_{z}^{E}(t)\rangle]$.
Thus the term $\mathrm{E}[\langle\sigma_{z}^{R}(t)\rangle]$ represents the
expectation value extracted from an ensemble averaged version of Eq. (III),
i.e., the evolution of the real system under the effect of dephasing. By
comparison to the stochastic filter equation we consider the following
nonstochastic filter equation,

$\displaystyle\frac{\rm{d}\rho_{E}}{\rm{d}t}$ $\displaystyle=$
$\displaystyle-\frac{i}{\hbar}\left[H_{\rm
DQD},\rho_{E}\right]+\kappa\mathcal{D}[\sigma_{z}]\rho_{E}$ $\displaystyle+$
$\displaystyle\sqrt{2\kappa}\mathcal{H}[\sigma_{z}]\rho_{E}\left[\frac{\Delta
I}{\sqrt{2J(0)}}\left\\{\langle\sigma_{z}^{R}(t)\rangle-\langle\sigma_{z}^{E}(t)\rangle\right\\}\right].$

Solving this equation directly is computationally trivial. To illustrate this
we plot the confidence as a function of time and measurement strength $\kappa$
in Fig. 4. We can easily see that as $\kappa$ increases the confidence
saturates quickly, but with an associated strong backaction, as expected.
Remarkably, if we inspect the density matrix elements of the estimated state
generated by Eq. (IV) to those generated by the ensemble average over
trajectories of Eq. (III.1) they coincide closely. This is illustrated in Fig.
4(c) and (d). Curiously we are unable to rigorously justify this
correspondence, though one can note that Eq. (IV) represents a valid solution
to Eq. (III.1) for the trajectory determined by $dW_{R}=0$. We also point out
that the fictitious nonlinear force in the second line of Eq. (IV) is not
physical, and the equation may not ensure positivity of the density matrix
$\rho_{E}$ at an arbitrary time. Why this works so well in reproducing results
from the ensemble averaged filter equation, at least for this case of a single
qubit measured in the $\sigma_{z}$ basis, is not clear, and represents a
possible avenue of future work.

Figure 5: (Color online) The epitome, $\mathcal{E}=S(\rho_{I}||\rho_{E})$, as
a function of the interaction time $t$ and measurement strength $\kappa$. This
figure is derived using the ensemble-averaged equation of motion (IV). The
superimposed white squares indicates the line of crossing points between the
confidence and backaction of figure 4, and naturally is close to the minimum
in the epitome.

Finally, in Fig. 5., we use this nonstochastic filter equation to plot the
epitome, $\mathcal{E}=S(\rho_{I}||\rho_{E})$. We see that at some intermediate
time, depending on the measurement strength, the epitome has an optimal minima
which coincides closely to the crossing point of the confidence and
backaction. Thus, as one would expect, with continous weak measurements there
is an optimal time at which our estimated state is closest to the original,
unperturbed, ideal state. In practise this optimization can be also discussed
in terms of goal programming (shown in Appendix A) Goal . Further methods to
improve the estimation, or minimize backaction, could include feedback or
other techniques from quantum control ueda .

##  V POVM and the discord as a bound on confidence

In quantum theory one can describe any measurement scenario as a positive-
operator valued measure (POVM). For example, weak measurement is sometimes
discussed in terms of a POVM on a combined system/measurement-apparatus, where
the measurement apparatus itself is also considered to be a quantum system. To
gain a more fundamental perspective on the confidence and backaction, as we
have defined them here, we consider an alternative measurement scheme of an
asymptotic POVM.

First, we retain our definition of the initial pure system-apparatus state
$\rho_{I}(0)=\rho_{s}(0)\otimes\rho_{A}(0)$. We then assume that a measurement
apparatus $A$ is allowed to interact with the system to produce the typically
entangled and correlated system-apparatus state $\rho_{s,A}$ (we suppress any
time argument here for complete generality). We define the “real” state of the
system then as this combined state $\rho_{s,A}$. Finally, we perform
measurements on $A$ associated with a POVM
$\\{\Pi_{j}^{A,\dagger}\Pi_{j}^{A}\\}$, where
$\sum_{j}\Pi_{j}^{A,\dagger}\Pi_{j}^{A}=1$. Our estimate of the combined
system-apparatus state given by the measurement is
$\rho_{m}=\sum_{j}\Pi_{j}^{A}\rho_{R}\Pi_{j}^{A,\dagger}$. Again we define the
confidence in terms of the relative entropy, so $C=S(\rho_{s,A}||\rho_{m})$.

The relative-entropy-based confidence has an interesting lower bound in the
case when the POVM becomes a projective valued measure. Writing,

$\displaystyle
C=-\mathrm{Tr}\rho_{s,A}\ln\sum_{j}\Pi_{j}^{A}\rho_{s,A}\Pi_{j}^{A,\dagger}-S(\rho_{s,A})$
(18)

we can substitute $\Pi_{j}^{A}=|j\rangle\langle j|$, where $|j\rangle$ is some
distinguishable orthonormal basis describing the measurement apparatus. Then
the confidence becomes,

$\displaystyle C$ $\displaystyle=$
$\displaystyle-\mathrm{Tr}\sum_{j}|j\rangle\langle
j|\rho_{s,A}\ln\sum_{j}|j\rangle\langle j|\rho_{s,A}|j\rangle\langle
j|-S(\rho_{s,A})$ (19) $\displaystyle=$
$\displaystyle-\mathrm{Tr}\sum_{j}|j\rangle\langle
j|\rho_{s,A}|j\rangle\langle j|\ln\sum_{j}|j\rangle\langle
j|\rho_{s,A}|j\rangle\langle j|$
$\displaystyle-S(\rho_{s,A})=S(\rho_{m})-S(\rho_{s,A})$ $\displaystyle\geq$
$\displaystyle\mathrm{min_{|j\rangle}}[S(\rho_{m})]-S(\rho_{s,A})=\mathcal{D},$

where $\mathcal{D}$ is the quantum discord Modi , when they assume
classicality in terms of only one sub-system. In their work the discord has
the meaning of the distance between a given state and the closest (system-
apparatus) separable state. In other words the lower bound on the confidence
is the distance between the real state and the closest separable state, as one
would expect with projective measurements.

In the case of a general POVM, one could argue that the minimization of $C$
over all possible POVMs is equivalent to a generalization of the definition of
Modi’s discord Modi . Finally, we note that there is a correspondence between
the estimator state $\rho_{E}$ we discussed in terms of the filter equation,
and the partial trace $\mathrm{Tr}_{A}(\rho_{m})$ over the state constructed
from asymptotic POVM measurements (and the same for the real state $\rho_{R}$
and $\mathrm{Tr}_{A}(\rho_{s,A})$ evolved in the measurement noise in the
filter equation example).

##  VI Conclusion

In many realistic quantum readout architectures the reliability of the quantum
measurement output is an important issue. In this article we discussed how to
measure the confidence and the backaction of a state reconstructed from
continuous weak quantum measurement. As a typical example, we considered a DQD
measured by a nearby QPC. Based on the theory of open quantum systems and the
quantum filter equation method we briefly discussed the trade-off between
measurement confidence and measurement-induced backaction. We also considered
a possible filter equation for ensemble averaged results. We finished by
discussing the case of a general POVM and how the confidence (when defined as
a relative entropy) has a lower bound related to the quantum discord.

###### Acknowledgements.

We thank Kurt Jacobs for helpful comments. W.C. is supported by the RIKEN FPR
Program. Y.O. is partially supported by the Special Postdoctoral Researchers
Program, RIKEN. F.N. is partially supported by the ARO, JSPS-RFBR Contract No.
12-02-92100, a Grant-in-Aid for Scientific Research (S), MEXT “Kakenhi on
Quantum Cybernetics,” and the JSPS via its FIRST program. J.Q.Y. is partly
supported by the National Basic Research Program of China Grant
No.2009CB929302, NSFC Grant No.91121015, and MOE Grant No.B06011. X.Y.L. is
supported by a JSPS Fellowship.

##  Appendix A Quantitative characterization of confidence and backaction via
goal programming

Figure 6: (Color online) Density profile of the optimization function
$\mathcal{O}$ in the goal programming model for several cases: (a)
$\eta_{1}=\eta_{2}=1$ and $\Delta_{C}=\Delta_{B}=0.1$; (b)
$\eta_{1}=\eta_{2}=1$ and $\Delta_{C}=\Delta_{B}=0.2$; (c)
$\eta_{1}=1,\eta_{2}=0.5$ and $\Delta_{C}=\Delta_{B}=0.1$; (d)
$\eta_{1}=0.5,\eta_{2}=1$ and $\Delta_{C}=\Delta_{B}=0.1$.

The results in Section IV show a clear trade-off relation between the
confidence and the backaction in the parameter space spanned by the
interaction time $t$ and the measurement strength $k$. Let us examine this
relation via a sophisticated method, goal programming Goal . We formulate our
problem setting more specifically; we determine $t$ and $k$ such that
$\mathcal{C}\leq\Delta_{C}$ and $\mathcal{B}\leq\Delta_{B}$ for given (small)
positive constants $\Delta_{C}$ and $\Delta_{B}$. The two parameters
$\Delta_{C}$ and $\Delta_{B}$ are regarded as, respectively, admissible
confidence error and permissible backaction. Thus, we can obtain a good
measurement scenario to increase the confidence (i.e., minimize $\mathcal{C}$)
while reducing the backaction (i.e., minimizing $\mathcal{B}$). Hereafter, the
confidence and the backaction are defined via the relative entropy, as seen in
Eq. (5).

Goal programming Goal  provides an optimization method to deal with two (or
more than two) conflicting objectives and it is widely used in mathematics,
information theory and engineering. Instead of finding solutions which
absolutely minimize or maximize objective functions, the task of goal
programming is to find solutions that, if possible, satisfy a set of goals, or
otherwise violate the goals minimally. This makes the approach more appealing
to practical designers, compared to other optimization methods (e.g., linear
programming models). Let us describe our measurement problem in terms of goal
programming:

$\displaystyle{\rm
Minimize}\quad\mathcal{O}\equiv\eta_{1}\delta_{1}^{+}+\eta_{2}\delta_{2}^{+},$
$\displaystyle{\rm
subject\,\,to}\,\left\\{\begin{array}[]{l}\mathcal{C}(\vec{x})-\delta_{1}^{+}+\delta_{1}^{-}=\Delta_{C}\\\
\mathcal{B}(\vec{x})-\delta_{2}^{+}+\delta_{2}^{-}=\Delta_{B}\\\
\delta_{1}^{\pm},\delta_{2}^{\pm}\geq 0\\\
\vec{x}=(t\Omega/2\pi,k\Omega)\in\mathcal{M}\end{array}\right..$ (24)

The weight factors $\eta_{1}$ and $\eta_{2}$ are given positive number, and
represent the relative priority of objectives. If $\eta_{1}>\eta_{2}$, the
condition for the confidence is prior to the one for the backaction, and vice
versa. The condition for the confidence ($\mathcal{C}\leq\Delta_{C}$) is
reformulated by the relation
$\mathcal{C}+\delta_{1}^{+}-\delta_{1}^{-}=\Delta_{C}$, with the deviations
between the admissible error and the actual value, $\delta_{1}^{+}$ and
$\delta_{1}^{-}$. When $\mathcal{C}>\Delta_{C}$ ($\mathcal{C}\leq\Delta_{C}$),
we set $\delta_{1}^{+}=\mathcal{C}-\Delta_{C}$ and $\delta_{1}^{-}=0$
($\delta_{1}^{+}=0$ and $\delta_{1}^{-}=\Delta_{C}-\mathcal{C}$). Similarly,
we can set $\delta_{2}^{\pm}$ via
$\mathcal{B}-\delta_{2}^{+}+\delta_{2}^{-}=\Delta_{B}$. The set
$\mathcal{M}=\\{t\Omega/2\pi,k\Omega\\}$ is the family of the measurement
scenarios. The smaller $\mathcal{O}$, the better performance of the
measurement scenario. The minimum value of $\mathcal{O}$ ($\mathcal{O}=0$)
corresponds to the best solution.

Figure 6 shows the contour profile of $\mathcal{O}$ for various cases, based
on the ensemble-averaged confidence $\mathcal{C}$ and backaction
$\mathcal{B}$. The optimization function $\mathcal{O}$ is a function of the
measurement scenarios: the interaction time $t\Omega/2\pi$ and measurement
strength $k\Omega$. The other parameters are the same as in Figs. 4(a) and
(b). In both Figs. 6(a) and (b), the confidence and the backaction have equal
importance ($\eta_{1}=\eta_{2}=1$). In Fig.6(a) we examine the case that
$\Delta_{C}=0.1$ and $\Delta_{B}=0.1$. We find that the best solution
($\mathcal{O}=0$) appears in the dark blue area. If we relax the restriction
to $\Delta_{C}=\Delta_{B}=0.2$, we find that more solutions for the
optimization $\mathcal{O}=0$, as seen in Fig. 6(b). We also examine the cases
when the confidence has a different importance, or weight, than the
backaction, as seen in Figs. 6(c) and (d). The solution for the case where the
measurement confidence is more important than the backaction
($\eta_{1}=1,\eta_{2}=0.5$) is given in Fig. 6(c). The solution for the
opposite situation ($\eta_{1}=0.5,\eta_{2}=1$) is given in Fig. 6(d). In the
above cases, we have considered a double-criterion goal programming problem
and we find that it is convenient for discussing the trade-off relation
between measurement confidence and measurement induced backaction.

## References

  * (1) H. M. Wiseman and G. J. Milburn, _Quantum Measurement and Control_ (Cambridge, New York, 2010).
  * (2) R. Alicki and M. Fannes, Rep. Math. Phys. 55, 47 (2005).
  * (3) A. G. Kofman, S. Ashhab, and F. Nori, Phys. Rep. 520.43 (2012).
  * (4) S. Ashhab, J. Q. You, and F. Nori, Phys. Rev. A 79, 032317 (2009).
  * (5) S. Ashhab, J. Q. You, and F. Nori, New J. Phys. 11, 083017 (2009).
  * (6) S. Croke, E. Andersson, S. M. Barnett, C. R. Gilson, and J. Jeffers, Phys. Rev. Lett. 96, 070401 (2006).
  * (7) V. R. Braginsky and F. Y. Khalili, _Quantum measurement_ (Cambridge, New York, 1992).
  * (8) A. N. Korotkov, Phys. Rev. B 60, 5737 (1999).
  * (9) A. N. Korotkov and D. V. Averin, Phys. Rev. B 64, 165310 (2001).
  * (10) A. N. Korotkov and A. N. Jordan, Phys. Rev. Lett. 97, 166805 (2006).
  * (11) I. Buluta, S. Ashhab, F. Nori, Rep. Prog. Phys, 74, 104401 (2011).
  * (12) J. Q. You and F. Nori, Nature 474, 589 (2011).
  * (13) J. Q. You and F. Nori, Phys. Today 58 (11), 42 (2005).
  * (14) A. F. Kockum, L. Tornberg and G. Johansson, Phys. Rev. A. 85, 052318 (2012).
  * (15) V. P. Belavkin, Rep. Math. Phys. 43 (3), A405 (1999).
  * (16) H. Carmichael, _An Open Systems Approach to Quantum Optics_ (Springer, New York, 1991).
  * (17) J. F. Ralph and N. P. Oxtoby, Phys. Rev. Lett. 107, 260503 (2011).
  * (18) A. C. Doherty, S. Habib, K. Jacobs, H. Mabuchi, and S. M. Tan, Phys. Rev. A 62, 012105 (2000).
  * (19) B. A. Chase, A. J. Landahl, and J. M. Geremia, Phys. Rev. A 77, 032304 (2008).
  * (20) J. Zhang, Y-X. Liu, R-B. Wu, C-W. Li, and T-J. Tarn , Phys. Rev. A 82, 022101(2010).
  * (21) A. Silberfarb, P. S. Jessen, and I. H. Deutsch, Phys. Rev. Lett. 95, 030402 (2005).
  * (22) J. F. Ralph, K. Jacobs, and C. D. Hill, Phys. Rev. A 84, 052119 (2011).
  * (23) S. Croke, E. Andersson, S. M. Barnett, C. R. Gilson, and J. Jeffers, Phys. Rev. Lett. 96, 070401(2006).
  * (24) M. Ozawa, Ann. Phys. (Leipzig) 311, 350 (2004).
  * (25) F. Buscemi, M. Hayashi, and M. Horodecki, Phys. Rev. Lett. 100, 210504 (2008).
  * (26) S. Luo and N. Li, Phys. Rev. A 84, 052309 (2011).
  * (27) K. Banaszek, Phys. Rev. Lett. 86, 1366 (2001).
  * (28) C. A. Fuchs and A. Peres, Phys. Rev. A 53, 2038 (2001).
  * (29) A. C. Doherty, K. Jacobs, and G. Jungman, Phys. Rev. A 63, 062306 (2001).
  * (30) K. Jacobs, Phys. Rev. A 67, 030301 (R) (2003).
  * (31) S. A. Gurvitz, Phys. Rev. B 56, 15215 (1997).
  * (32) J. R. Petta _et al._ , Science 309, 2180 (2005).
  * (33) N. Lambert, C. Emary, Y. N. Chen, and F. Nori, Phys. Rev. Lett. 105, 176801 (2010).
  * (34) A. N. Korotkov, arXiv:1111.4016v1, (unpublished).
  * (35) K. Modi, T. Paterek, W. Son, V. Vedral, and M. Williamson, Phys. Rev. Lett. 104, 080501, (2010).
  * (36) S. H. Ouyang, C. H. Lam, and J. Q. You, Phys. Rev. B 81, 075301 (2010).
  * (37) K. D. Petersson, J. R. Petta, H. Lu, and A. C. Gossard, Phys. Rev. Lett. 105, 246804 (2010).
  * (38) R. Aguado and T. Brandes, Phys. Rev. Lett. 92, 206601 (2004).
  * (39) R. Aguado and L. P. Kouwenhoven, Phys. Rev. Lett. 84, 1986 (2000).
  * (40) S. Gustavsson, M. Studer, R. Leturcq, T. Ihn, K. Ensslin, D. C. Driscoll and A. C. Gossard, Phys. Rev. Lett. 99, 206804 (2007).
  * (41) G. B. Lesovik, JETP Lett. 49, 592 (1989).
  * (42) H. W. Bode, C. E. Shannon, Proc. IRE 38(4), 417-425 (1950).
  * (43) S. K. Mitter, _Nonlinear filtering of diffusion processes_ , (Springer, New York, 1982).
  * (44) M. J. Schniederjans, _Goal Programming: Methodology and Applications_ (Kluwer, Boston, 1995).
  * (45) Y. Watanabe, T. Sagawa, M. Ueda, Phys. Rev. Lett. 104, 020401, (2010).


# Multivariate permutation entropy,
a Cartesian graph product approach ††thanks: J.S. Fabila-Carrasco and J.
Escudero were supported by the Leverhulme Trust via a Research Project Grant
(RPG-2020-158).

John Stewart Fabila-Carrasco  School of Engineering
Institute for Digital Communications
University of Edinburgh
Edinburgh, EH9 3FB, UK
<EMAIL_ADDRESS>Chao Tan  School of Electrical and
Information Engineering
Tianjin University
Tianjin, 300072, China
<EMAIL_ADDRESS>Javier Escudero  School of Engineering
Institute for Digital Communications
University of Edinburgh
Edinburgh, EH9 3FB, UK
<EMAIL_ADDRESS>
###### Abstract

Entropy metrics are nonlinear measures to quantify the complexity of time
series. Among them, permutation entropy is a common metric due to its
robustness and fast computation. Multivariate entropy metrics techniques are
needed to analyse data consisting of more than one time series. To this end,
we present a multivariate permutation entropy, $\operatorname{MPE_{G}}$, using
a graph-based approach.

Given a multivariate signal, the algorithm $\operatorname{MPE_{G}}$ involves
two main steps: 1) we construct an underlying graph ${G}$ as the Cartesian
product of two graphs $G_{1}$ and $G_{2}$, where $G_{1}$ preserves temporal
information of each times series together with $G_{2}$ that models the
relations between different channels, and 2) we consider the multivariate
signal as samples defined on the regular graph ${G}$ and apply the recently
introduced permutation entropy for graphs.

Our graph-based approach gives the flexibility to consider diverse types of
cross channel relationships and signals, and it overcomes with the limitations
of current multivariate permutation entropy.

###### Index Terms:

permutation entropy, graph signals, entropy metrics, complexity, multivariate
time series.

##  I Introduction

Entropy measurements are a common tool used in the analysis of time series to
describe the probability distribution of the states of a system. Based on this
concept, the seminal paper [1] introduced the so-called permutation entropy
($\operatorname{PE}$) as a measure to quantify complexity in time series, a
fundamental challenge in data analysis. This entropy involves calculating
permutation patterns, i.e., permutations defined by comparing neighbouring
values of the time series. $\operatorname{PE}$ has been applied in a wide
range of fields: physical systems [2], economics [3], and biomedicine [4, 5],
among many other applications.

There have been studies on the properties of permutation entropy, including
extensions to higher regular domains [6] and irregular domains or graphs [7].
Some modifications of $\operatorname{PE}$ consider nonlinear mappings to deal
with the differences between the amplitude values [8, 9], or weights in
permutation patterns [10]. Previous research also has extended
$\operatorname{PE}$ to different scales [11, 12], or studied its dependencies
with respect to random signals [13], autoregressive processes [14] or high-
order autoregressive processes [15].

Most biomedical and physical systems are multivariate. Therefore, univariate
entropy metrics have been generalised to a multivariate setting, including:
multivariate sample entropy [16], multivariate dispersion entropy [17], among
others. A multivariate multiscale permutation entropy ($\operatorname{MMPE}$)
to analyse physiological signals is proposed in [18]. However, such algorithm
extracts the permutation patterns from each channel separately regardless of
their cross-channel information. $\operatorname{MMPE}$ treated multichannel
signals as a unique block and without interactions between the channels. Thus,
it works appropriately when the components of a multivariate signal are
statistically independent but does not consider the spatial domain of time
series.

#### Contributions

We introduce a multivariate permutation entropy based on the Cartesian product
of graphs. Such approach would enable us to, first, overcome the limitation of
current multivariate permutation entropy and, second, give flexibility to
consider diverse types of cross channel relationships and signals.

#### Structure of the manuscript

The outline of the paper is as follows: Section II introduces the Cartesian
product of graphs and the permutation entropy: univariate, multivariate and
for graph signals. Section III presents the graph associated with a
multivariate signal, and Section IV presents the multivariate permutation
entropy. Section V shows how $\operatorname{MPE_{G}}$ applies to synthetic
signals. The conclusions and future lines of research are presented in Section
VI.

##  II Background: graph and permutation entropies

In this section, we introduce general background information, including the
definition of a graph and the Cartesian product (Section II-A), the original
permutation entropy (Section II-B) and the recently introduced permutation
entropy for graph signals (Section II-C).

###  II-A Graphs, Cartesian product and graph signal

An _undirected graph_ (or simply graph) $G$ is defined as the pair
$G=(\mathcal{V},\mathcal{E})$ which consists of a finite set of vertices or
nodes $\mathcal{V}=\\{1,2,3,\dots,n\\}$, an edge set
$\mathcal{E}\subset\\{(i,j):i,j\in\mathcal{V}\\}$. The adjacency matrix
$\mathbf{A}$ is the corresponding $N\times N$ symmetric matrix on edges with
entries $1=\mathbf{A}_{ij}=\mathbf{A}_{ji}$ if $(i,j)\in\mathcal{E}$ and $0$
otherwise.

A _directed graph_ or _digraph_ is a graph where each edge has an orientation
or direction.

The _Cartesian product_ of two graphs ${G}=(\mathcal{V},\mathcal{E})$ and
${G}^{\prime}=(\mathcal{V}^{\prime},\mathcal{E}^{\prime})$, denoted
${G}\square{G}^{\prime}$, is the graph defined by:

  1. 1.

the vertex set is given by:

$\mathcal{V}({G}\square{G}^{\prime})=\mathcal{V}\times\mathcal{V}^{\prime}=\\{\,(v,v^{\prime})\,|\,v\in\mathcal{V}\text{
and }v^{\prime}\in\mathcal{V}^{\prime}\,\\}\>;$

  2. 2.

two vertices $(v,v^{\prime})$ and $(u,u^{\prime})$ are adjacent in
${G}\square{G}^{\prime}$ if and only if either

     * •

$v=u$ and $v^{\prime}$ is adjacent to $u^{\prime}$ in $\mathcal{V}^{\prime}$,
or

     * •

$v^{\prime}=u^{\prime}$ and $v$ is adjacent to $u$ in $\mathcal{V}$.

A _graph signal_ is a real function defined on the vertices, i.e.,
$\mathbf{X}\colon\mathcal{V}\longrightarrow\mathbb{R}$. The graph signal
$\mathbf{X}$ can be represented as an $n$-dimensional column vector.

###  II-B Original permutation entropy: univariate and multivariate

#### Univariate permutation entropy

For a time series $\textbf{X}=\left\\{x_{i}\right\\}_{i=1}^{n}$, the algorithm
to compute $\operatorname{PE}$ is the following [1]:

$1)$ For $2\leq m\in\mathbb{N}$ the _embedding dimension_ and $L\in\mathbb{N}$
the _delay time_ , the _embedding vector_
$\textbf{x}_{i}^{m}(L)\in\mathbb{R}^{m}$ is given by
$\textbf{x}_{i}^{m}(L)=\left(x_{i+jL}\right)_{j=0}^{m-1}=\left(x_{i},x_{i+L},\dots,x_{i+(m-1)L}\right)$
for all $1\leq i\leq n-(m-1)L$.

$2)$ The embedding vector
$\textbf{x}_{i}^{m}(L)=\left(x_{i},x_{i+L},\dots,x_{i+(m-1)L}\right)$ is
arranged in the increasing order vector: $\left(x_{i+(k_{1}-1)L}\leq
x_{i+(k_{2}-1)L}\leq\dots\leq x_{i+(k_{m}-1)L}\right)$. We use the convention
in [4]. In the case of equal values, the order is given by the corresponding
$k^{\prime}s$. Therefore, any embedding vector $\textbf{x}_{i}^{m}(L)$ is
uniquely mapped onto the vector $(k_{1},k_{2},\dots,k_{m})\in\mathbb{N}^{m}$.

$3)$ The relative frequency for the distinct permutation
$\pi_{1},\pi_{2},\dots,\pi_{k}$, where $k=m!$, is denoted by
$p(\pi_{1}),p(\pi_{2}),\dots,p(\pi_{k})$. The permutation entropy
$\operatorname{PE}$ for the time series X is computed as the normalised
Shannon entropy for the $k$ distinct permutations as follows

$\operatorname{PE}=-\dfrac{1}{\ln(m!)}\sum_{i=1}^{m!}p(\pi_{i})\ln
p(\pi_{i})\;.$

#### Multivariate permutation entropy

$\operatorname{MMPE}$ is proposed in [18]. Let $\mathbf{U}$ be a multivariate
signal, $\operatorname{MMPE}$ applies steps $1)$ and $2)$ from the original
$\operatorname{PE}$ for each channel. The difference is step $3)$, where the
probability distribution aggregates the frequency of patterns from all
channels in the multivariate signal, but it does not account for inter-channel
relationships, i.e.

$3)$ The relative frequencies are denoted by $\\{\pi_{i,j}\\}$, then the
marginal relative frequencies describing the distributions of the patterns is
defined by: $P_{j}=\sum\pi_{i,j}$ for $j=1,2,\dots,m!$. The multivariate
$\operatorname{MMPE}$ is computed as the normalised Shannon entropy for the
marginal relative frequencies:

$\operatorname{MMPE}=-\dfrac{1}{\ln(m!)}\sum_{j=1}^{m!}P_{j}\ln P_{j}\;.$

###  II-C Permutation entropy for graph signals

Let $G=(\mathcal{V},\mathcal{E})$ be a graph, $\mathbf{A}$ its adjacency
matrix and $\textbf{X}=\left\\{x_{i}\right\\}_{i=1}^{n}$ be a signal on the
graph. The permutation entropy for the graph signals $\operatorname{PE_{G}}$
is defined in [7] as follows:

$1)$ For $2\leq m\in\mathbb{N}$ the _embedding dimension_ , $L\in\mathbb{N}$
the _delay time_ and for all $i=1,2,\dots,n$ we define
$y_{i}^{kL}=\frac{1}{\lvert\mathcal{N}_{kL}(i)\rvert}\sum_{j\in\mathcal{N}_{kL}(i)}x_{j}=\frac{1}{\lvert\mathcal{N}_{kL}(i)\rvert}(\mathbf{A}^{kL}\mathbf{X})_{i}\;$,
where $\mathcal{N}_{k}(i)=\\{\,j\in\mathcal{V}\,|\,\scalebox{0.9}[1.0]{it
exists a walk on $k$ edges joining $i$ and $j$}\,\\}\;$. Hence, we construct
the embedding vector $\textbf{y}_{i}^{m,L}\in\mathbb{R}^{m}$ given by
$\textbf{y}_{i}^{m,L}=\left(y_{i}^{kL}\right)_{k=0}^{m-1}=\left(y_{i}^{0},y_{i}^{L},\dots
y_{i}^{(m-1)L}\right)\;$.

$2)$ The embedding vector $\textbf{y}_{i}^{m,L}$ is arranged in increasing
order.

$3)$ The relative frequency for the distinct permutation
$\pi_{1},\pi_{2},\dots,\pi_{k}$, where $k=m!$, is denoted by
$p(\pi_{1}),p(\pi_{2}),\dots,p(\pi_{k})$. The permutation entropy
$\operatorname{PE_{G}}$ for the graph signal X is computed as the normalised
Shannon entropy

$\operatorname{PE_{G}}=-\dfrac{1}{\ln(m!)}\sum_{i=1}^{m!}p(\pi_{i})\ln
p(\pi_{i})\;.$

For time series, $\operatorname{PE_{G}}$ reduces to $\operatorname{PE}$. In
particular, if X is a time series and $G$ the directed path on $n$ vertices,
then for all $m$ and $L$, the equality holds:
$\operatorname{PE}(m,L)=\operatorname{PE}_{G}(m,L)$ (see [7, Prop. 3]).

##  III Construction of the graph

In this section, we will associate a graph for each multivariate signal. Let
$\mathbf{U}$ be a multivariate signal; we will construct a 2D graph (using the
Cartesian product). One dimension will preserve the temporal information, and
another will preserve the cross-channel information.

### Dimension 1. Temporal information

We associate the directed path with a time series, where a vertex represents
each sample time. A _directed path_ on $n$ vertices is a directed graph that
joins a sequence of different vertices with all the edges in the same
direction and is denoted by $\overrightarrow{P_{n}}$, i.e. its vertices are
$\\{1,2,\dots,n\\}$ and its arcs $(i,i+1)$ for all $1\leq i\leq k-1$. An
example is depicted in Fig. 1(a).

### Dimension 2. Relationships between channels

$\mathbf{U}=\\{U_{s}\\}^{s=1,2,\dots,p}$ is a multivariate signal consisting
on a set of $p$ time series (or channels). Let $I_{p}$ be the graph with $p$
vertices representing the interaction between different channels, i.e. $U_{i}$
and $U_{j}$ are adjacent in the graph $I_{p}$ if and only if they interact.

If we do not have any a priori information about the interactions between
channels, by default, we will consider equal interactions between all
channels. Complete graphs represent such relations, i.e., we will set
$I_{p}=K_{p}$ as the complete graph with $p$ vertices, see an example in Fig.
1(b).

######  Definition 1.

Let $\mathbf{U}=\\{u_{t,s}\\}_{t=1,2,\dots,n}^{s=1,2,\dots,p}$ be a
multivariate time series with $p-$channels of length $n$ and with $I_{p}$ the
graph of interactions between channels. We define $\mathcal{G}_{\mathbf{U}}$
as the graph associated with $\mathbf{U}$ and given by:

$\mathcal{G}_{\mathbf{U}}:=\overrightarrow{P_{n}}\square I_{p}\;.$

We set $I_{p}$ as the complete graph in case of not having additional
information about the interaction of the channels.

Figure 1: (a) Directed path with seven vertices, denoted by
$\overrightarrow{P_{7}}$. (b) Interactions between the four channels are
encoded with the complete graph on four vertices, denoted by $K_{4}$. (c) The
Cartesian product $\overrightarrow{P_{7}}\square K_{4}$.

By construction, the vertices in $\mathcal{G}_{\mathbf{U}}$ and the sample
points in the multivariate signal $\mathbf{U}$ are indexed by the same set,
i.e.,

${\mathcal{V}(\mathcal{G}_{\mathbf{U}})=\mathcal{V}(\overrightarrow{P_{n}}\square
I_{p})=\mathcal{V}(\overrightarrow{P_{n}})\times\mathcal{V}(I_{p})=\mathbb{Z}_{n}\times\mathbb{Z}_{p}}\;;$

hence, we can consider $\mathbf{U}$ as a graph signal defined on the regular
domain $\mathcal{G}_{\mathbf{U}}$.

The constructed $\mathcal{G}_{\mathbf{U}}$ is a 2D domain, time defines one
dimension, and cross-channel dependencies define another one. Hence, the graph
$\mathcal{G}_{\mathbf{U}}$ preserve the temporal/dependency structure of the
multivariate signal $\mathbf{U}$. This is not the case for the currently
available implementation of multivariate $\operatorname{PE}$.

######  Example 1.

Consider a multivariate time series with four channels and seven sample
points, i.e., $\mathbf{U}=\\{u_{t,s}\\}_{t=1,2,\dots,7}^{s=1,2,3,4}$, and we
do not have any additional information between channels. By default, we will
assume all channels interact with each other. Fig. 1(c) shows the graph
$\mathcal{G}_{\mathbf{U}}$ constructed in Def. 1.

An example of a multivariate signal where not all channels interact with each
others is shown in Fig. 2(a). The adjacency matrix associated with the graph
$\mathcal{G}_{\mathbf{U}}$ is shown in Fig. 2(b). The graph
$\mathcal{G}_{\mathbf{U}}$ constructed according to Def.1 is depicted in Fig.
2(c).

Figure 2: (a) An undirected graph ${G}$ represents the relationships between
channels. (b) The adjacency matrix corresponding to the graph ${G}$. (c) The
Cartesian product $\overrightarrow{P_{7}}\square G$.

#### More complex relations between channels

Observe that the only imposed condition on the graph $I_{p}$ is the number of
vertices, i.e., $I_{p}$ has as many vertices as the number of channels. Non-
complete graphs can model other dependencies between channels.

We will use undirected edges for bidirectional relationships between channels.
We also can use directed edges for unidirectional interaction and include
weighted edges for heterogeneous relations. Hence, in general, $I_{p}$ would
be a weighted (directed or undirected) graph.

#### Regular structure but not periodic structure

Covering graphs or periodic graphs are used as models of chemical compounds,
like graphene nanoribbons [19]. The graph $\mathcal{G}_{\mathbf{U}}$ is not a
covering graph, but it can be considered as a geometrical perturbation (see,
e.g., [20]) of a periodic graph. Hence, some properties of the periodic graphs
can be preserved in $\mathcal{G}_{\mathbf{U}}$, including spectral properties
[21]. Such properties are important in graph signal processing [22, 23] and
combinatorics [24]. This could be useful to formulate more general improved
multivariate signals entropies.

##  IV Multivariate permutation entropy ($\operatorname{MPE_{G}}$)

In this section, we define the multivariate permutation entropy
$\operatorname{MPE_{G}}$. We use the permutation entropy for graph signals
$\operatorname{PE_{G}}$ (Section II-C) and the graph construction described in
Def. 1.

######  Definition 2.

Multivariate permutation entropy ($\operatorname{MPE_{G}}$)

Let $\mathbf{U}=\\{u_{t,s}\\}_{t=1,2,\dots,n}^{s=1,2,\dots,p}$ be a
multivariate time series with interaction graph $I_{p}$ between channels

  1. 1.

Graph construction. Construct the graph $\mathcal{G}_{\mathbf{U}}$ described
in Def. 1, i.e.,

$\mathcal{G}_{\mathbf{U}}:=\overrightarrow{P_{n}}\square I_{p}\;.$

  2. 2.

Graph signal. Consider $\mathbf{U}$ as a signal defined on the graph
$\mathcal{G}_{\mathbf{U}}$, i.e.,

$\mathbf{U}\colon\mathcal{V}(\mathcal{G}_{\mathbf{U}})\longrightarrow\mathbb{R}\;.$

  3. 3.

$\mathbf{PE}$ for graph signals. The _multivariate permutation entropy (
$\operatorname{MPE_{G}}$)_ is defined as the permutation entropy for the graph
signal $\operatorname{PE_{G}}$ (see Section II-C) for the signal $\mathbf{U}$
and the graph $\mathcal{G}_{\mathbf{U}}$, i.e.,

$\operatorname{MPE_{G}}=\operatorname{PE_{G}}(\mathbf{U})\;.$

Proposition 1 proves some important relations between $\operatorname{MPE_{G}}$
and $\operatorname{PE}$ metrics presented in the literature (see Table I).

######  Proposition 1.

Let $\mathbf{U}=\\{u_{t,s}\\}_{t=1,2,\dots,n}^{s=1,2,\dots,p}$ be a
multivariate time series with interaction graph $I_{p}$ between channels:

  1. 1.

If $p=1$, then $\operatorname{MPE_{G}}(\mathbf{U})=\operatorname{PE}$.

  2. 2.

If $s=1$, then $\operatorname{MPE_{G}}(\mathbf{U})=\operatorname{PE_{G}}$.

  3. 3.

If $I_{p}$ is the graph defined by $p$ isolated vertices, then

$\operatorname{MPE_{G}}(\mathbf{U})=\operatorname{MMPE}\>.$

  4. 4.

If $I_{p}$ is a directed path on $p$ vertices, then

$\operatorname{MPE_{G}}(\mathbf{U})\approx\operatorname{PE}_{2D}\>.$

  5. 5.

If $\mathbf{U}=\\{U_{t}\\}_{t=1,2,\dots,n}^{s=1,2,\dots,p}$, then
$\operatorname{MPE_{G}}(\mathbf{U})\approx\operatorname{PE}(U_{t})$.

###### Proof.

$1)-3)$ are easy properties and follow from the entropy definitions (Table I).
$4)$ If $I_{p}$ is a directed path on $p$ vertices, then
$\mathcal{G}_{\mathbf{U}}$ is a directed grid graph with $n\times p$ vertices,
i.e.,
$\mathcal{G}_{\mathbf{U}}=\overrightarrow{P_{n}}\square\overrightarrow{P_{p}}$.
Hence, by definition of $\operatorname{MPE_{G}}$, the algorithm reduces to
apply $\operatorname{PE_{G}}$ to the signal $\mathbf{U}$ defined on the grid
$\mathcal{G}_{\mathbf{U}}$, and the performance of $\operatorname{PE_{G}}$ and
$\operatorname{PE}_{2D}$ are similar (see [7, Sec. IV]).
$5)$ We prove the case $m=2$ and $L=1$; the other cases are analogous. By
definition of $\operatorname{PE_{G}}$, each vertex belongs to the permutation
pattern $\pi_{1}$ or $\pi_{2}$. For every vertex
$v_{i}\in\mathcal{V}(\overrightarrow{P_{n}})$, where $v_{i}$ is not the last
vertex of the path, it is easy to show that the set of vertices
$\\{\,(v_{i},v_{j})\in\mathcal{V}(\mathcal{G}_{\mathbf{U}})\,|\,j=1,2,\dots,p\,\\}$
also belong to $\pi_{1}$ (similarly with $\pi_{2}$). Hence, the relative
frequencies are preserved in $\mathcal{G}_{\mathbf{U}}$ and
$\overrightarrow{P_{n}}$ (except for at most $s$ vertices corresponding to the
last vertex of each path). Then, for a large $n$, the values of its
corresponding Shannon entropies are close enough. ∎

TABLE I: Summary of some permutation entropy metrics. Entropy metric | Properties/Limitations
---|---
$\mathbf{PE}$ Permutation Entropy Reference: [1] -[5] |  Analyse univariate time series Simple and computationally fast Multiscale extension
$\mathbf{PE_{2D}}$: Regularly sampled 2D data $\mathbf{PE}$ Reference: [6] |  Analyse bidimensional data Multiscale extension Valuable for texture analysis
$\mathbf{PE_{G}}$ $\mathbf{PE}$ for Graph Signals Reference: [7] |  Analyse graph signals (including: time series and image) No multiscale extension yet
$\mathbf{MPE}$ Multivariate Multiscale $\mathbf{PE}$ Reference: [18] |  Analyse multivariate data but as a unique block (no interactions) Multiscale included
$\mathbf{MPE_{G}}$ Multivariate $\mathbf{PE}$, a graph product approach Reference: Definition 2 |  Analyse multivariate data including cross channel relationships Use $\mathbf{PE_{G}}$ for a graph Cartesian product No multiscale extension yet

##  V Experiments

In this section, we apply the algorithm to a set of multivariate synthetic
signals used in the study of dynamical systems.

###  V-A The Hénon map

In discrete-time dynamical systems, one of the most studied is the Hénon map
introduced in [25]. Using $\operatorname{MPE_{G}}$ we can detect dynamical
changes in the two-dimensional system defined by the equations:

$\displaystyle x_{n+1}$ $\displaystyle=1-ax_{n}^{2}+y_{n}$ $\displaystyle
y_{n+1}$ $\displaystyle=bx_{n}\;.$

The map depends on two parameters: $a$ and $b$. For the values $a=1.4$ and
$b=0.3$ indicate the existence of a strange attractor; hence the map is
chaotic. With $b=0.3$ and for other values of the parameter $a$, Fig. 4 shows
that the map may be periodic, chaotic or intermittent.

Figure 3: Orbit diagram for the Hénon map with $b=0.3$.

We will analyse the parameter $a$, we consider $a\in[1,1.4]$ with increments
in steps of $.0001$. For each iteration, we define the multivariate signal
$\mathbf{U}=\\{x_{t},y_{t}\\}_{t=1,2,\dots,n}$. The initial condition
considered are $x_{1}=0.5$, $y_{1}=0.1$ with $n=100$ (similar results are
obtained for other values). We apply $\operatorname{MPE_{G}}(\mathbf{U})$ for
detecting the dynamic between the two signals. We also apply
$\operatorname{PE}$ for each univariate signal and $\operatorname{MMPE}$ to
the multivariate signal, results obtained by $\operatorname{PE}$ and
$\operatorname{MMPE}$ are similar. Fig. 4 shows the entropy values for
$\operatorname{MPE_{G}}(\mathbf{U})$ and $\operatorname{PE}$ for $m=3$ and
$L=1$.

Figure 4: Entropy values computed for $m=3$ and $L=1$.

Our algorithm is able to detect chaotic behaviour and windows of stability.
Moreover, the wider gap between values of $\operatorname{MPE_{G}}$ indicates a
larger sensitivity.

###  V-B Lorenz system

The Lorenz system is an example of a system of ordinary differential
equations. This system has important applications in mechanics, biology, and
circuit theory [26]; and is given by three simultaneous equations:

$\displaystyle{x^{\prime}}$ $\displaystyle=\sigma(y-x),$
$\displaystyle{y^{\prime}}$ $\displaystyle=x(\rho-z)-y,$
$\displaystyle{z^{\prime}}$ $\displaystyle=xy-\beta z.$

Lorenz used the values $\sigma=10$ and $\beta=8/3$. It is well-known that for
$\rho=28$, the system showed chaotic behaviour. For $\rho<1$, the origin is a
global attractor, i.e., all orbits converge to a unique equilibrium point
[26]. $\operatorname{MPE_{G}}$ algorithm captures this fact. Table II shows
the entropy values computed $\operatorname{MPE_{G}}$ for $L=1$ and
$m=3,4,5,6,7$. Entropy values are larger when $\rho>1$, reflecting more
complexity, while for $\rho<1$, the system tends to the equilibrium; hence,
lower values are obtained.

TABLE II: Entropy values for the Lorenz system. | $m=3$ | $m=4$ | $m=5$ | $m=6$ | $m=7$
---|---|---|---|---|---
$\rho=0.8$ | 0.4524 | 0.2860 | 0.1981 | 0.1477 | 0.1166
$\rho=0.9$ | 0.4538 | 0.2878 | 0.1986 | 0.1489 | 0.1169
$\rho=1.2$ | 0.7258 | 0.6673 | 0.5564 | 0.4478 | 0.3787
$\rho=1.3$ | 0.7226 | 0.6872 | 0.5905 | 0.4815 | 0.4136

##  VI Conclusions and future work

We introduced a multivariate permutation entropy to quantify the complexity of
multivariate time series. The algorithm proposed use the Cartesian product of
graphs and the recently introduced permutation entropy for graph signals [7].
Our graph-based approach considers diverse type of cross channel relationships
and overcomes with the limitations of current multivariate permutation
entropy.

Future lines of research-based on the present work are:

####  VI-1 Multivariate dispersion entropy

Using a similar graph-technique presented in this paper, some univariate
metrics can be generalised to multivariate metrics, including Dispersion
Entropy.

####  VI-2 Multiscale permutation entropy for graph signals

$\operatorname{MPE_{G}}$ requires $\operatorname{PE_{G}}$ in its computation.
Multiscale entropy for time series involves downsampling or a coarse-graining
process. Such process is unclear for signals defined in graphs; hence a
multiscale $\operatorname{PE_{G}}$ is still an open issue.

####  VI-3 Interaction between channels changing with time

The graph $\mathcal{G}_{\mathbf{U}}$ (Def. 1) uses the Cartesian product.
Implicitly, we assumed the relationships between channels are preserved along
time. We will explore constructions involving changes in the interactions
between channels respect to time.

####  VI-4 Irregular domains and real-world data

The presented $\operatorname{MPE_{G}}$ deals with 2D constructions and
synthetic signals. We will explore similar entropies techniques for irregular
domains and apply to real-world data, including biomedical signals and phase-
flow patterns.

## References

  * [1] C. Bandt, and B. Pompe, “Permutation Entropy: A Natural Complexity Measure for Time Series”, _Physical Review Letters_ , 2002, 88(17), pp. 174102.
  * [2] R. Yan, Y. Liu, and R.X. Gao, “Permutation entropy: A nonlinear statistical measure for status characterization of rotary machines”, _Mechanical Systems and Signal Processing_ , 2012, 29, pp.474-484.
  * [3] L. Zunino, M. Zanin, B.M. Tabak, D.G. Pérez, O.A. Rosso, “Forbidden patterns, permutation entropy and stock market inefficiency”, _Physica A: Statistical Mechanics and its Applications_ , 2009, 388(14), pp.2854-2864.
  * [4] Y. Cao, W.W. Tung, J.B. Gao, V.A. Protopopescu, and L.M. Hively, “Detecting dynamical changes in time series using the permutation entropy”, _Physical review E_ , 2004, 70(4), p.046217.
  * [5] E. Olofsen, J.W. Sleigh, A. Dahan, “Permutation entropy of the electroencephalogram: a measure of anaesthetic drug effect”, _BJA: British Journal of Anaesthesia_ , 2008, 101(6), pp. 810–821.
  * [6] C. Morel, and A. Humeau-Heurtier, “Multiscale permutation entropy for two-dimensional patterns”, _Pattern Recognition Letters_ , 2021, 150, pp. 139-146.
  * [7] J.S. Fabila-Carrasco, C. Tan, and J. Escudero,“Permutation Entropy for Graph Signal”, _arXiv:2110.00628_ , 2021.
  * [8] H. Azami, and J. Escudero, “Amplitude- and fluctuation-based dispersion entropy”, _Entropy_ , 2018, 20(3), pp.1–21.
  * [9] M. Rostaghi, and H. Azami, “Dispersion entropy: A measure for time-series analysis”, _IEEE Signal Processing Letters_ , 2016, 23(5), pp. 610-614.
  * [10] I. Mitiche, G. Morison, A. Nesbitt, and B.G. Stewart, “ Classification of partial discharge EMI conditions using permutation entropy-based features”. In _25th European Signal Processing Conference (EUSIPCO)_ , August 2017, pp. 1375–1379.
  * [11] M. Costa, A.L. Goldberger, and C.K. Peng, “Multiscale entropy analysis of complex physiologic time series”, _Physical review letters_ , 2002. 89(6), p. 068102.
  * [12] H. Azami, and J. Escudero, “Improved multiscale permutation entropy for biomedical signal analysis: Interpretation and application to electroencephalogram recordings”, _Biomedical Signal Processing and Control_ , 2016, 23, pp.28-41.
  * [13] A. Dávalos, M. Jabloun, P. Ravier, and O. Buttelli, “Theoretical study of multiscale permutation entropy on finite-length fractional gaussian noise.” In _26th European Signal Processing Conference (EUSIPCO)_ , September 2018, pp. 1087–1091.
  * [14] A. Dávalos, M. Jabloun, P. Ravier, and O. Buttelli, “Multiscale permutation entropy: Statistical characterization on autoregressive and moving average processes”, In _27th European Signal Processing Conference (EUSIPCO)_ , September 2019, pp. 1–5.
  * [15] A. Dávalos, M. Jabloun, P. Ravier, and O. Buttelli, “On the Deterministic Estimaton of Multiscale Permutation Entropy of High-Order Autoregressive-Moving-Average Processes as a Function of ARMA Parameters”. In _28th European Signal Processing Conference (EUSIPCO)_ , January 2021, pp. 2225-2229.
  * [16] M.U. Ahmed and D.P. Mandic,“Multivariate multiscale entropy: A tool for complexity analysis of multichannel data,” _Physical Review E_ , 2011, 84, p. 61918.
  * [17] H. Azami, A. Fernández, and J. Escudero,“Multivariate Multiscale Dispersion Entropy of Biomedical Times Series,” _Entropy_ , 2019, 21(9), Sep. 2019, p. 913.
  * [18] F.C. Morabito, D. Labate, F. La Foresta, A. Bramanti, G. Morabito, and I. Palamara, “Multivariate Multi-Scale Permutation Entropy for Complexity Analysis of Alzheimer’s Disease EEG,” _Entropy_ , 2012, 14(7), pp. 1186-1202.
  * [19] J.S. Fabila-Carrasco, and F. Lledó, “Covering graphs, magnetic spectral gaps and applications to polymers and nanoribbons”, _Symmetry_ , 2019, 11(9), pp. 1–21.
  * [20] J.S. Fabila-Carrasco, F. Lledó, and O. Post, “Spectral preorder and perturbations of discrete weighted graphs,” _Mathematische Annalen_ , 2020, pp. 1-–49.
  * [21] J.S. Fabila-Carrasco, F. Lledó, and O. Post, “Spectral gaps and discrete magnetic Laplacians,” _Linear Algebra and its Applications_ , 2018, 547, pp. 183–216.
  * [22] A. Ortega, P. Frossard, J. Kovačević, J.M. Moura, and P. Vandergheynst, “Graph signal processing: Overview, challenges, and applications”, _Proceedings of the IEEE_ , 2018, 106(5), pp. 808-828.
  * [23] L. Stankovic, D. Mandic, M. Dakovic, M. Brajovic, B. Scalzo, and T. Constantinides, “Graph Signal Processing–Part I: Graphs, Graph Spectra, and Spectral Clustering”, 2019, _arXiv preprint_ arXiv:1907.03467.
  * [24] J.S. Fabila-Carrasco, F. Lledó, and O. Post, “Matching number, Hamiltonian graphs and magnetic Laplacian matrices,” _Linear Algebra and its Applications_ , 2022, 642, pp. 86-–100.
  * [25] M. Benedicks, and L. Carleson,“The dynamics of the Hénon map,” _Annals of Mathematics_ , 1991, 133(1), pp. 73–169.
  * [26] M.W. Hirsch, S. Smale, and R.L. Devaney, “Differential equations, dynamical systems, and an introduction to chaos”, _Academic press_ , 2012.