# On the Dynamics of Multiagent Nonlinear filtering and Learning

###### Abstract

Multiagent systems aim to accomplish highly complex learning tasks through
decentralised consensus seeking dynamics and their use has garnered a great
deal of attention in the signal processing and computational intelligence
societies. This article examines the behaviour of multiagent networked systems
with nonlinear filtering/learning dynamics. To this end, a general formulation
for the actions of an agent in multiagent networked systems is presented and
conditions for achieving a cohesive learning behaviour is given. Importantly,
application of the so derived framework in distributed and federated learning
scenarios are presented.

Index Terms—  Multiagent systems, nonlinear dynamics, distributed learning,
federated learning,

##  1 Introduction

Traditionally, signal processing and learning techniques have been concerned
with single agent operations [1, 2]. However, as early as the $1970$s, it
became clear that these agents will come into contact with one another and in
most cases cooperation will be beneficial [3, 4]. Moreover, nature provides
exquisite examples of achieving highly complex behaviour from a set of simple,
but collaborating agents [5, 6]. In broad terms, at the heart of most signal
processing and learning techniques lie iterative optimisation approaches that
aim to minimise a performance metric or cost function. Traditionally, these
solutions were formulated with a single information processing unit in mind, a
framework henceforth referred to as centralised, as the optimisation
procedures are performed at a central information processing unit. Although
shown to be effective, these techniques cannot accommodate modern applications
that require interaction of multiple agents over large-scale networks [7, 8].

Initial distributed signal processing and learning techniques were based on
the framework of Pareto optimisation [7, 9] and consensus seeking procedures
[10, 11]. In the learning arena, the federated learning approach, initiated at
Google [12], has gained traction. In this approach, agents update the
parameter set learnt using information available to them locally; then, they
pass these locally updated models to a fusion centre which integrates all the
updates and passes back to the agents a global update [12, 13]. In turn,
federated learning approaches are built upon distributed optimisation
techniques derived for multiagent signal processing tasks [14, 15, 9], where
agents establish a cohesive signal processing procedure among themselves,
using local communications only, without requiring a fusion centre.

The authors have addressed behaviour of linear distributed filtering
techniques in [16, 17], proving that convergence conditions of distributed and
centralised filtering techniques are similar. This article considers the
problem of federated/distributed learning with nonlinear dynamics. To this
end, a general class of adaptive learning techniques is formulated. Then,
operations of these formulated learning technique is decomposed to derive its
federated/distributed dual. Considering discrepancy between performance of
centralised and federated/distributed formulations as a metric, convergence
criteria of the federated/distributed adaptive learning techniques with
nonlinear dynamics is examined. The results significantly relax convergence
criteria for distributed iterative optimisation techniques; thus, opening the
opportunity for using distributed adaptive learning techniques in new
applications areas.

Nomenclature: Scalars, column vectors, and matrices are denoted by lowercase,
bold lowercase, and bold uppercase letters, while remainder of the
nomenclature is summarised as follows:

$\mathbb{N}$ | set of natural numbers
---|---
$\mathcal{N}$ | set of all agents in the network
$\mathcal{N}_{i}$ |  set of agents that can communicate with agent $i$
$|\cdot|$ | cardinality operator
$\|\cdot\|$ | Euclidean norm
$\mathsf{E}\left\\{\cdot\right\\}$ | statistical expectation operator
$(\cdot)^{\mathsf{T}}$ | transpose operator
$\delta(\cdot)$ | Kronecker delta function
$\mathbf{I}$ | identity matrix of appropriate size
$\mathbf{1}$ | vector of appropriate size with unit entries
$\otimes$ | Kronecker product
$\mathsf{diag}\left\\{\cdot\right\\}$ | constructs a block-diagonal matrix from its entries
$\mathsf{p}\left\\{\cdot\right\\}$ | spectral radius operator
$\mathsf{tr}\left\\{\cdot\right\\}$ | trace operator
$\nabla_{\chi}$ |  gradient operator with respect to $\chi$

###  1.1 Problem Formulation

The objective is to estimate/track a vector sequence
$\\{\mathbf{x}_{n}:n\in\mathbb{N}\\}$ representing the condition of a
system111This system can range from a physical system to the state of learning
machines, such as weights of a neural network., that evolves as

$\mathbf{x}_{n+1}=\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})$ (1)

where $\mathrm{f}(\cdot)$ is an evolution function with
$\\{\mathbf{v}_{n}:n\in\mathbb{N}\\}$ representing the evolution noise. To
this end, a set of agents are used that obtain observations modelled as

$\forall i\in\mathcal{N}:\hskip
4.55254pt\mathbf{y}_{i,n}=\mathrm{h}_{i}(\mathbf{x}_{n})+\mathbf{w}_{i,n}$ (2)

where $\mathbf{y}_{i,n}$ and $\mathbf{w}_{i,n}$ are the observation and
observation noise at agent $i$ at time instant $n$, whereas
$\mathrm{h}_{i}(\cdot)$ represents the observation function of agent $i$.

Assumption 1: Observation and state evolution noises are zero-mean white
Gaussian random sequences with joint covariance matrix given by

$\mathsf{E}\left\\{\begin{bmatrix}\mathbf{v}_{n}\\\
\mathbf{w}_{l,n}\end{bmatrix}\begin{bmatrix}\mathbf{v}_{m}\\\
\mathbf{w}_{k,m}\end{bmatrix}^{\mathsf{T}}\right\\}=\begin{bmatrix}\bm{\Sigma}_{\mathbf{v}}&\mathbf{0}\\\
\mathbf{0}&\bm{\Sigma}_{\mathbf{w}_{l}}\delta(l-k)\end{bmatrix}\delta(n-m).$
(3)

On the most fundamental level, this problem can be addressed by organising all
observations into a vector, so that

$\underbrace{\begin{bmatrix}\mathbf{y}_{1,n}\\\ \vdots\\\
\mathbf{y}_{|\mathcal{N}|,n}\end{bmatrix}}_{\mathbf{y}_{n}}=\underbrace{\begin{bmatrix}\mathrm{h}_{1}(\mathbf{x}_{n})\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}(\mathbf{x}_{n})\end{bmatrix}}_{\mathrm{h}(\mathbf{x}_{n})}+\underbrace{\begin{bmatrix}\mathbf{w}_{1,n}\\\
\vdots\\\ \mathbf{w}_{|\mathcal{N}|,n}\end{bmatrix}}_{\mathbf{w}_{n}}$

which allows for the use of centralised adaptation operations [18]

$\hat{\mathbf{x}}_{n}=\mathrm{f}(\hat{\mathbf{x}}_{n-1})+\mathbf{G}_{n}\Big{(}\mathbf{y}_{n}-\mathrm{h}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{n-1}\right)\right)\Big{)}$
(4)

where $\mathbf{G}_{n}$ is a gain matrix calculated to minimise an error
measure using a gradient-based format and $\hat{\mathbf{x}}_{n}$ denotes the
estimate of $\mathbf{x}_{n}$, that is,
$\hat{\mathbf{x}}_{n}=\mathsf{E}\left\\{\mathbf{x}_{n}|\mathbf{y}_{1},\ldots,\mathbf{y}_{n}\right\\}$,
representing the best approximation of the expected value of $\mathbf{x}_{n}$
given the available information, while $\mathrm{f}(\hat{\mathbf{x}}_{n-1})$ is
projection of the estimate form time instant $n-1$ onto time instant $n$.

Assumption 2: The estimation error of the centralised filter in (4) at time
instant $n$ is defined as

$\bm{\epsilon}_{n}=\mathbf{x}_{n}-\hat{\mathbf{x}}_{n}$ (5)

where it is assumed the matrix sequence $\\{\mathbf{G}_{n}:n\in\mathbb{N}\\}$
exist so that the filtering operations in (4) become exponentially bounded in
the mean square error (MSE) sense. That is, for some $\mu>0$ there is
$\eta,\nu>0$ and $\rho\in(0,1)$ so that

$\text{if}\hskip 3.41418pt\|\bm{\epsilon}_{1}\|^{2}\leq\mu\text{; then,}\hskip
3.41418pt\mathsf{E}\left\\{\|\bm{\epsilon}_{n}\|^{2}\right\\}\leq\eta\|\bm{\epsilon}_{1}\|^{2}\rho^{n}+\nu.$
(6)

###  1.2 From Centralised to Federated Learning

Focusing on the expression in (4) and partitioning the gain matrix into
appropriate sections so that one has

$\mathbf{G}_{n}=\begin{bmatrix}\mathbf{G}_{1,n},\ldots,\mathbf{G}_{|\mathcal{N}|,n}\end{bmatrix}$

allows (4) to be expressed as

$\displaystyle\hat{\mathbf{x}}_{n}=$
$\displaystyle\mathrm{f}(\hat{\mathrm{x}}_{n-1})+\sum_{\forall
i\in\mathcal{N}}\mathbf{G}_{i,n}\Big{(}\mathbf{y}_{i,n}-\mathrm{h}_{i}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{n-1}\right)\right)\Big{)}$
(7) $\displaystyle=$
$\displaystyle\underbrace{\frac{1}{|\mathcal{N}|}\sum_{\forall
i\in\mathcal{N}}}_{\text{information
fusion}}\underbrace{\left(\mathrm{f}(\hat{\mathbf{x}}_{n-1})+|\mathcal{N}|\mathbf{G}_{i,n}\Big{(}\mathbf{y}_{i,n}-\mathrm{h}_{i}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{n-1}\right)\right)\right)}_{\text{local
filtering operation}}.$

The expression in (7) essentially partitions the operation of (4) into two
processes; i) local filtering operations using the observation information
available at the individual agent, labelled as “local filtering operation” in
(7), ii) an information fusion operation that fusses local learning operations
into a cohesive global operation, labelled as “information fusion” in (7).

Remark 1: The formulation in (7) forms the cornerstone of federated learning
techniques [19, 14, 13]. In addition, note that only agents with new data sets
at each time instant have to send updates to the fusion centre, while agents
that have no new data sets only implement the projection operation. Moreover,
subsets of agents can be selected to communicate with the fusion centre to
avoid communication traffic. However, all these scenarios are special cases of
the presented framework.

###  1.3 From Federated to Distributed Learning

In distributed learning scenarios, agents are connected through an ad-hoc
network, that is modelled as a strongly connected graph. In this setting,
agent $l$ can only communicate with its neighbouring agents denoted as
$\mathcal{N}_{l}$, that includes agent $l$ itself. Thus, the information
fusion step in (7) is approximated in a distributed fashion as

$\displaystyle\bm{\phi}_{i,n}=$
$\displaystyle\mathrm{f}(\hat{\mathbf{x}}_{i,n-1})+\tilde{\mathbf{G}}_{i,n}\Big{(}\mathbf{y}_{i,n}-\mathrm{h}_{i}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{i,n-1}\right)\right)\Big{)}$
(8a) $\displaystyle\hat{\mathbf{x}}_{i,n}=$
$\displaystyle\bm{\phi}_{i,n}+\sum_{\forall
j\in\mathcal{N}_{i}}c_{i,j}\left(\bm{\phi}_{i,n}-\bm{\phi}_{j,n}\right)$ (8b)

where agent $i\in\mathcal{N}$ implements a local filtering operation, i.e.,
(8a), using its local gain matrix $\tilde{\mathbf{G}}_{i,n}$, yielding the
intermediate estimate $\bm{\phi}_{i,n}$; then, it aims to approximate the
information fusion operation in (7) whereby intermediate estimates are fused
with those of its neighbours in (8b) using combination weights
$\\{c_{i,j}:\forall j\in\mathcal{N}_{i}\\}$ which yields the final estimate,
$\hat{\mathbf{x}}_{i,n}$.

Assumption 3: In keeping with convention, for both diffusion [20] and
consensus [21, 22] frameworks, it is assumed that the network remains
connected222That is, there exists an path between any two agents. and

$\forall i\in\mathcal{N}:\hskip 3.41418pt\sum_{\forall
j\in\mathcal{N}_{i}}c_{i,j}=1.$ (9)

Therefore, matrix $\mathbf{C}$ with $i^{\text{th}}$ row and $j^{\text{th}}$
column element

$\mathbf{C}^{\\{i,j\\}}=\left\\{\begin{matrix}c_{i,j}&{}&\text{if}\hskip
3.41418ptj\in\mathcal{N}_{i}\\\ 0&{}&\text{otherwise}\end{matrix}\right.$

will be row-stochastic and primitive [7].

Remark 2: There are two main techniques for fully distributed information
fusion, diffusion [20] and consensus [22, 21] with the formulation in (7)
accommodating both.

##  2 Performance Discrepancy and Stability

The estimation error at agent $i$ is defined as

$\bm{\epsilon}_{i,n}=\mathbf{x}_{n}-\hat{\mathbf{x}}_{i,n}$ (10)

while the network-wide perspective of the estimation error of all agents is
given by

$\mathbfcal{E}_{n}=\begin{bmatrix}\bm{\epsilon}^{\mathsf{T}}_{1,n},\ldots,\bm{\epsilon}^{\mathsf{T}}_{|\mathcal{N}|,n}\end{bmatrix}^{\mathsf{T}}.$
(11)

Taking performance of the centralised adaptation process in (4) as the
benchmark, performance discrepancy between (4) and (8) can now be formulated
as

$\displaystyle\bm{\Delta}_{n+1}=$
$\displaystyle\mathbfcal{E}_{n+1}-\mathbf{1}\otimes\bm{\epsilon}_{n+1}$ (12)
$\displaystyle=$
$\displaystyle\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}\begin{bmatrix}\mathrm{f}(\hat{\mathbf{x}}_{n})\\\
\vdots\\\
\mathrm{f}(\hat{\mathbf{x}}_{n})\end{bmatrix}-\mathbfcal{C}\begin{bmatrix}\mathrm{f}(\hat{\mathbf{x}}_{1,n})\\\
\vdots\\\ \mathrm{f}(\hat{\mathbf{x}}_{|\mathcal{N}|,n})\end{bmatrix}$
$\displaystyle+\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}\mathbfcal{G}_{n+1}-\mathbfcal{C}\tilde{\mathbfcal{G}}_{n+1}\right)\mathbf{w}_{n+1}$
$\displaystyle+\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}\mathbfcal{G}_{n+1}-\mathbfcal{C}\tilde{\mathbfcal{G}}_{n+1}\right)\mathrm{h}\left(\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})\right)$
$\displaystyle-\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}\mathbfcal{G}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\end{bmatrix}+\mathbfcal{C}\tilde{\mathbfcal{G}}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\hat{\mathbf{x}}_{1,n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\hat{\mathbf{x}}_{|\mathcal{N}|,n})\right)\end{bmatrix}$

where $\textbf{\emph{1}}=\mathbf{1}\mathbf{1}^{\mathsf{T}}\otimes\mathbf{I}$,
$\mathbfcal{G}_{n}=\mathsf{diag}\left\\{|\mathcal{N}|\mathbf{G}_{i,n}:\forall
i\in\mathcal{N}\right\\}$, and
$\tilde{\mathbfcal{G}}_{n}=\mathsf{diag}\left\\{\tilde{\mathbf{G}}_{i,n}:\forall
i\in\mathcal{N}\right\\}$. Note that the derivation of (12) is given in
Supplementary Material SM-1 [23]. The expression in (12) shows performance
discrepancy between distributed and centralised approaches in relation to the
network typology, gain matrices, and observation/evolution functions.

###  2.1 The Case of Federated Learning

In order to formulate and isolate the effect of replacing federated
information fusion in (7) with that of the distributed information fusion step
in (8b), it is assumed that both approaches use the equivalent gain matrices,
that is,

$\forall i\in\mathcal{N}:\hskip
3.41418pt\tilde{\mathbf{G}}_{i,n}=|\mathcal{N}|\mathbf{G}_{i,n}.$

Now, from (12) and $\tilde{\mathbfcal{G}}_{n+1}$ one has

$\displaystyle\bm{\Delta}_{n+1}=$
$\displaystyle\mathbfcal{C}\left(\begin{bmatrix}\mathrm{f}(\hat{\mathbf{x}}_{n})\\\
\vdots\\\
\mathrm{f}(\hat{\mathbf{x}}_{n})\end{bmatrix}-\begin{bmatrix}\mathrm{f}(\hat{\mathbf{x}}_{1,n})\\\
\vdots\\\ \mathrm{f}(\hat{\mathbf{x}}_{|\mathcal{N}|,n})\end{bmatrix}\right)$
(13)
$\displaystyle-\mathbfcal{C}\tilde{\mathbfcal{G}}_{n+1}\left(\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\end{bmatrix}-\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\hat{\mathbf{x}}_{1,n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\hat{\mathbf{x}}_{|\mathcal{N}|,n})\right)\end{bmatrix}\right)$
$\displaystyle+\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}-\mathbfcal{C}\right)\tilde{\mathbfcal{G}}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})\right)\end{bmatrix}$
$\displaystyle-\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}-\mathbfcal{C}\right)\tilde{\mathbfcal{G}}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\end{bmatrix}$
$\displaystyle+\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}-\mathbfcal{C}\right)\tilde{\mathbfcal{G}}_{n+1}\mathbf{w}_{n+1}$

with the full derivation given in Supplementary Material SM-2 [23]. Now,
consider the case where at time instant $m\in\mathbb{N}$ there exists
$\mu^{\prime}>0$ so that
$\|\hat{\mathbf{x}}_{m}-\hat{\mathbf{x}}_{i,m}\|^{2}\leq\mu^{\prime}$
resulting in

$\forall i\in\mathcal{N}:\left\\{\begin{aligned}
&\mathrm{h}_{i}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{m}\right)\right)-\mathrm{h}_{i}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{i,m}\right)\right)=\mathbf{H}_{i}\mathbf{A}\bm{\Delta}_{m}+\mathrm{Res}_{\mathrm{h}_{i},\mathrm{f}}\\\
&\\\
&\mathrm{f}(\hat{\mathbf{x}}_{m})-\mathrm{f}(\hat{\mathbf{x}}_{i,m})=\mathbf{A}\bm{\Delta}_{m}+\mathrm{Res}_{\mathrm{f}}\end{aligned}\right.$
(14)

where $\mathbf{A}$ and $\mathbf{H}_{i}$ denote the first-order approximation
of $\mathrm{f}(\cdot)$ and $\mathrm{h}_{i}(\cdot)$, while
$\mathrm{Res}_{\mathrm{f}}$ and $\mathrm{Res}_{\mathrm{h}_{i},\mathrm{f}}$
denote the residuals. From (13) and (14) it follows that

$\displaystyle\mathsf{E}\left\\{\|\bm{\Delta}_{m+1}\|^{2}\right\\}\leq$
$\displaystyle\mathsf{p}\left\\{\mathbfcal{C}\left(\mathbf{I}-\tilde{\mathbfcal{G}}_{n+1}\mathbfcal{H}\right)\mathbfcal{A}\right\\}\mathsf{E}\left\\{\|\bm{\Delta}_{m}\|^{2}\right\\}$
(15)
$\displaystyle+\mathsf{p}\left\\{\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}-\mathbfcal{C}\right)\tilde{\mathbfcal{G}}_{n+1}\right\\}\mathsf{tr}\left\\{{\bm{\Sigma}}_{\mathbf{w}}\right\\}$
$\displaystyle+\|\mathrm{Res}_{\mathrm{f}}\|^{2}+\sum_{\forall
i\in\mathcal{N}}\|\mathrm{Res}_{\mathrm{h}_{i},\mathrm{f}}\|^{2}$
$\displaystyle+\mathsf{p}\left\\{\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}-\mathbfcal{C}\right)\tilde{\mathbfcal{G}}_{n+1}\right\\}\|\mathrm{Res}_{\text{ob}}\|^{2}$

where $\mathbfcal{H}=\mathsf{diag}\left\\{\mathbf{H}_{i}:\forall
i\in\mathcal{N}\right\\}$, $\mathbfcal{A}=\mathbf{I}\otimes\mathbf{A}$,
$\bm{\Sigma}_{\mathbf{w}}=\mathsf{E}\left\\{\mathbf{W}_{n}\mathbf{W}^{\mathsf{T}}_{n}\right\\}$,

$\mathrm{Res}_{\text{ob}}=\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})\right)\end{bmatrix}-\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\end{bmatrix}.$

This allows for the following results to be formulated

Lemma 1: If Assumption 3 holds; then,

$\mathsf{p}\left\\{\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}-\mathbfcal{C}\right\\}<1.$

Proof of Lemma 1: See the Supplementary Material, SM-3 [23].

Theorem 1: If $\forall i\in\mathcal{N}$ the information set
$\\{\mathbf{y}_{j,m}:j\in\mathcal{N}_{i},m\leq n\\}$ is sufficient to allow
agent $i$ to find a gain matrix that satisfies (6). Then,
$\|\bm{\Delta}_{n}\|^{2}$ as given in (15) is bounded.

Proof of Theorem 1: See Supplementary Material SM-4 [23].

Remark 3: The results in Theorem 1, expand on previous works in this area [24,
25].

Remark 3: In general, (15) indicates that combination weights and gains should
be selected in conjunction if convergence is to be achieved in wider
scenarios. Although the results form Theorem 1, expand on previous works in
this aria [24, 25], in the following section more general scenarios are
considered.

Remark 4: From (15), note that $\|\bm{\Delta}_{m}\|^{2}$ finds a minimum as
$\mathbfcal{C}\rightarrow\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}$. This can
be achieved through the a consensus framework [22, 26], with a communication
complexity penalty for any required level of accuracy [19, 27].

###  2.2 The Case of Distributed Learning

The aim of each agent, $i\in\mathcal{N}$, is to implement filtering operations
in (8a), where the matrix $\tilde{\mathbf{G}}_{i,n}$ is calculated so as to
implement the fundamental filtering process so that

$\bm{\phi}_{i,n+1}=\mathsf{E}\left\\{\mathbf{x}_{n}|\hat{\mathbf{x}}_{j,m},\mathbf{y}_{i,m}:m\leq
n-1,\forall j\in\mathcal{N}_{i}\right\\}.$ (16)

The formulation in (16) accommodates estimates shared with the agent from its
neighbourhood, that is, $\\{\hat{\mathbf{x}}_{j,m}:m\leq n-1,\forall
j\in\mathcal{N}_{i}\\}$, that according to (8) incorporate information from
estimates of their neighbours, as well as, their observations. Thus, modelling
how information is generated and flows through the network.

Theorem 2: If information set
$\\{\hat{\mathbf{x}}_{j,m},\mathbf{y}_{i,m}:m\leq n-1,\forall
j\in\mathcal{N}_{i}\\}$, yields gain matrix $\tilde{\mathbf{G}}_{i,n}$ so that

$\forall i\in\mathcal{N}:\hskip
9.10509pt\|\mathbf{x}_{n}-\bm{\phi}_{i,n}\|^{2}\leq\gamma_{i}\|\underbrace{\mathbf{x}_{n-1}-\hat{\mathbf{x}}_{i,n-1}}_{\bm{\epsilon}_{i,n-1}}\|^{2}$
(17)

where $\gamma_{i}\in(0,1)$; then, $\|\bm{\Delta}_{n}\|^{2}$ becomes bounded.

Proof of Theorem 2: See Supplementary Material SM-5 [23].

Remark 3: There are similar results available for linear [28] and nonlinear
[25, 24, 9] distributed filtering techniques, where convergence criteria are
either set at each agent having sufficient observation information to converge
[25, 24] or this information to be shared with them through neighbouring
agents [28]. However, Theorem 2 is more comprehensive, as its setting allows
introduction of more elaborate scenarios.

A more useful and widely encountered scenario is where the goal of each agent
remains as given in (16), whereas the condition in (17) only holds for a
subset of elements in $\bm{\phi}_{i,n}$. This setting is formulated as

$\bm{\phi}^{c}_{i,n}=\bm{\Upsilon}_{i}\bm{\phi}_{i,n}\hskip
6.82881pt\text{and}\hskip
6.82881pt\bm{\phi}^{c^{\prime}}_{i,n}=\bm{\Upsilon}^{\prime}_{i}\bm{\phi}_{i,n}$
(18)

where $\bm{\Upsilon}_{i}$ (cf. $\bm{\Upsilon}^{\prime}_{i}$) passes elements
of $\bm{\phi}_{i,n}$ that meet (cf. do not meet) the condition in (17)
unchanged, while zeroing other elements. From the definition in (18) it
follows that

$\forall i\in\mathcal{N}:\hskip
3.41418pt\bm{\phi}_{i,n}=\bm{\phi}^{c}_{i,n}+\bm{\phi}^{c^{\prime}}_{i,n}\hskip
6.82881pt\text{and}\hskip
6.82881pt\bm{\Upsilon}_{i}+\bm{\Upsilon}^{\prime}_{i}=\mathbf{I}$

and therefore, from (8b), we have

$\bm{\Delta}_{n+1}=\mathbfcal{C}\begin{bmatrix}\hat{\mathbf{x}}_{n+1}-\bm{\phi}^{c}_{1,n+1}\\\
\vdots\\\
\hat{\mathbf{x}}_{n+1}-\bm{\phi}^{c}_{|\mathcal{N}|,n+1}\end{bmatrix}+\mathbfcal{C}\begin{bmatrix}\hat{\mathbf{x}}_{n+1}-\bm{\phi}^{c^{\prime}}_{1,n+1}\\\
\vdots\\\
\hat{\mathbf{x}}_{n+1}-\bm{\phi}^{c^{\prime}}_{|\mathcal{N}|,n+1}\end{bmatrix}.$
(19)

Assumption 4: Since, the elements of $\bm{\phi}^{c^{\prime}}_{i,n}$ are not
observable at agent $i$, it is assumed that the local filtering operation in
(8a) can be simplified to

$\bm{\phi}^{c^{\prime}}_{i,n}=\bm{\Upsilon}^{{}^{\prime}}_{i}\mathrm{f}\left(\hat{\mathbf{x}}_{i,n-1}\right)$

that is, these elements are only assumed to be projections from the previous
time instant, as there is not enough information in the observations to update
them.

Theorem 3: Given
$g_{\bm{\Delta}_{m}}=\mathsf{p}\left\\{\mathbfcal{C}\left(\mathbf{I}-\bm{\Upsilon}\tilde{\mathbfcal{G}}_{m+1}\mathbfcal{H}\right)\mathbfcal{A}\right\\}$,
if

$\exists k\in\mathbb{N}\hskip 3.41418pt\text{so~{}that}\hskip
3.41418pt\left(\prod^{k}_{i=0}g_{\bm{\Delta}_{m-i}}\right)<1$

then, $\|\bm{\Delta}_{n}\|^{2}$ becomes bounded.

Proof of Theorem 3: See Supplementary Material SM-6 [23].

Remark 4: Theorem 3 relates convergence to the information available to the
agent, the dynamics of the system, that is, $\mathrm{f}(\cdot)$ and
$\\{\mathrm{h}_{i}(\cdot):i\in\mathcal{N}\\}$, and the mechanism with which
information flows through the network. In essence, this condition states that
for each element of $\bm{\Phi}^{c^{\prime}}_{i,n}$ that cannot be learnt from
the information available to agent $i$, there must exist agents in the $k$-hop
neighbourhood of agent $i$ to which this information is available. The
parameter $k$ itself, determined by the state-space dynamics.

##  3 Numerical Example

For simulations, the network of $120$ agents were used with structure of the
network as given in Supplementary Material SM-7 [23].

###  3.1 Statistical Learning

The most fundamental learning paradigm is considered, that is, each agent
$i\in\mathcal{N}$ aims to learn the parameter $\hat{\mathbf{x}}_{i,n}$, so
that

$\hat{\mathbf{y}}_{i,n}=\mathrm{h}_{i,n}(\hat{\mathbf{x}}_{i,n})=\varphi(\mathbf{Z}_{i,n}\hat{\mathbf{x}}_{i,n})$
(20)

is a best fit for observed information $\\{\mathbf{y}_{i,m}:\forall
i\in\mathcal{N},m<n\\}$, where $\mathbf{Z}_{i,n}$ is a stochastic signal at
agent $i$. This is achieved via solutions of the form in (8), where
$\mathrm{f}(\cdot)$ is an identity function, $\mathrm{h}_{i,n}(\cdot)$ is
formulated in (20), and
$\tilde{\mathbf{G}}_{i,n}=\zeta\nabla_{\hat{\mathbf{x}}_{i,n}}\left(\|\mathbf{y}_{i,n}-\hat{\mathbf{y}}_{i,n}\|^{2}\right)$
with $\zeta\in(0,1)$. Agents seek the Pareto solution to the optimisation task

$\min\mathsf{E}\left\\{\frac{1}{|\mathcal{N}|}\sum_{\forall
i\in\mathcal{N}}\|\mathbf{y}_{i,n}-\hat{\mathbf{y}}_{i,n}\|^{2}\right\\}\cdot$

For simulations, $\\{\mathbf{Z}_{i,n}:\forall
i\in\mathcal{N},n\in\mathbb{N}\\}$ was a set of zero-mean white Gaussian
process and the nonlinearity, $\varphi(\cdot)$, used was the hyperbolic
tangent function. As a worst case scenario, it was considered that information
required for convergence was only available to the collective of the agents
That is, no one agent or neighbourhood of an gent had enough information for
convergence. This falls beyond the convergence criteria set previously;
however, the conditions satisfy Theorem 3\. The MSE performance of each agent
is benchmarked agents that of centralised learning in Fig. 1. Notice that once
more all agents exhibited convergent behaviour, validating the introduced
concepts in this article.

Fig. 1: MSE performance of distributed and centralised learning solutions.
Mean averages were calculated over $2000$ independent realisations. The region
in blue indicates where MSE of each agent lies for the distributed approach.

###  3.2 State-Space Filtering

The problem of tracking a state-space system, with the nonlinear dynamics
given in [29], where the main objective was to track a target manoeuvring
inside a two-dimensional box, was considered. The problem is referred to as
“particle-in-the-box” usually assuming measurements of the target location is
available at each agent. However, in the simulation example here, it was
assumed that only one agent had position measurements on the vertical axis,
while other agents could only acquire position measurements on the horizontal
axes. Note, that this setting presents a worst case scenario for distribute
filtering. For the sake of simulations, the sampling rate was set at $25$ Hz,
while the state evolution and observation noise covariances were
$0.04\times\mathbf{I}$ and $0.16$. The distributed extended Kalman filter in
[25, 24] was used to collaboratively track the target. Although these
conditions fall outside the convergence requirements originally proposed in
[25, 24], the results in Fig. 2 indicate all agents tracked the target with
reasonable accuracy validating the theoretical concepts in Theorem 3.

Fig. 2: Estimation error for the state of particle in the box. Estimation
error of all elements of the state vector at all agents are shown.

##  4 Conclusion

The dynamics of distributed filtering and learning have been revisited from a
novel perspective. The focus has been on the amount of information available
to each agent, implementing iterative optimisation process in order to learn
desired parameters from observed signals. This has allowed the formulation of
a number of scenarios for distributed learning and filtering the convergence
and stability of which falling beyond the scope of former studies.

## References

  * [1] T. Kailath, A. H. Sayed, and B. Hassibi,  Linear estimation,  Prentice Hall, 2000.
  * [2] D. P. Mandic and J. A. Chambers,  Recurrent Neural Networks for Prediction: Learning Algorithms and Stability,  John Wiley & Sons, LTD, 2001.
  * [3] J. Speyer,  “Computation and transmission requirements for a decentralized linear-quadratic-Gaussian control problem,”  IEEE Transactions on Automatic Control, vol. AC-24, no. 2, pp. 266–269, April 1979.
  * [4] N. Sandell, P. Varaiya, M. Athans, and M. Safonov,  “Survey of decentralized control methods for large scale systems,”  IEEE Transactions on Automatic Control, vol. 23, no. 2, pp. 108–128, 1978.
  * [5] A. H. Sayed, S. Y. Tu, J. Chen, X. Zhao, and Z. J. Towfic,  “Diffusion strategies for adaptation and learning over networks: An examination of distributed strategies and network behavior,”  IEEE Signal Processing Magazine, vol. 30, no. 3, pp. 155–171, 2013\.
  * [6] R. Olfati-Saber,  “Flocking for multi-agent dynamic systems: Algorithms and theory,”  IEEE Transactions on Automatic Control, vol. 51, no. 3, pp. 401–420, March 2006.
  * [7] Ali H. Sayed,  “Adaptation, learning, and optimization over networks,”  Foundations and Trends${}^{\text{\tiny{\textregistered}}}$ in Machine Learning, vol. 7, no. 4-5, pp. 311–801, 2014.
  * [8] L. Busoniu, R. Babuska, and B. De Schutter,  “A comprehensive survey of multiagent reinforcement learning,”  IEEE Transactions on Systems, Man, and Cybernetics, Part-C (Applications and Reviews), vol. 38, no. 2, pp. 156–172, 2008.
  * [9] S. Vlaski and A. H. Sayed,  “Distributed learning in non-convex environments-part i: Agreement at a linear rate,”  IEEE Transactions on Signal Processing, vol. 69, pp. 1242–1256, 2021.
  * [10] M. H. DeGroot,  “Reaching a consensus,”  Journal of American Statististical Association, vol. 69, no. 345, pp. 118–121, 1974.
  * [11] Ondrej Sluvciak,  Convergence Analysis of Distributed Consensus Algorithms,  Ph.D. Thesis, Vienna University of Technology, 2013.
  * [12] J. Konecny, H. B. McMahan, D. Ramage, and P. Richtarik,  “Federated optimization: Distributed machine learning for on-device intelligence,”  arXiv:1610.02527, 2016.
  * [13] T. Li, A. K. Sahu, A. Talwalkar, and V. Smith,  “Federated learning: Challenges, methods, and future directions,”  IEEE Signal Processing Magazine, vol. 37, no. 3, pp. 50–60, 2020\.
  * [14] Angelia Nedic,  “Distributed gradient methods for convex machine learning problems in networks: Distributed optimization,”  IEEE Signal Processing Magazine, vol. 37, no. 3, pp. 92–101, 2020\.
  * [15] S. Hosseinalipour, C. G. Brinton, V. Aggarwal, H. Dai, and Mung M. Chiang,  “From federated to fog learning: Distributed machine learning over heterogeneous wireless networks,”  IEEE Communications Magazine, vol. 58, no. 12, pp. 41–47, 2020\.
  * [16] Sayed Pouria Talebi, Stefan Werner, Yih-Fang Huang, and Vijay Gupta,  “Distributed algebraic Riccati equations in multi-agent systems,”  In Procceedings of European Control Conference, pp. 1810–1817, 2022\.
  * [17] Sayed Pouria Talebi, Stefan Werner, Vijay Gupta, and Yih-Fang Huang,  “On stability and convergence of distributed filters,”  IEEE Signal Processing Letters, vol. 28, pp. 494–498, 2021.
  * [18] S. Särkkä,  Bayesian Filtering and Smoothing,  Cambridge University Press, 2013.
  * [19] S. P. Talebi and S. Werner,  “Distributed Kalman filtering and control through embedded average consensus information fusion,”  IEEE Transactions on Automatic Control, vol. 64, no. 10, pp. 4396–4403, 2019.
  * [20] A. H. Sayed and C. G. Lopes,  “Adaptive processing over distributed networks,”  IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences, vol. E90-A, no. 8, pp. 1504–1510, August 2007.
  * [21] Lin Xiao and Stephen Boyd,  “Fast linear iterations for distributed averaging,”  Systems & Control Letters, vol. 53, no. 1, pp. 65–78, 2004.
  * [22] R. Olfati-Saber and J. S. Shamma,  “Consensus filters for sensor networks and distributed sensor fusion,”  In Proceedings of IEEE Conference on Decision and Control, December 2005.
  * [23] S. P. Talebi and D. P. Mandic,  “On the dynamics of mulitagent nonlinear filtering and learning,”  arXiv.
  * [24] S. F. Cattivelli and A. H. Sayed,  “Distributed nonlinear Kalman filtering with applications to wireless localization,”  In Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, pp. 3522–3525, 2010.
  * [25] S. Kanna and D. P. Mandic,  “Stability of distributed extended Kalman filters,”  In Proceedings of International Conference on Digital Signal Processing (DSP), pp. 1–5, 2017.
  * [26] L. Xiao and S. Boyd,  “Fast linear iterations for distributed averaging,”  Systems & Control Letters, vol. 53, no. 1, pp. 65–78, 2004.
  * [27] S. P. Talebi, S. Werner, and D. P. Mandic,  “Quaternion-valued distributed filtering and control,”  IEEE Transactions on Automatic Control, vol. 65, no. 10, pp. 4246–4257, 2020.
  * [28] F. S. Cattivelli and A. H. Sayed,  “Diffusion strategies for distributed Kalman filtering and smoothing,”  IEEE Transactions on Automatic Control, vol. 55, no. 9, pp. 2069–2084, September 2010.
  * [29] R. Olfati-Saber and Parisa Jalalkamali,  “Collaborative target tracking using distributed Kalman filtering on mobile sensor networks,”  In Proceedings American Control Conference, pp. 1100–1105, 2011\.
  * [30] K. Reif, S. Gunther, E. Yaz, and R. Unbehauen,  “Stochastic stability of the discrete-time extended Kalman filter,”  IEEE Transactions on Automatic Control, vol. 44, no. 4, pp. 714–728, 1999.

## Supplementary Material

### SM-1:Derivation of (12)

From substituting (1) and (8) into (10), one has

$\displaystyle\bm{\epsilon}_{i,n+1}=$
$\displaystyle\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})-\sum_{\forall
j\in\mathcal{N}_{i}}c_{i,j}\mathrm{f}(\hat{\mathbf{x}}_{j,n})$
$\displaystyle-\sum_{\forall
j\in\mathcal{N}_{i}}c_{i,j}\tilde{\mathbf{G}}_{j,n+1}\Big{(}\mathbf{y}_{j,n+1}-\mathrm{h}_{j}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{j,n}\right)\right)\Big{)}.$

that can be reformulated using (9) to yield

$\displaystyle\bm{\epsilon}_{i,n+1}=$ $\displaystyle\sum_{\forall
j\in\mathcal{N}_{i}}c_{i,j}\Big{(}\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})-\mathrm{f}(\hat{\mathbf{x}}_{j,n})\Big{)}$
(21) $\displaystyle-\sum_{\forall
j\in\mathcal{N}_{i}}c_{i,j}\tilde{\mathbf{G}}_{j,n+1}\Big{(}\mathrm{h}_{j}\left(\mathrm{f}\left(\mathbf{x}_{n},\mathbf{v}_{n}\right)\right)-\mathrm{h}_{j}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{j,n}\right)\right)\Big{)}$
$\displaystyle-\sum_{\forall
j\in\mathcal{N}_{i}}c_{i,j}\tilde{\mathbf{G}}_{j,n+1}\mathbf{w}_{j,n+1}.$

Thus, the network-wide error in (11) can be expressed as

$\displaystyle\mathbfcal{E}_{n+1}=$
$\displaystyle\mathbfcal{C}\begin{bmatrix}\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})-\mathrm{f}(\hat{\mathbf{x}}_{1,n})\\\
\vdots\\\
\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})-\mathrm{f}(\hat{\mathbf{x}}_{|\mathcal{N}|,n})\end{bmatrix}-\mathbfcal{C}\tilde{\mathbfcal{G}}_{n+1}\mathbfcal{W}_{n+1}$
(22)
$\displaystyle-\mathbfcal{C}\tilde{\mathbfcal{G}}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}\left(\mathbf{x}_{n},\mathbf{v}_{n}\right)\right)-\mathrm{h}_{1}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{1,n}\right)\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}\left(\mathbf{x}_{n},\mathbf{v}_{n}\right)\right)-\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{|\mathcal{N}|,n}\right)\right)\end{bmatrix}$

where $\mathbfcal{C}=\mathbf{C}\otimes\mathbf{I}$,
$\tilde{\mathbfcal{G}}_{n}=\mathsf{diag}\left\\{\tilde{\mathbf{G}}_{i,n}:\forall
i\in\mathcal{N}\right\\}$, and

$\mathbfcal{W}_{n}=\begin{bmatrix}\mathbf{w}^{\mathsf{T}}_{1,n},\ldots,\mathbf{w}^{\mathsf{T}}_{|\mathcal{N}|,n}\end{bmatrix}^{\mathsf{T}}.$

Now, deducting (5) from (22) gives

$\displaystyle\bm{\Delta}_{n+1}=$
$\displaystyle\mathbfcal{E}_{n+1}-\mathbf{1}\otimes\bm{\epsilon}_{n+1}$
$\displaystyle=$
$\displaystyle\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}\begin{bmatrix}\mathrm{f}(\hat{\mathbf{x}}_{n})\\\
\vdots\\\
\mathrm{f}(\hat{\mathbf{x}}_{n})\end{bmatrix}-\mathbfcal{C}\begin{bmatrix}\mathrm{f}(\hat{\mathbf{x}}_{1,n})\\\
\vdots\\\ \mathrm{f}(\hat{\mathbf{x}}_{|\mathcal{N}|,n})\end{bmatrix}$
$\displaystyle+\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}\mathbfcal{G}_{n+1}-\mathbfcal{C}\tilde{\mathbfcal{G}}_{n+1}\right)\mathbf{w}_{n+1}$
$\displaystyle+\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}\mathbfcal{G}_{n+1}-\mathbfcal{C}\tilde{\mathbfcal{G}}_{n+1}\right)\mathrm{h}\left(\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})\right)$
$\displaystyle-\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}\mathbfcal{G}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\end{bmatrix}$
$\displaystyle+\mathbfcal{C}\tilde{\mathbfcal{G}}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\hat{\mathbf{x}}_{1,n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\hat{\mathbf{x}}_{|\mathcal{N}|,n})\right)\end{bmatrix}$

with
$\mathbfcal{G}_{n}=\mathsf{diag}\left\\{|\mathcal{N}|\mathbf{G}_{i,n}:\forall
i\in\mathcal{N}\right\\}$ and
$\textbf{\emph{1}}=\mathbf{1}\mathbf{1}^{\mathsf{T}}\otimes\mathbf{I}$.

### SM-2: Derivation of (13)

Substituting, all gain related matrices in (12) with
$\tilde{\mathbfcal{G}}_{n+1}$ yields

$\displaystyle\bm{\Delta}_{n+1}=$
$\displaystyle\mathbfcal{C}\left(\begin{bmatrix}\mathrm{f}(\hat{\mathbf{x}}_{n})\\\
\vdots\\\
\mathrm{f}(\hat{\mathbf{x}}_{n})\end{bmatrix}-\begin{bmatrix}\mathrm{f}(\hat{\mathbf{x}}_{1,n})\\\
\vdots\\\ \mathrm{f}(\hat{\mathbf{x}}_{|\mathcal{N}|,n})\end{bmatrix}\right)$
$\displaystyle-\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}\tilde{\mathbfcal{G}}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\end{bmatrix}$
$\displaystyle+\mathbfcal{C}\tilde{\mathbfcal{G}}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\hat{\mathbf{x}}_{1,n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\hat{\mathbf{x}}_{|\mathcal{N}|,n})\right)\end{bmatrix}$
$\displaystyle+\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}-\mathbfcal{C}\right)\tilde{\mathbfcal{G}}_{n+1}\mathbf{w}_{n+1}$
$\displaystyle+\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}-\mathbfcal{C}\right)\tilde{\mathbfcal{G}}_{n+1}\mathrm{h}\left(\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})\right)$

which can be rearranged into

$\displaystyle\bm{\Delta}_{n+1}=$
$\displaystyle\mathbfcal{C}\left(\begin{bmatrix}\mathrm{f}(\hat{\mathbf{x}}_{n})\\\
\vdots\\\
\mathrm{f}(\hat{\mathbf{x}}_{n})\end{bmatrix}-\begin{bmatrix}\mathrm{f}(\hat{\mathbf{x}}_{1,n})\\\
\vdots\\\ \mathrm{f}(\hat{\mathbf{x}}_{|\mathcal{N}|,n})\end{bmatrix}\right)$
(23)
$\displaystyle-\mathbfcal{C}\tilde{\mathbfcal{G}}_{n+1}\left(\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\end{bmatrix}-\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\hat{\mathbf{x}}_{1,n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\hat{\mathbf{x}}_{|\mathcal{N}|,n})\right)\end{bmatrix}\right)$
$\displaystyle+\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}-\mathbfcal{C}\right)\tilde{\mathbfcal{G}}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\mathbf{x}_{n},\mathbf{v}_{n})\right)\end{bmatrix}$
$\displaystyle-\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}-\mathbfcal{C}\right)\tilde{\mathbfcal{G}}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}(\hat{\mathbf{x}}_{n})\right)\end{bmatrix}$
$\displaystyle+\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}-\mathbfcal{C}\right)\tilde{\mathbfcal{G}}_{n+1}\mathbf{w}_{n+1}.$

### SM-3: Proof of Lemma 1

On one hand, note that

$\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}-\mathbfcal{C}=\left(\frac{1}{|\mathcal{N}|}\mathbf{1}\mathbf{1}^{\mathsf{T}}-\mathbf{C}\right)\otimes\mathbf{I}$

and on the other hand, one has

$\left(\frac{1}{|\mathcal{N}|}\mathbf{1}\mathbf{1}^{\mathsf{T}}-\mathbf{C}\right)^{n}=(-1)^{n-1}\left(\frac{1}{|\mathcal{N}|}\mathbf{1}\mathbf{1}^{\mathsf{T}}-\mathbf{C}\right)\mathbf{C}^{n-1}.$

Moreover, from Assumption 2, it follows that

$\forall\mathbf{a}:\exists a\hskip 5.69046pt\text{so that}\hskip
5.69046pt\lim_{n\rightarrow\infty}\mathbf{C}^{n}\mathbf{a}=a\mathbf{1}$

Therefore,

$\displaystyle\forall\mathbf{a}:\hskip 5.69046pt\lim_{n\rightarrow\infty}$
$\displaystyle\left(\frac{1}{|\mathcal{N}|}\mathbf{1}\mathbf{1}^{\mathsf{T}}-\mathbf{C}\right)^{n}\mathbf{a}$
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}(-1)^{n-1}\left(\frac{1}{|\mathcal{N}|}\mathbf{1}\mathbf{1}^{\mathsf{T}}-\mathbf{C}\right)\mathbf{C}^{n-1}\mathbf{a}$
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}(-1)^{n-1}\left(\frac{1}{|\mathcal{N}|}\mathbf{1}\mathbf{1}^{\mathsf{T}}-\mathbf{C}\right)(a\mathbf{1})=0$

indicating
$\mathsf{p}\left\\{\frac{1}{|\mathcal{N}|}\mathbf{1}\mathbf{1}^{\mathsf{T}}-\mathbf{C}\right\\}<1$,
and hence,
$\mathsf{p}\left\\{\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}-\mathbfcal{C}\right\\}<1$.
Thus, concluding the proof.

### SM-4: Proof of Theorem 1:

From Lemma 1 and given (6), it follows that $\|\mathrm{Res}_{\text{ob}}\|^{2}$
will be bounded. Moreover, conventionally, it is assumed that
$\|\mathrm{Res}_{\mathrm{f}}\|^{2}$,
$\|\mathrm{Res}_{\mathrm{h}_{i},\mathrm{f}}\|^{2}$, and
$\mathsf{tr}\left\\{\bm{\Sigma}_{\mathbf{w}}\right\\}$ are also bounded [30,
24, 25]. Moreover, from the conditions of this theorem and by definition,

$\forall i\in\mathcal{N}:\mathsf{p}\left\\{\sum_{\forall
j\in\mathcal{N}_{i}}\left(\mathbf{I}-\tilde{\mathbf{G}}_{j,n}\mathbf{H}_{j}\right)\mathbf{A}\right\\}<1.$

Thus, given Assumption 2, one has

$\mathsf{p}\left\\{\mathbfcal{C}(\mathbf{I}-\tilde{\mathbfcal{G}}_{n}\mathbfcal{H})\mathbfcal{A}\right\\}<1$

concluding the proof.

### SM-5: Proof of Theorem 2:

Substituting (8) into (12) gives

$\displaystyle\bm{\Delta}_{n+1}=$
$\displaystyle\mathbfcal{E}_{n+1}-\mathbf{1}\otimes\bm{\epsilon}_{n+1}$
$\displaystyle=$ $\displaystyle\begin{bmatrix}\mathbf{x}_{n+1}-\sum_{\forall
j\in\mathcal{N}_{1}}\bm{\phi}_{j,n+1}\\\ \vdots\\\
\mathbf{x}_{n+1}-\sum_{\forall
j\in\mathcal{N}_{|N|}}\bm{\phi}_{j,n+1}\end{bmatrix}-\begin{bmatrix}\mathbf{x}_{n+1}-\hat{\mathbf{x}}_{n+1}\\\
\vdots\\\ \mathbf{x}_{n+1}-\hat{\mathbf{x}}_{n+1}\end{bmatrix}$
$\displaystyle=$
$\displaystyle\mathbfcal{C}\begin{bmatrix}\hat{\mathbf{x}}_{n+1}-\bm{\phi}_{1,n+1}\\\
\vdots\\\ \hat{\mathbf{x}}_{n+1}-\bm{\phi}_{|\mathcal{N}|,n+1}\end{bmatrix}$

which yields

$\|\bm{\Delta}_{n+1}\|^{2}=\mathsf{p}\left\\{\mathbfcal{C}\right\\}\sum_{\forall
i\in\mathcal{N}}\|\hat{\mathbf{x}}_{n+1}-\bm{\phi}_{i,n+1}\|^{2}.$ (24)

Substituting (17) into (24) gives

$\|\bm{\Delta}_{n+1}\|^{2}\leq\mathsf{p}\left\\{\mathbfcal{C}\right\\}\gamma\sum_{\forall
i\in\mathcal{N}}\|\hat{\mathbf{x}}_{n+1}-\bm{\phi}_{i,n+1}\|^{2}$

with $\gamma=\max\\{\gamma_{i},\forall i\in\mathcal{N}\\}$, and thus, it would
transpire that

$\|\bm{\Delta}_{n+1}\|^{2}\leq\mathsf{p}\left\\{\mathbfcal{C}\right\\}\gamma\|\bm{\Delta}_{n}\|^{2}.$
(25)

Given that $\mathsf{p}\left\\{\mathbfcal{C}\right\\}=1$ and $\gamma\in(0,1)$,
(25) indicates all agents will converge and $\|\bm{\Delta}_{n}\|^{2}$ becomes
bounded. Hence, concluding the proof.

### SM-6: Proof of Theorem 3:

From (19) and Assumption 4 one has

$\displaystyle\bm{\Delta}_{n+1}=$
$\displaystyle\mathbfcal{C}\begin{bmatrix}\mathrm{f}(\hat{\mathbf{x}}_{n})-\mathrm{f}(\hat{\mathbf{x}}_{1,n})\\\
\vdots\\\
\mathrm{f}(\hat{\mathbf{x}}_{n})-\mathrm{f}(\hat{\mathbf{x}}_{|\mathcal{N}|,n})\end{bmatrix}$
(26)
$\displaystyle+\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}\mathbfcal{G}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}\left(\mathbf{x}_{n},\mathbf{v}_{n}\right)\right)-\mathrm{h}_{1}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{n}\right)\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}\left(\mathbf{x}_{n},\mathbf{v}_{n}\right)\right)-\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{n}\right)\right)\end{bmatrix}$
$\displaystyle-\mathbfcal{C}\bm{\Upsilon}\tilde{\mathbfcal{G}}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}\left(\mathbf{x}_{n},\mathbf{v}_{n}\right)\right)-\mathrm{h}_{1}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{1,n}\right)\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}\left(\mathbf{x}_{n},\mathbf{v}_{n}\right)\right)-\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{|\mathcal{N}|,n}\right)\right)\end{bmatrix}$

where
$\bm{\Upsilon}=\mathsf{diag}\left\\{\bm{\Upsilon}_{i}:i=1,\ldots,|\mathcal{N}|\right\\}$.
The expression in (26) can be rearranged to give

$\displaystyle\bm{\Delta}_{n+1}=$
$\displaystyle\mathbfcal{C}\begin{bmatrix}\mathrm{f}(\hat{\mathbf{x}}_{n})-\mathrm{f}(\hat{\mathbf{x}}_{1,n})\\\
\vdots\\\
\mathrm{f}(\hat{\mathbf{x}}_{n})-\mathrm{f}(\hat{\mathbf{x}}_{|\mathcal{N}|,n})\end{bmatrix}+\mathrm{Res}_{\bm{\Delta}_{m}}$
(27)
$\displaystyle-\mathbfcal{C}\bm{\Upsilon}\tilde{\mathbfcal{G}}_{n+1}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{n}\right)\right)-\mathrm{h}_{1}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{1,n}\right)\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{n}\right)\right)-\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{|\mathcal{N}|,n}\right)\right)\end{bmatrix}$

with
$\mathbf{S}=\left(\frac{1}{|\mathcal{N}|}\textbf{\emph{1}}\mathbfcal{G}_{n+1}-\mathbfcal{C}\bm{\Upsilon}\tilde{\mathbfcal{G}}_{n+1}\right)$
so that

$\mathrm{Res}_{\bm{\Delta}_{n}}=\mathbf{S}\begin{bmatrix}\mathrm{h}_{1}\left(\mathrm{f}\left(\mathbf{x}_{n},\mathbf{v}_{n}\right)\right)-\mathrm{h}_{1}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{n}\right)\right)\\\
\vdots\\\
\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}\left(\mathbf{x}_{n},\mathbf{v}_{n}\right)\right)-\mathrm{h}_{|\mathcal{N}|}\left(\mathrm{f}\left(\hat{\mathbf{x}}_{n}\right)\right)\end{bmatrix}.$
(28)

Once more, consider $\exists m\in\mathbb{N}$ so that
$\|\hat{\mathbf{x}}_{m}-\hat{\mathbf{x}}_{i,m}\|^{2}\leq\mu^{\prime}$
resulting in the approximations in (14) to hold, and hence, yielding

$\mathsf{E}\left\\{\|\bm{\Delta}_{m+1}\|^{2}\right\\}\leq
g_{\bm{\Delta}_{m}}\mathsf{E}\left\\{\|\bm{\Delta}_{m}\|^{2}\right\\}+\|\mathrm{Res}_{\bm{\Delta}_{m}}\|^{2}.$
(29)

Through $k$ iteration, (29) gives

$\displaystyle\mathsf{E}\left\\{\|\bm{\Delta}_{m+1}\|^{2}\right\\}\leq$
$\displaystyle\left(\prod^{k}_{i=0}g_{\bm{\Delta}_{m-i}}\right)\mathsf{E}\left\\{\|\bm{\Delta}_{m-k}\|^{2}\right\\}$
(30)
$\displaystyle+\sum^{k}_{i=0}\left(\prod^{i-1}_{j=0}g_{\bm{\Delta}_{m-j}}\right)\|\mathrm{Res}_{\bm{\Delta}_{m-k}}\|^{2}$
$\displaystyle+\|\mathrm{Res}_{\bm{\Delta}_{m}}\|^{2}.$

From (28), it follows that $\|\mathrm{Res}_{\bm{\Delta}_{m}}\|^{2}$ is
bounded. Thus, from Assumption 2 and 4 it follows that the expression in (30)
converges, if

$\exists k\in\mathbb{N}:\left(\prod^{k}_{i=0}g_{\bm{\Delta}_{m-i}}\right)<1$

which concludes the proof.

### SM-7: Network used in Simulation Examples

Fig. 3 shows typology of the network used in all simulation examples. In
keeping with the literature conversion [7, 9, 14] all connections are assumed
bidirectional.

Fig. 3: Network of $120$ nodes and $352$ connections used in simulation
examples.


print-report

Pseudo Complexity of Purification for free scalar field theories

Aranya<EMAIL_ADDRESS>Arpan
<EMAIL_ADDRESS>Sabyasachi
<EMAIL_ADDRESS>

${}^{a}\,$Centre for High Energy Physics, Indian Institute of Science, C.V.
Raman Avenue,
Bangalore 560012, India. ${}^{b}\,$Indian Institute of Technology,
Gandhinagar, Gujarat 382355, India. ${}^{c}\,$Theory Division, Saha Institute
of Nuclear Physics, 1/AF, Bidhannagar 700064,
West Bengal, India. ${}^{d}\,$Homi Bhabha National Institute, Training School
Complex, Anushakti Nagar,
Mumbai 400094, India.

###### Abstract

We compute the _pseudo complexity of purification_ corresponding to the
reduced transition matrices for free scalar field theories with an arbitrary
dynamical exponent. We plot the behaviour of complexity with various
parameters of the theory under study and compare it with the complexity of
purification of the reduced density matrices of the two states
$|\psi_{1}\rangle$ and $|\psi_{2}\rangle$ that constitute the transition
matrix. We first find the transition matrix by reducing to a small number ($1$
and $2$) of degrees of freedom in lattice from a lattice system with many
lattice points and then purify it by doubling the degrees of freedom ($2$ and
$4$ respectively) for this reduced system. This is a primary step towards the
natural extension to the idea of the complexity of purification for reduced
density matrices relevant for the studies related to postselection.

##  1 Introduction

Quantum entanglement continues to grab eyeballs across multiple disciplines.
Recently, a generalization of entanglement entropy: called pseudo entropy was
proposed in [1, 2]. There, instead of starting from a usual density matrix,
the authors define a matrix using two pure quantum states $\ket{\psi_{1}}$ and
$\ket{\psi_{2}}$

$\tau^{1|2}=\frac{\ket{\psi_{1}}\bra{\psi_{2}}}{\langle\psi_{2}|\psi_{1}\rangle},$
(1.1)

which is dubbed the transition matrix. This has a two-fold interest to begin
with. The first is that in the case of post-selection in quantum experiments,
the transition matrix becomes important once one specifies the initial state
as $|\psi_{1}\rangle$ and the final state as $|\psi_{2}\rangle$, given they
are not orthogonal to each other. In the post-selection experiment, the
transition matrix then plays the role of the density matrix while computing
the weak expectation value
$\langle\mathcal{O}\rangle=\text{Tr}(\mathcal{O}\tau^{1|2})$ of an observable
$\mathcal{O}$. The second motivation comes from the AdS/CFT correspondence [3,
4], where pseudo entropy is proposed to be the dual CFT quantity to the area
of a minimal hypersurface in Euclidean AdS spacetime. In holography, this is
constructed in the path integral technique by dividing a Euclidean timeslice
into two regions corresponding to the two states. Dividing the entire system
into $A$ and its complement $A^{C}$, the pseudo entropy of $A$ (relative to
$A^{C}$) is defined as

$S\left(\tau^{1|2}\right)=-\text{Tr}_{A}\left[\tau_{A}^{1|2}\log\tau_{A}^{1|2}\right],$
(1.2)

where $\tau_{A}^{1|2}=\text{Tr}_{A^{C}}\left(\tau^{1|2}\right)$ is called the
reduced transition matrix, in analogy with existing literature on entanglement
entropy derived from reduced density matrix. Once the two states are taken to
be the same, the transition matrix reduces to the usual density matrix and
pseudo entropy reduces to usual entanglement entropy.

Another quantum information (QI) theoretic quantity that has been at the
centre of interest for the last few years is quantum circuit complexity. In
recent times, circuit complexity has been explored in the context of quantum
field theory [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23, 24, 25] 111This list is by no means exhaustive. Interested readers are
referred to the reviews [26, 27], and references therein for more details..
Complexity is usually defined as the number of elementary structural
components needed to construct either an evolved state or operator from a
simple initial state or operator. Here we will try to compute the circuit
complexity of the reduced transition matrix. This reduced transition matrix is
an analogue of a mixed state corresponding to a reduced density matrix,
although it can be non-hermitian in general. The methods for computing circuit
complexity we apply are those for the complexity of purification [28, 29, 30,
31, 32, 33, 34].

First, given a transition matrix, which is written in density matrix form, a
state can be associated with this operator by writing the transition matrix in
the vectorized form. Then we compute the circuit complexity for this state
using Nielsen’s method [35, 36, 37]. This vectorization means changing the bra
part of the transition matrix to a ket which amounts to a doubling of the
Hilbert space[38, 39, 40]. Hence the corresponding transition state looks like

$\psi_{1|2}=\frac{|\psi_{1}\rangle|\psi_{2}\rangle}{\langle\psi_{1}|\psi_{2}\rangle},$
(1.3)

which, as mentioned in [2], for a transition matrix, turns the corresponding
state into a tensor product of the two states. Then the corresponding
complexity, following the approach of [5], of the state dual to the transition
matrix being

$C(\tau^{1|2})=C(|\psi_{1}\rangle)\oplus C(|\psi_{2}\rangle).$ (1.4)

Note that, the complexity defined using canonical purification in (1.4) is
simply the direct sum of two complexities of two individual states
$|\psi_{1,\,2}\rangle$. The direct sum is representative of the fact that the
vectorized version is a tensor product. It is quite unlikely that this will
have any nontrivial holographic interpretation, as typically holographic duals
involve some sort of minimization. Motivated by this, we want to define some
notion of pseudo complexity of purification ($C(\tau_{A}^{1|2})$) in this
paper extending the notion of complexity of purification [29, 30, 31, 32, 33,
34] for this case. We first purify the reduced transition matrix and compute
the minimal complexity among all possible purifications 222The state in (1.3),
can be obtained by using “operator-state mapping,” a particular example of
purification, namely, “canonical purification”.. We do these computations for
free scalar field theories and Lifshitz field theories in $\left(1+1\right)$
dimensions.

One can ask what pseudo complexity means physically. This is related to the
idea of post-selection as suggested in [1, 2, 41] while defining pseudo
entropy. The idea is that if, along some physical process, one starts with a
state $|\psi_{1}\rangle$ and somehow the final state is post-selected to be
$|\psi_{2}\rangle$, which is _not_ necessarily a result of a simple unitary
evolution with a simple Hamiltonian, the process can be still approximated as
a result of an operator $|\psi_{2}\rangle\langle\psi_{1}|$ operating on the
initial state. In general, $|\psi_{2}\rangle$ might be a result of many
unitary and non-unitary steps along the way 333By non-unitary, we mean that
there might be measurements made along an evolution due to which the state
might collapse to a different state which results in an overall non-unitarity.
It is worth noting that the process of purification by itself is a non-unitary
one, which however is different from the non-unitarity caused by the post-
selection.. Hence, this operator $|\psi_{2}\rangle\langle\psi_{1}|$ might as
very well be a non-hermitian operator. This operator is considered a
transition matrix, similar to a density matrix in this treatment. The reason
is that this operator does not have any nontrivial information of the system
other than the outer product of the initial and the post-selected state,
similar to a density matrix. Then the reduced transition matrix, defined after
a reduction in degrees of freedom, can associate a notion of entropy similar
to the entanglement entropy of a reduced density matrix. Similarly, the pseudo
complexity is supposed to measure the complexity of this operator, whereas the
pseudo complexity of purification is supposed to measure the complexity of
forming such a reduced transition matrix once the two states
$|\psi_{1,2}\rangle$ are given. Another point of view of looking at this
problem is to relate it to the averaged number of maximally entangled pair of
qubits to be distilled from the intermediate state once the final state is
postselected [1]. While pseudo entropy measures this number, pseudo complexity
of purification can be thought to be measuring the amount of work (of course,
in terms of resources or gates) needed for this distillation starting from the
reduced transition matrix.

The rest of this article is organized as follows: in Sec. (2) we first review
a path-integral calculation of Gaussian transition matrix following [41] and
propose a definition of pseudo-complexity of purification. We then use our new
definition and calculate this quantity in a simple setting of coupled harmonic
oscillators on a lattice. The results of our numerical analyses are collected
in Sec. (3). We conclude in Sec. (4) with discussions and outlook.

##  2 Pseudo Complexity of Purification: A Proposal

The concept of purification originates from quantum information theory (can be
considered to be a close cousin of Schmidt decomposition as a process), and
the process of purification refers to constructing a pure state from a mixed
reduced density (/transition) matrix such that if one traces out the auxiliary
degrees of freedom added while constructing the purified state, one gets back
the original reduced density (/transition) matrix. While both the actual and
purified state generates the same reduced density matrix, the purified state
only has important information about the density matrix because the extra
parameters added in the process of purification are random. The only
constraints on them are that the purified state should follow the standard
quantum mechanical postulates and properties of a pure state. As mentioned
already, we will be considering random purification by all possible values of
parameters for which the purified state can be consistently formed. However,
we will choose a specific purification among them by minimizing the complexity
functional [30]. This is motivated by the definition of complexity, where
minimal resources play an important part. In terms of quantum gates, this
means choosing the minimum number of gates so that resources needed are
minimized. Now we describe our main setup.

###  2.1 Setup

We consider a system of scalar fields which is governed by the Hamiltonian

$H=\frac{1}{2}\int
dx\left[\pi^{2}+\left(\partial_{x}^{z}\right)^{2}+m^{2z}\phi^{2}\right].$
(2.1)

This theory exhibits an anisotropic (Lifshitz) scaling symmetry
$t\to\lambda^{z}t$, $x\to\lambda x$ in the $m\to 0$ limit; for $z=1$, this is
an ordinary relativistic scalar field theory. For practial purposes, it is
useful to consider the discretized Hamiltonian

$H=\sum_{i=1}^{N}\left[\frac{\pi_{i}^{2}}{2}+\frac{m^{2z}}{2}\phi_{n}^{2}+\frac{1}{2}\left(\sum_{k=0}^{z}\left(-1\right)^{z+k}\begin{pmatrix}z\\\
k\end{pmatrix}\phi_{i-1+k}\right)^{2}\right],$ (2.2)

N being the total number of points in the lattice.

We consider a collection of coupled linear harmonic oscillators on a lattice
of space points, labelled by capital Latin indices. Any two Gaussian states
$\ket{\psi_{\alpha}}$ $\left(\alpha=1,2\right)$ of this system can be
generically expressed in their position representation as

$\langle
q_{A}|\psi_{\alpha}\rangle=N^{\left(\alpha\right)}\exp\left[-\frac{1}{2}q_{A}W_{AB}^{\left(\alpha\right)}q_{B}\right],$
(2.3)

where $q_{A}$ denotes displacement of the $A$-th oscillator,
$W_{AB}^{\left(\alpha\right)}$ is a positive definite matrix and $N^{\alpha}$
is the normalization constant

$N^{\left(\alpha\right)}=\left(\text{det}\left(\frac{W_{AB}^{\left(\alpha\right)}}{\pi}\right)\right)^{\frac{1}{4}}.$

We divide the entire system into a part $A$ and its complement $A^{C}$, the
lattice points within the sub-system $\Omega$ are labelled by lowercase Latin
letters, while lowercase Greek indices label those outside it. We adopt the
notation.

$W_{AB}^{\left(\alpha\right)}=\begin{pmatrix}W_{ab}^{\left(\alpha\right)}&W_{a\beta}^{\left(\alpha\right)}\\\
W_{\alpha
b}^{\left(\alpha\right)}&W_{\alpha\beta}^{\left(\alpha\right)}\end{pmatrix}=\begin{pmatrix}A^{\left(\alpha\right)}&B^{\left(\alpha\right)}\\\
B^{\left(\alpha\right)T}&C^{\left(\alpha\right)}\end{pmatrix},$ (2.4)

to denote sub-matrices of $W_{AB}^{\left(\alpha\right)}$. The reduced
transition matrix for the sub-system $\Omega$ in this representation is found
by integrating over the rest

$\langle
q_{a}|\text{Tr}_{A^{C}}\left(\ket{\psi_{1}}\bra{\psi_{2}}\right)|q_{b}\rangle=\int
dq_{\alpha}\langle
q_{a},q_{\alpha}|\psi_{1}\rangle\langle\psi_{2}|q_{b},q_{\alpha}\rangle.$
(2.5)

This is a Gaussian integral. Finally, the matrix elements of the reduced
transition matrix can be expressed in a simple form [2]

$\langle
q^{\left(1\right)}|\tau_{A}^{1|2}|q^{\left(2\right)}\rangle=\frac{N^{\prime}}{\sqrt{\text{det}\frac{\bar{C}}{\pi}}}\times\exp\left[-\frac{1}{2}\begin{pmatrix}q^{\left(1\right)T}&q^{\left(2\right)T}\end{pmatrix}M\begin{pmatrix}q^{\left(1\right)}\\\
q^{\left(2\right)}\end{pmatrix}\right]\,,$ (2.6)

where

$M=\begin{pmatrix}X^{\left(1\right)}&2Y\\\
2Y^{T}&X^{\left(2\right)}\end{pmatrix}$ (2.7)

and

$\displaystyle
X^{\left(\alpha\right)}=A^{\left(\alpha\right)}-\frac{1}{2}B^{\left(\alpha\right)}\bar{C}^{-1}B^{\left(\alpha\right)T},$
$\displaystyle
Y=-\frac{1}{4}B^{\left(1\right)}\bar{C}^{-1}B^{\left(2\right)T},$
$\displaystyle\bar{C}=\frac{1}{2}\left(C^{\left(1\right)}+C^{\left(2\right)}\right),$
$\displaystyle N^{\prime}=\sqrt{\text{det}\frac{\bar{W}}{\pi}},$
$\displaystyle\bar{W}=\frac{1}{2}\left(W^{\left(1\right)}+W^{\left(2\right)}\right)\,.$
(2.8)

###  2.2 Auxiliary parameters and purification:

We use Nielsen’s method for the computation of circuit complexity [35, 36, 37,
5] and choose the $F^{1}$ norm for the complexity functional [5, 42]. For a
pure (or purified) state, this is written in terms of the (unentangled)
reference state frequency and the frequency of the normalized version of the
coupled oscillator system in the lattice pictures. We start from a various
number of initial oscillators and reduce to transition matrices after tracing
out all but one or two of the oscillators. Afterwards, we purify this reduced
transition matrix $\tau_{A}^{1|2}$ by adding parameters corresponding to
adding one or two more oscillators (similar to the doubling of degrees of
freedom in qubits purifications). However, we do not choose canonical
purification since it is always returns the original transition state defined
in (1.3), the complexity of which is trivial, as explained before 444As
mentioned in (1.4), for canonical purification the complexity of the purified
state is simply the sum of the complexity of $|\psi_{1}\rangle$ and
$|\psi_{2}\rangle$.. For the reduced transition matrix, the parameters that
come into the picture while matching the original reduced transition matrix
with the reduced density matrix of the purified state take care of all
possible purifications. In this optimal picture of purification, the number of
unknown arbitrary parameters is much higher than the number of unknowns needed
to purify $\rho_{1}^{A}$ or $\rho_{2}^{A}$. This is because the information of
not one but two pure states goes into the transition matrix and hence the
reduced transition matrix.

In order to define pseudo-complexity, we shall begin by introducing
purification of the transition matrix $\tau_{A}^{1|2}$. We consider a
(fictional) auxiliary system $\tilde{A}$ whose Hilbert space is of the same
dimension as our original subsystem $A$ and consider two Gaussian states in
the enlarged Hilbert space $\mathcal{H}_{A}\otimes\mathcal{H}_{\tilde{A}}$

$\langle
q_{A},q_{\tilde{A}}|\Psi_{\alpha}\rangle=N_{A\tilde{A}}\exp\left[-\frac{1}{2}\begin{pmatrix}q_{A}^{T}&q_{\tilde{A}}^{T}\end{pmatrix}\begin{pmatrix}J^{\left(\alpha\right)}&K^{\left(\alpha\right)}\\\
K^{\left(\alpha\right)T}&L^{\left(\alpha\right)}\end{pmatrix}\begin{pmatrix}q_{A}\\\
q_{\tilde{A}}\end{pmatrix}\right],\hskip 14.22636pt\left(\alpha=1,2\right)$
(2.9)

then proceeding exactly as above, we calculate a transition matrix in the
enlarged Hilbert space. This particular choice of Gaussian purification is
primarily motivated by the study done in [43] in the context of entanglement
of purification. This simplifies the numerical analysis as well as the
computation of circuit complexity. We consider that this is a purification of
our original reduced transition matrix (2.6) for the subsystem $A$; in this
case, the following constraints must be obeyed.

$\displaystyle
J^{\left(\alpha\right)}-\frac{1}{2}K^{\left(\alpha\right)}\bar{L}^{-1}K^{\left(\alpha\right)T}=A^{\left(\alpha\right)}-\frac{1}{2}B^{\left(\alpha\right)}\bar{C}^{-1}B^{\left(\alpha\right)T},$
(2.10) $\displaystyle
K^{\left(1\right)}\bar{L}^{-1}K^{\left(2\right)T}=B^{\left(1\right)}\bar{C}^{-1}B^{\left(2\right)T},$
(2.11)

where $\bar{L}=\frac{1}{2}\left(L^{\left(1\right)}+L^{\left(2\right)}\right)$.
In our numerical analysis, we adopt the following conditions to satisfy these
constraints

$\displaystyle J^{\left(\alpha\right)}=A^{\left(\alpha\right)},$ (2.12a)
$\displaystyle\bar{L}=\left(K^{\left(2\right)-1}B^{\left(2\right)}\bar{C}^{-1}\left(K^{\left(2\right)}B^{\left(2\right)}\right)^{T}\right)^{-1},$
(2.12b) $\displaystyle
K^{\left(2\right)}=\left(\bar{L}K^{\left(1\right)-1}B^{\left(1\right)}\bar{C}^{-1}B^{\left(2\right)T}\right)^{T}\,.$
(2.12c)

This particular choice leaves only $K^{\left(1\right)}$ and
$L^{\left(1\right)}$ undetermined.

For the particular model of our choice, the matrices
$W_{AB}^{\left(\alpha\right)}$ are given by

$W_{AB}^{\left(\alpha\right)}=\frac{1}{N}\sum_{C=1}^{N}\sqrt{m_{\alpha}^{2z_{\alpha}}+\left(2\left(1-\cos\frac{2\pi
n}{N}\right)\right)^{z_{\alpha}}}\exp\left(\frac{2\pi
iC\left(A-B\right)}{N}\right),\hskip 14.22636pt\left(\alpha=1,2\right),$
(2.13)

where $N$ is the total number of lattice points. $m_{\alpha}$ and $z_{\alpha}$
are the mass and dynamical exponent of the respective theory. The matching of
$J,\,K,\,L^{(\alpha)}$ matrices with their corresponding counterparts
introduce a few arbitrary parameters in the purified state. These parameters
denote the infinitely many possible Gaussian purifications of the reduced
transition matrix. Afterwards, we define the complexity of purification as the
complexity of the specific purification for which the complexity functional is
minimized in terms of all the parameters of the purified version of the
reduced transition matrix. This involves purification parameters for both
$|\psi_{1,\,2}\rangle$ states, constituting the transition matrix. Finally, we
compare the complexity of purification of the reduced transition matrix
$C(\tau_{A}^{(1|2)})$ with the complexity of purifications for the individual
reduced density matrices ($\rho_{A}^{1}$ and $\rho_{A}^{2}$) derived from the
initial and final states $|\psi_{1}\rangle$ and $|\psi_{2}\rangle$
respectively.

We propose that the pseudo complexity of purification is given by

$\mathcal{C}_{P}=\underset{\tilde{A}}{\mathrm{min}}\frac{1}{2}\sqrt{\sum_{j=1}^{N}\sum_{\alpha=1}^{2}\log\left(\frac{\Omega^{\left(\alpha\right)}_{j}}{\omega_{2}}\right)},$
(2.14)

where $\Omega^{\left(\alpha\right)}_{j}$ happens to denote the $j$-th
eigenvalue of the matrix
$\begin{pmatrix}J^{\left(\alpha\right)}&K^{\left(\alpha\right)}\\\
K^{\left(\alpha\right)T}&L^{\left(\alpha\right)}\end{pmatrix}$, the functional
is minimized over the parameters associated with the auxiliary system
$\tilde{A}$. This means choosing a particular purification associated with
minimizing parameter values.

Let us point out a few attributes of our method of computation. Most of our
numerical analyses are performed for the simplest case of a reduced system
made of only one oscillator; where the above functional is to be minimized
over two unknown parameters. In general, when the reduced system consists of
$n$ no. of linear harmonic oscillators; the $J$, $K$, and $L$ matrices are of
dimension $n\times n$. The choice (2.12) always leaves two of them
undetermined, and therefore we end up with a total of $2n^{2}$ unknown
parameters. On the other hand, to calculate the more well-known complexity of
purification from a reduced density matrix, the no. of unknown parameter
values over which one is required to perform an minimization scales as $n^{2}$
with the reduced system size. It is easy to conclude that the calculation of
the new functional (2.14) is at least twice as hard since the transition
matrices have information of both the initial as well as the final state.

There, however, exists an way to purify the system mode-by-mode [30]. That way
it might be possible to reduce the no. of unknown parameters to $2n$. Even
then, the problem is harder (twice) than computing complexity of purification
from density matrices, where the associated no. scales as $n$.

It may also be of relevance to note that any purification we consider takes
into account only Gaussian states and the auxiliary system always has the same
dimension as the original one, motivated by the usual doubling of Hilbert
space typically done in purification using Schmidt decomposition. These are
restrictions we invoke in order to keep the numerical job simple. There may
exist ways to relax the conditions.

##  3 Results and discussion

###  3.1 On the pseudo complexity of purification:

In this subsection, we report our numerical results after the minimization
with varying different variables of the theory. In the primary Hamiltonian
(2.2), we chose a general version which applies to both the ordinary free
scalar QFTs $(z_{1}=z_{2}=1)$ as well as Lifshitz scalar field theories
$(z_{1}\neq 1\neq z_{2})$. The variables considered are the dynamical
exponents $z_{1,2}$, the masses $m_{1,2}$, the number of oscillators before
tracing out and the number of oscillators of the reduced system and the
frequency of the unentangled reference state.

  1. 1.

Varying $z_{2}$ and $m_{2}$: From the plots of pseudo complexity of
purification against reference frequency in Figs. (1(a) )and (1(b)), we notice
that as the reference frequency is increased, the pseudo complexity of
purification decreases primarily and then shows a polynomial growth and
saturation. When the scaling factor $z_{2}$ is increased in the same plots,
the pseudo complexity of purification for higher values of reference frequency
decreases for higher $z_{2}$. However, the behaviour is opposite for small
reference frequency values. Hence, there is a crossover in complexity after
which the small to large behaviour changes. Similar change is observed from
the Fig. (1(b)) where the mass parameter $m_{2}$ is varied. We find that at
large frequencies, the pseudo complexity of purification decreases as the mass
is increased, whereas, at very small frequencies, the behaviour is again
opposite.

(a)

(b)

Figure 1: Behaviour of $\mathcal{C}_{P}$ with increasing $\omega_{2}$ for (a)
$m_{1}=m_{2}=0.5,z_{1}=1$ and different $z_{2}$ and (b)
$m_{1}=1.0,z_{1}=1,z_{2}=2$ and different $m_{2}$.

  2. 2.

Pseudo complexity of purification vs $z_{2}$ plots: From Fig. (2(a)) we
observe that as $z_{2}$ is increased, the pseudo complexity of purification
primarily goes down and then grows linearly. By looking at multiple plots for
various values of $z_{1}$, we find that with increasing $z_{1}$ as $z_{2}$ is
varied, the pseudo complexity of purification decreases. In this case, there
is no change of behaviour or crossover between different plots, which we saw
earlier.

  3. 3.

Pseudo complexity of purification vs $m_{2}$ plots: As shown in Fig. (2(b)),
we find that one of the mass parameters is varied from zero to one, and the
pseudo complexity of purification does not change much. Although the plots are
not completely parallel to the x-axis, their variance is relatively much
smaller than what we have observed for other parameters. On the other hand, if
we vary one of the scaling parameters $z_{2}$, we find that with increasing
values of $z_{2}$, pseudo complexity of purification increases.

(a)

(b)

Figure 2: Behaviour of $\mathcal{C}_{P}$ with increasing (a) $z_{2}$ and (b)
$m_{2}$ for $\omega_{2}=25$.

  4. 4.

Varying number of oscillators before tracing out: All the above plots are
given for the case where the total number of oscillators associated with the
states $\psi_{1,2}$ before tracing out degrees of freedom is $2$. So we trace
out one of the two oscillators in each of the states. This begs the question
of how the results might possibly change in case the number of oscillators are
increased, which is more relevant when thinking about quantum field theories,
which can be written as the lattice of infinitely many coupled harmonic
oscillators. We, however, find that the pseudo complexity of purification does
not vary much at all when the number of oscillators of the actual system is
varied. This can be seen from Fig. (3).

Figure 3: Comparison of $\mathcal{C}_{P}$ for different total number of
lattice sites $\left(N\right)$. In all graphs $m_{1}=m_{2}=0.5$ and $z_{1}=1$,
$z_{2}=2$.

  5. 5.

Varying number of oscillators in the reduced system: All of our results so far
are for the case when the number of oscillators in the reduced transition
matrix is one. The reason is two-fold: first and foremost, as we have
mentioned before, the number of unknown parameters for pseudo complexity of
purification in terms of which the complexity functional is minimized is much
more than the case of usual complexity of purification. This makes the problem
much more challenging numerically as the number of parameters grow very
quickly as one increases the number of oscillators in the reduced system.
Secondly, from a physical point of view, with the motivation of relating the
entanglement entropy to some thermal notion of entropy, one typically chooses
the number of oscillators in the reduced system to be small. Although there is
no such clear motivation for complexity, we stick to this particular choice.

The qualitative behaviour of the pseudo complexity of purification with the
varying reference frequency remains similar when the subsystem size is
increased to include two oscillators as shown in Fig. (4). Hence, we expect
all of our results to be universal and hold for larger subsystem size.

Figure 4: Pseudo complexity of purification $\mathcal{C}_{P}$ when the reduced
system contains two LHOs.

###  3.2 On the difference $\Delta C$: Mutual pseudo complexity of
purification

The reduced transition matrix carries complicated mixed information of not one
but two pure states $\psi_{1}$ and $\psi_{2}$. Hence, while studying the
complexity of purification for the reduced transition matrix, it is natural to
ask how different it is from the complexities of individual reduced density
matrices of those two states. Inspired by the definition of mutual complexity
$\left(\triangle C\left(\rho_{AB}\right)\right)$ [44, 45, 30], which is
defined as the difference of complexity between the full state $\rho_{AB}$ and
the sum of complexities of the two reduced density matrices $\rho_{A}$ and
$\rho_{B}$, we define the difference as the _mutual pseudo complexity of
purification_. In the following, we report some properties of the mutual
pseudo complexity of purification based on our numerical analysis.

  1. 1.

Subadditivity: From figure 5, it is easy to observe that

$\triangle
C\left(\tau_{A}^{1|2}\right)=C\left(\tau_{A}^{1|2}\right)-\frac{C\left(\rho_{A}^{1}\right)+C\left(\rho_{A}^{2}\right)}{2}\leq
0.$ (3.1)

Therefore, we find that the pseudo complexity of purification is always sub-
additive to the sum of the individual complexities of purification of the two
relevant states. Our numerical results suggest that the term $\triangle
C\left(\tau_{A}^{1|2}\right)$ is always negative or zero. The exact equality
is expected for the case when the two states become the same and hence
$\rho_{A}^{1}=\rho_{A}^{2}$.

(a)

(b)

Figure 5: Comparison of $\mathcal{C}_{P}$ with individual complexities of
purification of the two relevant states $\psi_{1}$ and $\psi_{2}$, we let
$z_{1}=1$, $z_{2}=2$ and $m_{1}=m_{2}=0.5$.

  2. 2.

Near saturation behaviour: It is observed that both the complexity of
purification and the pseudo complexity of purification tend to saturate at a
very large reference frequency. By looking at the saturation values of the
corresponding plots in Fig. (5), we can make the following comment on the
saturation of pseudo complexity of purification in terms of complexities of
purification of the individual reduced density matrices,

$C(\tau_{A}^{1|2})\sim\frac{C(\rho_{A}^{1})+C(\rho_{A}^{2})}{2}$ (3.2)

This means that the quantity $\Delta C(\tau_{A}^{1|2})$ approaches zero at
large reference frequencies. This indicates that the dependence of
$C(\tau_{A}^{1|2})$ on the reference frequency $\omega$ is such that for large
enough values of $\omega$, the complexity does not distinguish between a
transition matrix and the density matrix.

  3. 3.

Varying masses: We also study the behaviour of
$\frac{C(\rho_{A}^{1})+C(\rho_{A}^{2})}{2}-C(\tau_{A}^{1|2})$ ($=\triangle
C(\tau_{A}^{1|2})$) with change in one of the masses $m_{2}$ corresponding to
the state $\psi_{2}$ while keeping the other parameters
($m_{1},\,z_{1},\,z_{2}$) fixed. We look at the behaviour both in a small mass
difference as well as a large mass difference regime. In both cases, we find
that the difference decreases as one increases the difference between mass
parameters, see Fig. (6). In the small mass difference regime, the slope of
the plot is linear. This means $\left|\triangle
C(\tau_{A}^{1|2})\right|\propto-a\left(m_{2}-m_{1}\right)$ in this range. On
the other hand in the large mass difference regime, the decay of the slope is
close to exponential inspiring us to write $\left|\triangle
C(\tau_{A}^{1|2})\right|\propto e^{-a\left(m_{2}-m_{1}\right)}$, with lower
saturation close to zero as $\left(m_{2}-m_{1}\right)\sim 1$.

(a)

(b)

Figure 6: Behaviour of $\frac{1}{2}\left(C_{1}+C_{2}\right)-\mathcal{C}_{P}$
(a) small mass difference and (b) large mass difference. In all results,
$\omega_{2}=50$.

  4. 4.

Varying $z_{2}$: In the same figure, we see that as the difference between
$z_{2}-z_{1}$ is increased, the absolute value of $\triangle
C(\tau_{A}^{1|2})$ increases. However, the lower saturation $\left|\triangle
C\left(\tau_{A}^{1|2}\right)\right|\sim 0$, near $(m_{2}-m_{1})\sim 1$,
remains unchanged. Hence, the numerics definitely suggests that the saturation
is universal and independent of $z_{2}-z_{1}$.

  5. 5.

$\left|\triangle C\right|\,\text{vs}\,z_{2}$ plot: As shown in Fig. (7), if we
keep $z_{1},\,m_{1},\,\text{and}\,m_{2}$ fixed while varying $z_{2}$,
$|\triangle C|$ increases linearly.

$|\triangle C|\propto\alpha z_{1}+\beta z_{2},$ (3.3)

the above expression is written in this way because our treatment is symmetric
with respect to $z_{1,2}$. In the same plot, we also notice that once $z_{1}$
is chosen to be $2$ instead of $1$, $|\triangle C|$ increases, but the linear
behaviour remains unchanged.

Figure 7: Behaviour of $\frac{1}{2}\left(C_{1}+C_{2}\right)-\mathcal{C}_{P}$
with $z_{2}$ in the small mass regime.

##  4 Conclusion

We have studied the pseudo complexity of purification of the reduced
transition matrix (2.14) numerically for free scalar and Lifshitz scalar field
theories. We also defined the mutual pseudo complexity of purification by
comparing the pseudo complexity of purification with the complexities of
purification of the individual reduced density matrices of the states
$|\psi_{1}\rangle$ and $|\psi_{2}\rangle$. The main findings of this work are
listed below.

  1. 1.

Although the pseudo complexity is supposed to be a simple sum of the two
states involved in the construction of the transition matrix, the pseudo
complexity of purification can show non-trivial behaviour compared to the
individual complexities of purification associated with the reduced density
matrices of the two states involved.

  2. 2.

We also find that the mutual pseudo complexity of purification satisfies the
generally expected inequality $\Delta C(\tau_{A}^{1|2})\leq 0$ [30] 555Note
that the equality sign is different from [30] due to our choice of
definition.. Also, for all scaling exponents and masses corresponding to the
two states $|\psi_{1}\rangle$ and $|\psi_{2}\rangle$, the saturation of the
pseudo complexity of purification is similar to the saturation of half the sum
of individual complexities of purification for the reduced density matrices
$\rho_{A}^{1}$ and $\rho_{A}^{2}$ associated to the initial and final states.
More details on the behaviour of mutual pseudo complexity of purification are
there in Sec. (3.2).

  3. 3.

From the behaviour of pseudo complexity of purification, it appears that the
qualitative behaviour remains similar to that of the usual complexity of
purification of a reduced density matrix. Hence, we expect the general
dependence of the quantity in terms of the system parameters to be similar to
that of the complexity of purification of the reduced density matrices.
However, the pseudo complexity of purification, since it carries information
of both the states $|\psi_{1}\rangle$ and $|\psi_{2}\rangle$, is expected to
depend on parameters of two states instead of one. This also results in an
increment in the number of auxiliary parameters to be added and minimized for
finding the most optimal purification. As we find, the dependence on the two
states should be symmetric. It is not surprising since the transition matrix,
or more specifically, the reduced transition matrix by structure, does not
differentiate in any way between the two states used to define it.

  4. 4.

For the most general purification of the reduced transition matrix with $n$
oscillators, the number of unknown parameters scales as $2n^{2}$. This may be
understood from eqs. (2.12); for a reduced system made of $n$ linear harmonic
oscillators, the $J$, $K$, and $L$ matrices would be square matrices of rank
$n$. The choice (2.12) always leave any two of the six matrices undetermined,
and thus the number of unknown parameters in the final optimization problem is
$2n^{2}$.

Our analysis is mostly done for a reduced transition matrix of one oscillator
after integrating the rest. This analysis can be extended to more number of
reduced oscillators in many ways. However, due to the increase of unknown
parameters, the complexity optimization problem becomes numerically more
challenging, in general. However, one can perform a mode-by-mode purification
as done in [30] on a physical basis. There, instead of optimizing all the
modes together, one purifies each of the modes individually and then optimizes
the sum of the complexities for all the modes. This decreases the number of
unknowns significantly. For mode-by-mode purifications, the number of unknowns
will scale as just $2n$ as each mode involves $2$ auxiliary parameters.
However, since we get similar quantitative behaviour for pseudo complexity of
purification as usual complexity of purification for the reduced system of one
oscillator, we expect this to hold for the reduced system with more number of
modes. With this assumption, we can therefore use a similar form for the
pseudo complexity of purification, in terms of variables of two states instead
of one

$\mathcal{C}_{P}\left(\tau_{A}^{1|2}\right)=\frac{N}{2}\log\left(\omega_{2}\,\delta\right)+\frac{1}{2}\log\left(\frac{1}{f_{1}\left(m_{1},z_{1},m_{2},z_{2}\right)N\delta}\right)-\frac{f_{2}\left(m_{1},z_{1},m_{2},z_{2}\right)N^{2}\delta^{2}}{48}+\cdots,$
(4.1)

where $\delta$ is the gap between two points on the lattice 666This variable
does not come in our study since we did not consider any disjoint systems.,
this expression should work for very small mass, and the ellipses denote
higher order terms. Usually, $f_{1}\sim\mathcal{O}\left(m_{1,2}\right)$, and
$f_{2}\sim\mathcal{O}\left(m_{1,2}^{2}\right)$. The functions $f_{1}$ and
$f_{2}$ are symmetric under the exchange $1\leftrightarrow 2$, but arbitrary
otherwise 777One possibility is
$f_{1}=\frac{1}{2}\sqrt{m_{1}^{2z_{1}}+m_{2}^{2z_{2}}}$ and
$f_{2}=f_{1}^{2}$.. Apart from the symmetric dependence, the pseudo complexity
of purification has to follow all the necessary properties as mentioned in
Sec. (3) for both small and large mass regime. Also, the nontrivial dependence
on the scaling exponents have to abide by the plots we have found and reported
in the previous section. A more detailed study for subsystems with considering
larger size (by taking mode by mode purification), disjoint nature etc. can
yield more exact behaviour which is a very natural extension to this work.

There are a bunch of interesting questions one can ask in this direction.
Since the post-selected states might be, in general result of a non-unitary
evolution, it might be related to studies of complexity for open systems where
the Lindbladian evolution is non-unitary. One recent study in this direction
is [46], where initial steps towards studying Krylov complexity have been
taken. It would be interesting if one could relate the notion of Nielsen
complexity with the Krylov complexity along these lines. Another natural
direction is to study the pseudo complexity of purification for spin systems,
e.g. in the transverse field Ising model and check if the complexity can probe
phase transitions like pseudo entropy 888This can be further compared with the
study of entanglement of purification for transverse field Ising model as done
in [47] as well as with the corresponding study of phase transitions using
circuit complexity. Interested readers are referred to [48] and the citations
of it for the application of circuit complexity to detect phase transitions in
various quantum systems.. It would also be interesting to study conformal
systems as done for usual complexity of purification in [32]. Finally, some
holographic notions of pseudo entropy by computing minimal surfaces in
Euclidean setups were proposed in [1]. Therefore, it would be natural to
compute the volumes below those minimal surfaces and compute a holographic
pseudo subregion complexity (similar to holographic subregion complexity [49,
50, 51, 52, 53]) in those geometries. We hope to address some of these
problems later.

## References

  * [1] Y. Nakata, T. Takayanagi, Y. Taki, K. Tamaoka and Z. Wei, _New holographic generalization of entanglement entropy_ , _Phys. Rev. D_ 103 (2021) 026005 [2005.13801].
  * [2] A. Mollabashi, N. Shiba, T. Takayanagi, K. Tamaoka and Z. Wei, _Pseudo Entropy in Free Quantum Field Theories_ , _Phys. Rev. Lett._ 126 (2021) 081601 [2011.09648].
  * [3] J. M. Maldacena, _The Large N limit of superconformal field theories and supergravity_ , _Adv. Theor. Math. Phys._ 2 (1998) 231 [hep-th/9711200].
  * [4] E. Witten, _Anti-de Sitter space and holography_ , _Adv. Theor. Math. Phys._ 2 (1998) 253 [hep-th/9802150].
  * [5] R. Jefferson and R. C. Myers, _Circuit complexity in quantum field theory_ , _JHEP_ 10 (2017) 107 [1707.08570].
  * [6] S. Chapman, M. P. Heller, H. Marrochio and F. Pastawski, _Toward a Definition of Complexity for Quantum Field Theory States_ , _Phys. Rev. Lett._ 120 (2018) 121602 [1707.08582].
  * [7] L. Hackl and R. C. Myers, _Circuit complexity for free fermions_ , _JHEP_ 07 (2018) 139 [1803.10638].
  * [8] R. Khan, C. Krishnan and S. Sharma, _Circuit Complexity in Fermionic Field Theory_ , _Phys. Rev. D_ 98 (2018) 126001 [1801.07620].
  * [9] A. Bhattacharyya, A. Shekar and A. Sinha, _Circuit complexity in interacting QFTs and RG flows_ , _JHEP_ 10 (2018) 140 [1808.03105].
  * [10] A. Bhattacharyya, P. Nandy and A. Sinha, _Renormalized Circuit Complexity_ , _Phys. Rev. Lett._ 124 (2020) 101602 [1907.08223].
  * [11] T. Ali, A. Bhattacharyya, S. Shajidul Haque, E. H. Kim and N. Moynihan, _Time Evolution of Complexity: A Critique of Three Methods_ , _JHEP_ 04 (2019) 087 [1810.02734].
  * [12] J. M. Magán, _Black holes, complexity and quantum chaos_ , _JHEP_ 09 (2018) 043 [1805.05839].
  * [13] P. Caputa and J. M. Magan, _Quantum Computation as Gravity_ , _Phys. Rev. Lett._ 122 (2019) 231302 [1807.04422].
  * [14] A. Bhattacharyya, W. Chemissany, S. Shajidul Haque and B. Yan, _Towards the web of quantum chaos diagnostics_ , _Eur. Phys. J. C_ 82 (2022) 87 [1909.01894].
  * [15] T. Ali, A. Bhattacharyya, S. S. Haque, E. H. Kim, N. Moynihan and J. Murugan, _Chaos and Complexity in Quantum Mechanics_ , _Phys. Rev. D_ 101 (2020) 026021 [1905.13534].
  * [16] A. Bhattacharyya, S. Das, S. Shajidul Haque and B. Underwood, _Cosmological Complexity_ , _Phys. Rev. D_ 101 (2020) 106020 [2001.08664].
  * [17] A. Bhattacharyya, S. Das, S. S. Haque and B. Underwood, _Rise of cosmological complexity: Saturation of growth and chaos_ , _Phys. Rev. Res._ 2 (2020) 033273 [2005.10854].
  * [18] A. Bhattacharyya, W. Chemissany, S. S. Haque, J. Murugan and B. Yan, _The Multi-faceted Inverted Harmonic Oscillator: Chaos and Complexity_ , _SciPost Phys. Core_ 4 (2021) 002 [2007.01232].
  * [19] J. Erdmenger, M. Gerbershagen and A.-L. Weigel, _Complexity measures from geometric actions on Virasoro and Kac-Moody orbits_ , _JHEP_ 11 (2020) 003 [2004.03619].
  * [20] M. Flory and M. P. Heller, _Geometry of Complexity in Conformal Field Theory_ , _Phys. Rev. Res._ 2 (2020) 043438 [2005.02415].
  * [21] M. Flory and M. P. Heller, _Conformal field theory complexity from Euler-Arnold equations_ , _JHEP_ 12 (2020) 091 [2007.11555].
  * [22] J. Erdmenger, M. Flory, M. Gerbershagen, M. P. Heller and A.-L. Weigel, _Exact Gravity Duals for Simple Quantum Circuits_ , 2112.12158.
  * [23] N. Chagnet, S. Chapman, J. de Boer and C. Zukowski, _Complexity for Conformal Field Theories in General Dimensions_ , _Phys. Rev. Lett._ 128 (2022) 051601 [2103.06920].
  * [24] R. d. M. Koch, M. Kim and H. J. R. Van Zyl, _Complexity from spinning primaries_ , _JHEP_ 12 (2021) 030 [2108.10669].
  * [25] A. Bhattacharyya, G. Katoch and S. R. Roy, _Complexity of warped conformal field theory_ , 2202.09350.
  * [26] S. Chapman and G. Policastro, _Quantum Computational Complexity – From Quantum Information to Black Holes and Back_ , 2110.14672.
  * [27] A. Bhattacharyya, _Circuit complexity and (some of) its applications_ , _Int. J. Mod. Phys. E_ 30 (2021) 2130005.
  * [28] M. Ghodrati, X.-M. Kuang, B. Wang, C.-Y. Zhang and Y.-T. Zhou, _The connection between holographic entanglement and complexity of purification_ , _JHEP_ 09 (2019) 009 [1902.02475].
  * [29] C. A. Agon, M. Headrick and B. Swingle, _Subsystem complexity and holography_ , _Journal of High Energy Physics_ 2019 (2019) .
  * [30] E. Caceres, S. Chapman, J. D. Couch, J. P. Hernandez, R. C. Myers and S.-M. Ruan, _Complexity of Mixed States in QFT and Holography_ , _JHEP_ 03 (2020) 012 [1909.10557].
  * [31] H. A. Camargo, P. Caputa, D. Das, M. P. Heller and R. Jefferson, _Complexity as a novel probe of quantum quenches: universal scalings and purifications_ , _Phys. Rev. Lett._ 122 (2019) 081601 [1807.07075].
  * [32] H. A. Camargo, L. Hackl, M. P. Heller, A. Jahn, T. Takayanagi and B. Windt, _Entanglement and complexity of purification in ( 1+1 )-dimensional free conformal field theories_ , _Phys. Rev. Res._ 3 (2021) 013248 [2009.11881].
  * [33] A. Bhattacharyya, S. S. Haque and E. H. Kim, _Complexity from the reduced density matrix: a new diagnostic for chaos_ , _JHEP_ 10 (2021) 028 [2011.04705].
  * [34] A. Bhattacharyya, T. Hanif, S. S. Haque and M. K. Rahman, _Complexity for an open quantum system_ , _Phys. Rev. D_ 105 (2022) 046011 [2112.03955].
  * [35] M. A. Nielsen, _A geometric approach to quantum circuit lower bounds_ , _arXiv e-prints_ (2005) quant [quant-ph/0502070].
  * [36] M. A. Nielsen, M. R. Dowling, M. Gu and A. C. Doherty, _Quantum computation as geometry_ , _Science_ 311 (2006) 1133.
  * [37] M. R. Dowling and M. A. Nielsen, _The geometry of quantum computation_ , _arXiv e-prints_ (2006) quant [quant-ph/0701004].
  * [38] A. Jamiołkowski, _Linear transformations which preserve trace and positive semidefiniteness of operators_ , _Reports on Mathematical Physics_ 3 (1972) 275 .
  * [39] M.-D. Choi, _Completely positive linear maps on complex matrices_ , _Linear Algebra and its Applications_ 10 (1975) 285 .
  * [40] M. Jiang, S. Luo and S. Fu, _Channel-state duality_ , _Phys. Rev. A_ 87 (2013) 022310.
  * [41] A. Mollabashi, N. Shiba, T. Takayanagi, K. Tamaoka and Z. Wei, _Aspects of pseudo entropy in field theories_ , _Phys. Rev. Res._ 3 (2021) 033254 [2106.03118].
  * [42] M. Guo, J. Hernandez, R. C. Myers and S.-M. Ruan, _Circuit Complexity for Coherent States_ , _JHEP_ 10 (2018) 011 [1807.07677].
  * [43] A. Bhattacharyya, T. Takayanagi and K. Umemoto, _Entanglement of Purification in Free Scalar Field Theories_ , _JHEP_ 04 (2018) 132 [1802.09545].
  * [44] M. Alishahiha, K. Babaei Velni and M. R. Mohammadi Mozaffar, _Black hole subregion action and complexity_ , _Phys. Rev. D_ 99 (2019) 126016 [1809.06031].
  * [45] E. Cáceres, J. Couch, S. Eccles and W. Fischler, _Holographic Purification Complexity_ , _Phys. Rev. D_ 99 (2019) 086016 [1811.10650].
  * [46] A. Bhattacharya, P. Nandy, P. P. Nath and H. Sahu, _Operator growth and Krylov construction in dissipative open quantum systems_ , 2207.05347.
  * [47] A. Bhattacharyya, A. Jahn, T. Takayanagi and K. Umemoto, _Entanglement of Purification in Many Body Systems and Symmetry Breaking_ , _Phys. Rev. Lett._ 122 (2019) 201601 [1902.02369].
  * [48] T. Ali, A. Bhattacharyya, S. Shajidul Haque, E. H. Kim and N. Moynihan, _Post-Quench Evolution of Complexity and Entanglement in a Topological System_ , _Phys. Lett. B_ 811 (2020) 135919 [1811.05985].
  * [49] A. Bhattacharya, K. T. Grosvenor and S. Roy, _Entanglement Entropy and Subregion Complexity in Thermal Perturbations around Pure-AdS Spacetime_ , _Phys. Rev. D_ 100 (2019) 126004 [1905.02220].
  * [50] A. Bhattacharya, A. Chanda, S. Maulik, C. Northe and S. Roy, _Topological shadows and complexity of islands in multiboundary wormholes_ , _JHEP_ 02 (2021) 152 [2010.04134].
  * [51] A. Bhattacharya, A. Bhattacharyya, P. Nandy and A. K. Patra, _Islands and complexity of eternal black hole and radiation subsystems for a doubly holographic model_ , _JHEP_ 05 (2021) 135 [2103.15852].
  * [52] A. Bhattacharya, A. Bhattacharyya, P. Nandy and A. K. Patra, _Partial islands and subregion complexity in geometric secret-sharing model_ , _JHEP_ 12 (2021) 091 [2109.07842].
  * [53] A. Bhattacharya, A. Bhattacharyya, P. Nandy and A. K. Patra, _Bath deformations, islands, and holographic complexity_ , _Phys. Rev. D_ 105 (2022) 066019 [2112.06967].