# Online Cluster Validity Indices for
Streaming Data

Masud Moshtaghi, James C. Bezdek, Sarah M. Erfani, Christopher Leckie, James
Bailey M. Moshtaghi, J.C. bezdek, S. Erfani, C. Leckie and J.Bailey were with
the School of Computing and Information Systems, University of Melbourne,
Melbourne, Australia, 3010. M. Moshtaghi is now at Amazon.com
E-mail<EMAIL_ADDRESS>
###### Abstract

Cluster analysis is used to explore structure in unlabeled data sets in a wide
range of applications. An important part of cluster analysis is validating the
quality of computationally obtained clusters. A large number of different
internal indices have been developed for validation in the offline setting.
However, this concept has not been extended to the online setting. A key
challenge is to find an efficient incremental formulation of an index that can
capture both cohesion and separation of the clusters over potentially infinite
data streams. In this paper, we develop two online versions (with and without
forgetting factors) of the Xie-Beni and Davies-Bouldin internal validity
indices, and analyze their characteristics, using two streaming clustering
algorithms (sk-means and online ellipsoidal clustering), and illustrate their
use in monitoring evolving clusters in streaming data. We also show that
incremental cluster validity indices are capable of sending a distress signal
to online monitors when evolving clusters go awry. Our numerical examples
indicate that the incremental Xie-Beni index with forgetting factor is
superior to the other three indices tested.

##  1 Introduction

The intrinsic nature of streaming data requires algorithms that are capable of
fast data analysis to extract knowledge. Online clustering algorithms provide
a way to extract patterns from continuous data streams. Therefore, online
clustering has gained popularity in applications involving massive streams of
data, such as router packet analysis and environmental sensing [1, 2]. In
these applications, the velocity and volume of data is too high for the
processing unit to access each data sample more than once. A category of fast
online clustering algorithms, also referred to as sequential clustering,
process data quickly and efficiently by receiving samples one at a time and
updating cluster statistics (such as prototypes) with each sample [3, 4, 5].

An important aspect of any clustering algorithm is assessment of the quality
of the resulting clusters, i.e., how well does any partition match the input
data? Cluster validity indices (CVIs) comprise computational models and
algorithms whose job is to identify the “best” member among different
partitions of the input data. Most CVIs are max-optimal or min-optimal,
meaning that the partition preferred by the index is indicated by the maximum
(minimum) value of the index on the partitions being evaluated. To date, CVIs
have been used exclusively in a static (or batch) setting, being applied to
sets of partitions generated by different parameter settings of the clustering
algorithm applied to a collected set of data. In online clustering, the
question of how well a clustering algorithm matches the data becomes even more
important as we need to know how the clustering algorithm is performing at any
point in time and how it reacts to the changes in the partitions over time.
The key assumption in the online context is that data are processed once, and
historical data will not be available for a retrospective analysis.

While batch clustering algorithms are supported by a wide range of cluster
validation methods [6], no effort has been made to extend cluster validation
methods to online clustering algorithms. This study concerns itself with the
use of internal incremental cluster validity indices (iCVIs) (we drop the word
“internal” in the rest of the paper), that are computed online (on the fly)
and used to control/interpret clustering for streaming data.

Most of existing batch indices assess two basic characteristics of a set of
clusters: compactness (or density) and separation of the clusters [6].
Compactness is usually calculated based on the observations while separation
is often measured by the distance between cluster prototypes. In the online
setting, where each observation can be accessed only once, an
incremental/recursive calculation of compactness is necessary. In this paper,
we propose two incremental methods to estimate within cluster dispersion: (1)
an exact incremental update of a batch formula; and (2) an online formula
incorporating an exponential forgetting factor. Using these two methods we
then derive online versions of two well-known CVIs namely, the Xie-Beni (XB)
[7] and the Davies-Bouldin (DB) [8] indices. These indices can be applied to
both hard and soft partitions.

This paper offers four main contributions: (1) we propose a new concept of
incremental validation of online clustering algorithms which provides new
insights into these algorithms; (2) we propose two incremental versions of
within cluster dispersion, viz., with and without forgetting, a facility that
enables the iCVI to gracefully forget earlier inputs; (3) we propose
incremental versions of two well-known batch CVIs allowing exact calculation
of the two indices with fast sequential processing of data; and (4) we analyze
and discuss the properties of the proposed iCVIs within the context of two
online clustering algorithms. Our results demonstrate useful insights produced
by the iCVIs about the partitions produced by online clustering algorithms.
Moreover, the proposed iCVIs can indicate learning difficulties experienced by
the clustering algorithm and can signal the appearance of new clusters in the
evolution of the data stream.

The next section summarizes related work. In Section 3, we present definitions
and notation needed in this paper. Section 4 contains background information
on two important online clustering algorithms. In Section 5, we derive two
versions of the iXB, iDB indices and analyze their characteristics. Section 7
contains numerical examples used to evaluate the proposed models. A summary
and conclusions are given in Section 8.

##  2 Background and Related Work

In this section, we briefly describe related work in cluster validation and
online clustering algorithms. It is especially important to record our
definition of an online clustering algorithm as our goal is to validate these
algorithms using cluster validity indices.

According to Guha et al. [1], clustering in data streams can be divided into
two main strategies: (1) buffering a window of streamed inputs, finding
clusters in the window using a batch algorithm such as the classic k-means
algorithm; and then merging clusters in adjacent windows to obtain a final
clustering. This strategy is espoused, for example, in [9, 10, 11]; (2) using
incremental learning techniques to find clusters in the evolving data stream.
Examples of this strategy include [3, 5, 12, 13, 14, 15]. We refer to this
second approach as online or incremental clustering. Algorithms for online
clustering can themselves be divided into two categories. The first category
is general clustering algorithms for any sequence of data (we call this type
2a data). These algorithms do not assume any ordering in the data stream and
require the number of clusters to be specified in advance, for example
sequential k-means, or more briefly, sk-means, and sequential agglomerative
clustering [3]. A second category of online clustering algorithms assume a
natural ordering in the data (time-series) and operate on the assumption that
close observations in time will be closely related. These algorithms use this
assumption to dynamically create clusters in evolving data streams (we call
this type 2b data). In this paper, we propose cluster validation methods for
online clustering algorithms of types 2a and 2b.

Cluster validity indices (CVIs) can be grouped into two categories: internal
and external indices. Internal indices use only the information available from
the algorithmic outputs and the observed unlabeled data. In contrast, external
CVIs use additional external information about substructure in the data,
usually in the form of a reference partition (a ground truth partition), so
the data processed are labeled. External CVIs are used to compare partitions
obtained by a clustering algorithm to ground truth labels. Another use of
external CVIs is to correlate external and internal assessments of labeled
data [6]. In this application the external CVI becomes a tool for selection of
the “best” internal CVI when unlabeled data are to be clustered. The use of
external CVIs to choose a “good” internal CVI is comprehensively discussed in
[6]. Our focus in this paper is on internal CVIs. Milligan and Cooper’s 1985
paper is widely regarded as a landmark study for the comparison of internal
CVIs [16].

Internal cluster validity indices fall under the broad umbrella of goodness-
of-fit measures such as likelihood ratio tests and Root Mean Squared Error
[17]. The majority of goodness-of-fit measures target parametric models while
CVIs provide a non-parametric mechanism to evaluate clustering outputs.
Another differentiating point between internal CVIs and goodness-of-fit
measures is the meaning of the fitness term. Most internal CVI models have
components that attempt to capture cohesion and separation, while goodness of
fit indices usually assess the fit of a model to the data that generates it.

There are two categories of internal CVIs based on the way that they measure
cohesion and separation. The CVIs in the first category use only the
partitions generated by the clustering to determine the quality of the
partition. Measures of this type include the partition coefficient and
partition entropy [18]. Indices such as these often appear in the context of
fuzzy cluster validity. However, most CVIs fall into the second category, that
is, they use both the data and the partition information.

We develop incremental CVIs by deriving an incremental formula for the data-
dependent part of two well-known indices, i.e., XB and DB. After determining
the incremental update formula for cluster cohesion, we derive one step update
formulae for these two indices. We then investigate their application to
cluster analysis in data streams. While the results of any online clustering
algorithm can be analyzed using the proposed iCVIs, we focus on two clustering
algorithms - a crisp clustering algorithm from the general data stream
clustering category (type 2a), viz., sk-means; and the online elliptical
clustering (OEC) clustering algorithm from the time-series clustering category
(type 2b), which is a soft/fuzzy clustering algorithm.

##  3 Problem Statement and Definitions

Traditional batch clustering algorithms aim to find crisp or
fuzzy/probabilistic k-partitions of a collection of $n$ samples of static
data, viz.,
$X=\left\\{\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{n},\ldots,\bm{x}_{n}\right\\}\subset\Re^{p}$.
All vectors in this article are column vectors. Crisp and fuzzy partitions of
$X$ are conveniently represented by matrices in the following sets:

$\begin{split}M_{\textit{fkn}}=&\left\\{\right.U\in\Re^{kn}:\textit{for }1\leq
i\leq k,1\leq j\leq n:\\\ &0\leq u_{ij}\leq 1:\sum_{i=1}^{k}u_{ij}=1\;\forall
j;\;\sum_{j=1}^{n}u_{ij}>0\;\forall i\left.\right\\};\end{split}$ (1)
$M_{\textit{hkn}}=\left\\{U\in
M_{\textit{fkn}}:u_{ij}\in\left\\{0,1\right\\}\forall i,j\right\\}.$ (2)

Now suppose that $n$ inputs have arrived sequentially in the streaming data,
and we have found, from these $n$ inputs $U_{n}\in M_{\textit{fkn}}$, a set of
soft/fuzzy clusters of $X$, together with a set
$V_{n}=\left\\{\bm{v}_{1n},\ldots,\bm{v}_{kn}\right\\}\subset\Re^{kp}$ of
cluster centers. We can use $\left(U_{n},V_{n}\right)$ to calculate various
CVIs.

When input $\bm{x}_{n+1}$ arrives, it is used by an online clustering
algorithm to find the membership, $u_{i,n+1}$, of the new point in the
$i^{\textit{th}}$ cluster, $1\leq i\leq k$. Let the vector
$\bm{u}_{n+1}=\left\\{u_{i,n+1}|i=1,\ldots,k\right\\}$ be the label vector of
$\bm{x}_{n+1}$ in the set of ($k$) clusters. The clustering algorithm will
also use $\bm{x}_{n+1}$ to update $V_{n}\rightarrow V_{n+1}$.

Given the new input $\bm{x}_{n+1}$, its cluster assignment $\bm{u}_{n+1}$ and
the updated cluster centers $V_{n+1}$, we will derive one-step update formula
for a particular CVI.

The calculation of $\bm{u}_{n+1}$ and updates to $V$ are done using specific
incremental clustering algorithms such as sk-means [3] or OEC [5]. The
question posed here is how the value of the chosen CVI changes incrementally
with this update. Fig. 1 illustrates the overall process. The two time series
at the top of the figure form the input to the online clustering algorithm on
the bottom right. When $\bm{x}_{n+1}$ becomes available, the online clustering
algorithm finds the membership values $\bm{u}_{n+1}$, and updates the cluster
centers to produce $V_{n+1}$. Then, $\bm{u}_{n+1},V_{n+1},\bm{x}_{n+1}$ are
passed to the incremental cluster validation process. The objective in this
paper is to answer the question “iCVI(n+1)=?” posed in the bottom left panel
of Fig. 1.

Figure 1: One-step update of incremental CVI in an online setting.

To answer this question, we first need to describe how $\bm{u}_{n+1}$ and
$V_{n+1}$ are calculated incrementally by online clustering algorithms.
Section 4, briefly describes two algorithms to solve this problem.

##  4 Incremental Clustering Algorithms

The calculation of iCVIs depends on the information provided by a particular
clustering algorithm. Therefore, we start by providing a brief overview of two
different types of incremental clustering algorithms.

###  4.1 Sequential k-means

The sk-means algorithm shown below as Algorithm 1 and its variants have been
studied in various forms in the literature of self organizing maps (SOMs) [3,
19]. Algorithm 1 records the basic sequential k-means method studied by
Macqueen [20].

Data: $X$ \- set of data points

Input : $k$ \- number of clusters

Note : $\left\|.\right\|$ is the Euclidean norm

Initialize $V_{k}$ with the first $k$ data points
$V_{k}=\left\\{\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{k}\right\\}$;

$U_{k}$ is the $k\times k$ identity matrix;

Initialize the counter for each cluster
$\left\\{n_{1},n_{2},\ldots,n_{k}\right\\}$ with 1;

foreach  _$\bm{x}_{n}$ in the stream_ do

$m=\operatorname*{argmin}_{m\in\left\\{1\ldots
k\right\\}}\left\|\bm{x}_{n}-\bm{v}_{m}\right\|$;

$n_{m}=n_{m}+1$;

$\bm{v}_{m}=\bm{v}_{m}+\left(\bm{x}_{n}-\bm{v}_{m}\right)/n_{m}$;

$\bm{u}_{n}=\left(0,0,\ldots,1,\ldots,0\right)^{T}\in\left\\{0,1\right\\}^{k}$,
where the 1 appears in the $m-th$ place, is appended to the current $U$ as its
$n-th$ column;

end foreach

Algorithm 1 The basic sequential k-means algorithm (Macqueen [20]).

Macqueen’s algorithm requires pre-specification of $k$, the number of clusters
it builds. After initializing the $k$ clusters by designating the first $k$
points in $X$ as the initial cluster centers, when the next input arrives, sk-
means computes the distance from the next input to the $k$ cluster centers,
and assigns the input point to the cluster of the nearest prototype. Then the
winning prototype is updated as shown in the third line of the for part of
Algorithm 1. So at step $n+1$ of the algorithm the $k$ vector $\bm{u}_{n+1}$
has only one element with the value of 1 (the winner of a nearest prototype
competition), and the remaining $k-1$ values are 0.

###  4.2 Online Ellipsoidal Clustering

This section briefly reviews the online ellipsoidal clustering (OEC) algorithm
developed in [5]. This clustering algorithm is specific to time-series and
does not require the number of clusters to be pre-specified. This algorithm
has ellipsoidal cluster prototypes defined by $k$ sets of means and covariance
matrices $\left\\{(\bm{m}_{i,n},S^{-1}_{i,n})|i=1,\ldots,k\right\\}$ at any
step $n$. Two different radii are considered for each ellipsoid, one called
the effective boundary, i.e., the smaller radius, and the other one is called
the outlier boundary. The outlier boundary prevents the cluster prototype from
being affected by large outliers. These boundary thresholds are selected from
the inverse of the chi-squared distribution, $(\chi^{2})^{-1}_{p}(\gamma)$,
where $\gamma$ is the probability that a cluster member falls inside the
ellipsoid. In the OEC clustering algorithm, cluster center $\bm{v}_{i,n+1}$ is
formed from the sample mean $\bm{m}_{i,n+1}$ of each cluster, and
$\bm{u}_{n+1}$ is formed by the values obtained from the soft assignments of
the input point to the clusters, $\bm{u}_{i,n+1},\;1\leq i\leq k$, defined as

$\bm{u}_{i,n+1}=\left[\sum_{j=1}^{k}\left(F_{i,n+1}/F_{j,n+1}\right)^{2/(m-1)}\right]^{-1},\;\;\;m\in(1,\infty),$
(3)

where
$F_{i,n+1}=(\bm{x}_{n+1}-\bm{m}_{i,n})^{T}S^{-1}_{i,n}(\bm{x}_{n+1}-\bm{m}_{i,n})$.
Readers may recognize this as the membership update formula required by fuzzy
k-means using the sample-based Mahalanobis norm for the inner product induced
distance [21]. The OEC algorithm uses two other parameters that we need to set
during the experiments.

Stabilization period ($n_{s}$) -  An ellipsoidal prototype in $\Re^{p}$ is
created by $p+1$ consecutive distinct points, but to obtain a reliable
estimate of $(\bm{m}_{i,n},S^{-1}_{i,n})$, more data from the input stream is
required. The integer $n_{s}$ stabilizes these incremental estimates by
temporarily disabling the OEC guard zone and new cluster detection tests until
the current cluster contains $n_{s}$ points.

Forgetting factor ($\lambda_{\textit{OEC}}$)-  OEC has a special forgetful
prototype. When this prototype does not overlap with any of the existing
clusters in the system, a new cluster is formed. This prototype use an
exponential forgetting factor ($\lambda_{\textit{OEC}}$) for this purpose.

Now that we have described how $\bm{v}_{n+1}$ and $\bm{u}_{n+1}$ can be
updated using two different clustering algorithms, we can describe how iCVIs
use this information to assess the evolving clusters incrementally. In Section
5, we introduce an incremental calculation of a common cohesion measure in
CVIs, then we discuss how two particular CVIs which use this cohesion measure
can be updated using its previous value at time $n$ with the additional
information provided by $\bm{u}_{n+1},V_{n+1}$ and $\bm{x}_{n+1}$.

##  5 The Incremental CVIs

Normally there is no external ground truth information available with
streaming inputs, hence we explore the usage of internal CVIs to assess the
evolving performance of online clustering algorithms that generate both
memberships and cluster prototypes. We consider two well-known indices in this
group that work for both hard and soft partitions.

The first index of this type is the Xie-Beni (XB) index [7]. A general
formulation of this index for a batch collection of $n$ inputs is [22], where
$\mathscr{A}$ is a positive-definite weight matrix which induces the inner
product norm $\left\|x\right\|^{2}_{\mathscr{A}}=x^{T}\mathscr{A}x$:

$\begin{split}&\textit{XB}_{m\mathscr{A}}\left(U,V;X\right)=\frac{J_{m\mathscr{A}}\left(U,V;X\right)}{n\left(\min_{i\neq
j}\left\\{\left\|\bm{v}_{i}-\bm{v}_{j}\right\|_{\mathscr{A}}^{2}\right\\}\right)}\\\
&\text{, where }m\in[1,\infty)\text{, and }\\\
&J_{m\mathscr{A}}\left(U,V;X\right)=\sum_{j=1}^{n}\sum_{i=1}^{k}(u_{ij})^{m}\left\|\bm{x}_{n}-\bm{v}_{i}\right\|_{\mathscr{A}}^{2}.\end{split}$
(4)

The parameter $m$ is called the fuzzifier of the model. For simplicity, we
consider only $\mathscr{A}=I_{p}$ (the Euclidean norm)and $m=2$, and drop
subscripts, writing $J_{2I_{p}}$ as $J$ in the sequel.

The second index is a relative of the DB index [8] introduced by Araki  _et
al._ [23],

$\begin{split}&\textit{DB}\left(U,V;X\right)=\frac{1}{k}\sum_{i=1}^{k}{\max_{j,j\neq
i}{\frac{L_{i}+L_{j}}{\left\|\bm{v}_{i}-\bm{v}_{j}\right\|^{2}}}}\\\
&L_{i}=\frac{\sum_{j=1}^{n}{u_{ij}^{2}\left\|\bm{x}_{j}-\bm{v}_{i}\right\|^{2}}}{\sum_{j=1}^{n}{u_{ij}^{2}}}.\end{split}$
(5)

This is not a true generalization of the the DB index because square roots are
missing for this choice of $p=q=2$ in [8]. but the values of (5) are closely
related to the true DB index in the crisp case. The common part of these two
(and many other) indices is the way that they capture the within-cluster
dispersion. We define the fuzzy within cluster dispersion, for $U\in M_{fkn}$,
as,

$\begin{split}C_{i,n}=\sum_{j=1}^{n}{(u_{ij})^{2}\left\|\bm{x}_{j}-\bm{v}_{i,n}\right\|^{2}}.\end{split}$
(6)

This is the only part of either index that depends directly on the input data.
Therefore, we first develop an incremental calculation of this part of the
index, and then use it to propose formulae for incremental calculation of
these two indices.

###  5.1 Incremental Cluster Dispersion Measure

In this section, we derive a formula for (6) at time step $n+1$ based on its
value at time step $n$. We assume that at time $n+1$, all of the previous
input values $\bm{x}_{1},...,\bm{x}_{n}$ have been discarded, so the only data
point we have to work with is $\bm{x}_{n+1}$.The goal is to write the update
formula in terms of the value in the previous step and a change to this value
on seeing the $n+1^{\textit{st}}$ input, i.e., $C_{i,n+1}=C_{i,n}+\Delta
C_{i,n}$. We begin with

$\begin{split}C_{i,n+1}=\sum_{j=1}^{n+1}{(u_{ij})^{2}\left\|\bm{x}_{j}-\bm{v}_{i,n+1}\right\|^{2}}\end{split}$
(7)

First we isolate the effect of the last point in the summation to obtain

$\begin{split}C_{i,n+1}&=\sum_{j=1}^{n}{(u_{ij})^{2}\left\|\bm{x}_{j}-\bm{v}_{i,n+1}\right\|^{2}}+\\\
&\underbrace{(u_{i,n+1})^{2}\left\|\bm{x}_{n+1}-\bm{v}_{i,n+1}\right\|^{2}}_{A_{i,n+1}}.\end{split}$
(8)

At this point if we could assume that
$\left\|\bm{v}_{i,n+1}-\bm{v}_{i,n}\right\|\approx 0$, we would have
$C_{i,n+1}=C_{i,n}+A_{i,n+1}$. Here, we are looking for an exact calculation
of $C_{i,n+1}$ so we need to compute the effect of the change in the cluster
centers. To isolate the change of the centers, we add and subtract
$\bm{v}_{i,n}$ inside the Euclidean norm in the first term in (8), and rewrite
the norm in terms of the (Euclidean) inner product,

$\begin{split}&C_{i,n+1}=\sum_{j=1}^{n}{(u_{ij})^{2}<\bm{x}_{j}-\bm{v}_{i,n}+\bm{v}_{i,n}-\bm{v}_{i,n+1},\bm{x}_{j}-\bm{v}_{i,n}}\\\
&+\bm{v}_{i,n}-\bm{v}_{i,n+1}>+A_{i,n+1}.\end{split}$ (9)

After few steps and some simplification we obtain

$\begin{split}&C_{i,n+1}=C_{i,n}+\overbrace{A_{i,n+1}+M_{i,n}B_{i,n+1}+2Q_{i,n+1}}^{\Delta
C_{i,n}}\end{split}$ (10)

where

$\begin{split}&Q_{i,n+1}=\sum_{j=1}^{n}{(u_{ij})^{2}\left\langle\bm{x}_{j}-\bm{v}_{i,n},\bm{v}_{i,n}-\bm{v}_{i,n+1}\right\rangle}\end{split}$
(11)
$\begin{split}&B_{i,n+1}=\left\|\bm{v}_{i,n}-\bm{v}_{i,n+1}\right\|^{2}\end{split}$
(12)

$\begin{split}&M_{i,n+1}=M_{i,n}+u^{2}_{i,n+1},\;M_{i,1}=u^{2}_{i,1}\end{split}$
(13)

The term $Q_{i,n+1}$ in equation (11) depends on the previous (discarded)
values of $\bm{x}_{j}$ and $u_{ij}$ where $j=1,\ldots,n$, so we cannot yet
make the incremental calculation of $\Delta C_{i,n}$. Since the second part of
the dot product in (11), i.e., $\left(\bm{v}_{i,n}-\bm{v}_{i,n+1}\right)$,
does not depend on $j$, we can write $Q_{i,n+1}$ as

$\begin{split}Q_{i,n+1}=\left[\bm{v}_{i,n}-\bm{v}_{i,n+1}\right]^{T}\overbrace{\left[\sum_{j=1}^{n}{(u_{ij})^{2}\left(\bm{x}_{j}-\bm{v}_{i,n}\right)}\right]}^{G_{i,n}}\end{split}$
(14)

Using the same trick of adding and subtracting $\bm{v}_{i,n}$ in (14) we can
write an incremental update formula for $G_{i,n+1}$,

$\begin{split}&G_{i,n+1}=G_{i,n}+\Delta G_{i,n},G_{i,1}=\overrightarrow{0}\\\
&\Delta
G_{i,n}=M_{i,n}(\bm{v}_{i,n}-\bm{v}_{i,n+1})+(u_{i,n+1})^{2}\left(\bm{x}_{n+1}-\bm{v}_{i,n+1}\right).\end{split}$
(15)

We now have all the components needed to calculate $\Delta C_{i,n}$. The two
terms $A_{i,n+1}$ and $B_{i,n+1}$ are calculated directly and $Q_{i,n+1}$ and
$M_{i,n+1}$ are incrementally calculated. The Algorithm 2 depicts a function
that incrementally calculates the compactness.

Data: $\bm{v}_{i,n},\bm{v}_{i,n+1},u_{i,n+1}$, $\bm{x}_{n+1}$

Input : $G_{i,n},M_{i,n},C_{i,n}$

Output : $G_{i,n+1},M_{i,n+1},C_{i,n+1}$

/* note $i=1,\ldots,k$

foreach  _$i\in\left\\{1,\ldots,k\right\\}$_ do

$Q_{i,n+1}=\left[\bm{v}_{i,n}-\bm{v}_{i,n+1}\right]^{T}G_{i,n}$;

$B_{i,n+1}=\left\|\bm{v}_{i,n}-\bm{v}_{i,n+1}\right\|^{2}$;

$A_{i,n+1}=(u_{i,n+1})^{2}\left\|\bm{x}_{n+1}-\bm{v}_{i,n+1}\right\|^{2}$;

$C_{i,n+1}=C_{i,n}+A_{i,n+1}+M_{i,n}B_{i,n+1}+2Q_{i,n+1}$;

$G_{i,n+1}=G_{i,n}+M_{i,n}(\bm{v}_{i,n}-\bm{v}_{i,n+1})+$
$(u_{i,n+1})^{2}\left(\bm{x}_{n+1}-\bm{v}_{i,n+1}\right)$;
$M_{i,n+1}=M_{i,n}+u^{2}_{i,n+1}$;

end foreach

Algorithm 2 The incremental compactness function for data point
$\bm{x}_{n+1}$.

As $n\rightarrow\infty$ the effect of each new sample on the total value of
$C_{i,n}$ is expected to become small. In data streaming applications, which
have, in theory, an infinite data stream, we are interested in how well a
clustering algorithm keeps up with the evolution of cluster structure of
points in the stream. Our objective is to use iCVIs to capture the quality of
the partitions over a window of the most recent observations at any point in
time.

Exponential fading memory is a common approach in online learning methods. In
this approach a forgetting factor $0<\lambda<1$ is included in the incremental
estimations so that the data sample from $f$ steps before the current sample
is weighted by $\lambda^{f}$. In this way, older samples become less and less
relevant to the current estimation. The batch representation of $C_{i,n}$ with
a forgetting is shown in (16).

$C_{\lambda
i,n}=\sum_{j=1}^{n}{\lambda^{n-j}u_{ij}^{2}\left\|\bm{x}_{j}-\bm{v}_{i,n}\right\|^{2}}$
(16)

An argument similar to the one used to derive equation (10) leads to an
incremental update formula for $C_{\lambda i,n+1}$.

$\begin{split}&Q_{\lambda
i,n+1}=\left(\bm{v}_{i,n}-\bm{v}_{i,n+1}\right)G_{\lambda i,n}\end{split}$
(17)

$\begin{split}B_{i,n+1}=\left\|\bm{v}_{i,n}-\bm{v}_{i,n+1}\right\|^{2}\end{split}$
(18)

$\begin{split}A_{i,n+1}=u^{2}_{i,n+1}\left\|\bm{x}_{n+1}-\bm{v}_{i,n+1}\right\|^{2}\end{split}$
(19)

$\begin{split}&C_{\lambda i,n+1}=\lambda C_{\lambda i,n}+\Delta C_{\lambda
i,n}\\\ &\Delta C_{\lambda i,n}=2\lambda Q_{\lambda i,n+1}+\lambda M_{\lambda
i,n}B_{i,n+1}+A_{i,n+1}\end{split}$ (20)

$\begin{split}&G_{\lambda i,n+1}=\lambda G_{\lambda i,n}+\Delta G_{\lambda
i,n}\\\ &\Delta G_{\lambda i,n}=\lambda M_{\lambda
i,n}(\bm{v}_{i,n}-\bm{v}_{i,n+1})+(u_{i,n+1})^{2}\left(\bm{x}_{n+1}-\bm{v}_{i,n+1}\right).\end{split}$
(21)

$\begin{split}M_{\lambda i,n+1}=\lambda
M_{i,n}+u_{i,n}^{2},M_{i,1\lambda}=u_{i,1}^{2}\end{split}$ (22)

In the next section, we use these formulas to derive two incremental versions
of the XB and DB indices at equations (4) and (5).

###  5.2 Incremental Xie-Beni Index

Let XB($n+1$) denote the value of the XB index we seek when $\bm{x}_{n+1}$
arrives after $n1$ inputs have been processed. Let $J_{n+1}$ denote the value
of $J$ at step $n+1$. To compute XB incrementally, i.e., to compute XB($n+1$),
we need an incremental update for $J_{n+1}$ and for the denominator of (4).
The numerator of XB is updated using the value of $C_{i,n+1}$ from Algorithm
2, and the denominator is calculated with the updated centers $V_{n+1}$.
Equation (23) shows the one step update of $J_{n+1}$ at step $n+1$.

$J_{n+1}=\sum_{i=1}^{k}{C_{i,n+1}}.$ (23)

Let

$h_{n+1}=\min_{i\neq j,\bm{v}_{i},\bm{v}_{j}\in
V_{n+1}}\left\\{\left\|\bm{v}_{i}-\bm{v}_{j}\right\|^{2}\right\\}.$ (24)

Then value of the incremental XB index, is

$\begin{split}\textit{XB}(n+1)=\frac{J_{n+1}}{(n+1)h_{n+1}}.\end{split}$ (25)

For batch clustering, the case $k=1$ is usually not considered. However, in
the streaming environment and in algorithms like OEC, the number of clusters
dynamically changes and starts from $k=1$. With one cluster, $h_{n+1}$ is
undefined in (24). When $k=1$, we replace (24) with
$h_{n+1}=\max{\left\\{h_{n},\left\|\bm{v}_{1}-\bm{x}_{n+1}\right\|^{2}\right\\}}$.

The one step update of iXB with forgetting is obtained by replacing
$C_{i,n+1}$ with $C_{\lambda i,n+1}$ in (23) to obtain in $J_{\lambda,n+1}$,

$\begin{split}\textit{XB}_{\lambda}(n+1)=\frac{(1-\lambda)J_{\lambda,n+1}}{h_{n+1}}.\end{split}$
(26)

Please note that XB($n+1$) and XBλ($n+1$) are the values of the incremental XB
indices without and with the forgetting factor after $\bm{x}_{n+1}$ is
processed, while iXB and iXBλ are the names of the incremental XB models.

###  5.3 Incremental DB Index

Let DB($n$) denote the value of the Davies-Bouldin index after $n$ inputs. We
want to compute incrementally updated values DB($n+1$) and DBλ($n+1$) when
input $\bm{x}_{n+1}$ arrives. We need to normalize $C_{i,n+1}$ and $C_{\lambda
i,n+1}$ with the number of data points in the $i^{\textit{th}}$ cluster. The
index DB($n+1$) at time step $n+1$ can be written as

$\begin{split}\textit{DB}(n+1)=\frac{1}{k}\sum_{i=1}^{k}{\max_{j,j\neq
i}{\frac{L_{i,n+1}+L_{j,n+1}}{\left\|\bm{v}_{i,n+1}-\bm{v}_{j,n+1}\right\|^{2}}}},\end{split}$
(27)

where

$\begin{split}L_{i,n+1}=\frac{C_{i,n+1}}{M_{i,n+1}}.\end{split}$ (28)

To calculate the index with the forgetting factor, DB${}_{\lambda}(n+1)$, we
need to define $L_{\lambda i,n+1}$ to be used instead of $L_{i,n+1}$ and we
control the decay of the denominator by clamping the decay at 1 with $\max$
function.

$\begin{split}L_{\lambda i,n+1}=\frac{C_{\lambda
i,n+1}}{\max{\left\\{1,M_{\lambda i,n+1}\right\\}}}.\end{split}$ (29)

(a) S1: Synthetic Locally Linear Processes

(b) S2: Shifting Gaussian Distribution

(c) S3: Circular Shifting Gaussian Distribution with Systematic Noise

Figure 2: Scatter plots of three synthetic datasets. S2 and S3 show the progression of time with colors, starting with black and becoming lighter with time. TABLE I: Summary characteristics of the datasets used in the evaluations. Dataset | # samples ($n$) | # Labels ($k$) | # dimensions ($p$) | Noise Attributes | Labeling
---|---|---|---|---|---
S1 | 1955 | 2 | 2 | No Noise - Drift between Clusters | Exact
S2 | 2727 | 11 | 2 | 1% - Uniform Random | Exact
S3 | 2000 | 10 | 2 | 1% - Systematic Random | Exact
LG | 2016 | 3∗ | 2 | Unknown | Our guess∗
GSA | 9969 | 3∗ | 8 | Unknown | Our guess∗
∗ These are unlabeled data sets, so we use physical arguments to justify the approximate |
number of clusters shown in column 2 for these two data sets. |

As with the XB case, DB($n$) and DB($n+1$) are values of the incremental
indices; iDB and iDBλ𝛌 are the incremental models that produce them.

##  6 Computational Protocols

In this section, we first describe the synthetic and real-life datasets used
in our evaluations and then we study the characteristics of iXB,
$\text{XB}_{\lambda}(n)$, iDB and $\text{DB}_{\lambda}(n)$.

###  6.1 Datasets and Parameters

The synthetic dataset $S1$ consists of two-dimensional vectors
$(x_{q},y_{q}),q>2$, which are generated using two modes, $M_{1}$ and $M_{2}$,
with different dynamic functions and input signals ($x_{n},q=3,4,\ldots$).
Values of the independent variable ($x$) are random i.i.d. samples from a
Gaussian distribution with $\mu=\sigma=1$. Values of the dependent variable
($y$) from $M_{1}$ and $M_{2}$ are then computed according to (30) or (31)
respectively.

$\begin{split}&y_{n}=1.018x_{n-1}+1.801y_{n-1}-0.8187y_{n-2}\\\
&y_{0}=y_{1}=y_{2}=0\end{split}$ (30)

$\begin{split}&y_{n}=x_{n-1}+0.5x_{n-2}+1.5y_{n-1}-0.7y_{n-2}\\\
&y_{0}=y_{1}=y_{2}=0\end{split}$ (31)

To build $S1$, we considered 4 mode changes between the two modes at uniform
random intervals between 200 and 500 samples starting with $M_{1}$. Instead of
a sudden shift between the modes, we gradually change the individual
parameters of one mode to the other mode in 5 equal steps, during which we
generate 10 samples in each intermediate mode. Fig. 2(a) shows a scatter plot
of the $S1$ dataset with the input (Feature 1 = $x$) and output (Feature 2 =
$y$).

The second synthetic dataset $S2$ (shown in Fig. 2(b)), is generated by
considering two modes, $M_{1}$ and $M_{2}$, with different two-dimensional
normal distributions $N(\mu_{1},\Sigma_{1})$ and $N(\mu_{2},\Sigma_{2})$ and 9
intermediate modes. The parameter values of the modes $M_{1}$ and $M_{2}$ are:
$\Sigma_{1}=\begin{pmatrix}3.8418&-2.6474\\\ -2.6474&4.8478\end{pmatrix}$,
$\Sigma_{2}=\begin{pmatrix}1.5239&-0.5390\\\ -0.5390&1.6467\end{pmatrix}$, and
$\mu_{1}=(95,75)$ and $\mu_{2}=(5,5)$. $M_{1}$ is the initial mode, and
$M_{2}$ is the final mode. $M_{1}$ is transformed as follows. First, 500
samples $\left\\{k=1\ldots 500\right\\}$ are drawn from $M_{1}$. Sampling
continues as each individual value of the covariance matrix and the mean are
changed in 10 equal steps from their values in $M_{1}$ to those in $M_{2}$.
After the first step, 200 samples $\left\\{n=501\ldots 700\right\\}$ are taken
from the new normal distribution. After each new step 200 more samples are
added to the dataset. The final step ends at mode $M_{2}$. The squares show 1%
of the samples from each normal distribution, which are perturbed by uniform
noise from $[-10,\;10]$. A small level of noise is added to this dataset to
investigate how the algorithms react to noise.

The third synthetic data set, S3, is generated by drawing samples from two-
dimensional Gaussian distributions that rotate around a circle with 10 equal
shifts. Before each shift, 200 samples are generated using the current
Gaussian. In this data set, the noise (the blue cluster in the center of Fig.
2(c)) is generated by a Gaussian at the center of the circle. At each of 10
steps a random number of samples between 1 and 20 are removed from the outer
distribution and then this number of samples are drawn from and added to the
inner noise distribution. Table I summarizes the characteristics of the three
synthetic datasets.

Column 2 of Table I specifies the number of physically labeled subsets in each
data set. These subsets may, or may not, correspond to visually apparent or
computationally acquired clusters. For example, if you imagine Figure 2(b)
without the colors (which show the labels), most observers would assert that
this data set has only 1 cluster.

We also use two real-life datasets. The LG dataset is from a collection of
weather station nodes in the Le Genepi (LG) region in Switzerland [24]. Two
weeks of data at node 18 starting from October 10th 2007 are used in the
evaluation. We use average surface temperature (T) and humidity (H) readings
at node 18 over 10-minute intervals to form two dimensional input vectors
$\left\\{x_{n}=\left(T_{n},H_{n}\right)\right\\}$.

(a) Scatterplot of LG

(b) Time-series view of the data

Figure 3: Scatterplot and time-series views of the LG data. The scatterplot
shows three expected clusters.

Fig. 3(a) shows a scatter plot of the LG data. The imagery information from
the site shows that there is a snowy day during the two weeks of data
collection and the data confirms that a cold and windy day precedes the snow.
Fig. 3(b) shows this change of weather in time-series plots of temperature and
relative humidity data. Therefore, we assume that there are three physical
events, and that these may correspond to three clusters in the data: sunny
days before and after the snow, cold front moving in, and the snowy day and
label data accordingly, i.e., all of the points for days 1-6 and 10-14 are in
cluster 1 (blue), day 7 in cluster 2 (green), and days 8 and 9 in cluster 3
(red).

The second real-life data set is a gas sensor array under dynamic gas mixtures
from the UCI repository that contains conductivity samples at 100Hz obtained
from 16 chemical sensors (4 unique sensors). The conductivity of these types
of sensors changes in the presence of different gas mixture concentrations.
More information on the generated data set can be found in Fonollosa et al.
[25]. We select two pairs of each unique sensor (8 sensors) and take the mean
of their 100 samples per second over 5 minutes of the experiment as our
evaluation data set (a stream of 300 8-dimensional input vectors), which we
refer to as GSA. During this five minutes the sensors are exposed to two
different concentrations of gases CO and Ethylene. This will lead to three
distinct behaviours in the dataset: no gas being present, presence of CO and
presence of Ethylene.

###  6.2 Initialization

Macqueen’s sk-means algorithm has only one parameter to choose, $k$. We use
the recommended parameter values for OEC from [5]. These parameters are
forgetting factor $\lambda=0.9$, effective cluster boundary and outlier
boundary threshold of 0.99 and 0.999 respectively, a stabilization period of
$n_{s}=20$, and and an OEC forgetting factor 𝛌$\lambda_{\textit{OEC}}=0.9$.

The two clustering algorithms have slightly different initialization
procedures. In sk-means the first $k$ points are the initial cluster
prototypes and the index calculations start at the $k+1$-th point. In the OEC
clustering algorithm, the first $p+1$ points are used to calculate a single
cluster prototype and the cluster evaluation starts from point $p+2$ with only
one cluster in the system.

Upon start of the evaluation at step $n$, the iCVIs are initialized with
$C_{i,n}=0$, $M_{i,n}=n$ and $G_{i,n}=\overrightarrow{0}$ (zero vector in
$\Re^{p}$). After initialization, the clustering algorithms and iCVIs process
data one sample at a time.

##  7 Numerical Experiments

We first study the effect of forgetting in the indices. In these experiments
we specify $k$ to be the correct number of labels for sk-means for the
synthetic datasets. We show that both indices with forgetting reveal more
information about the data streams than when forgetting is not in use.

###  7.1 Forgetting or Not

In the iCVI indices (similar to their batch counterparts) each new data point
at time $n$ affects the overall value of the index with the weight $1/n$ (see
(25)). Therefore, over a long data stream, we expect $\textit{XB}(n)$ and
$\textit{DB}(n)$ to become saturated by data points and lose the sensitivity
they need to reflect changes due to new data inputs.

Fig. 4 shows the values of the two indices with and without forgetting for
online clustering in $S2$ with the OEC algorithm. The times when the
distribution of samples changes are indicated with red vertical lines. You can
see a sudden jump in the indices at these times. Even in a fairly small data
set such as $S2$, the jumps in the indices with no forgetting as $n$ increases
(bottom of Fig. 4) becomes very small. A data distribution change occurs
around 1750, but neither DB($n$) nor XB($n$) show this change. This suggests
that these two indices are not as suitable for monitoring evolving clusters in
streaming data as their forgetting factor counterparts, which are shown in the
top part of Fig. 4. The jumps in the indices with forgetting factor are larger
and much clearer than those in with no forgetting, and do not seem to depend
on the number of samples processed by the index.

Figure 4: Values of XB($n$), DB($n$), $\text{XB}_{\lambda}(n)$ and
$\text{DB}_{\lambda}(n)$ with $\lambda=0.9$ in OEC clusters of the S2 dataset.

(a) S1: Synthetic Locally Linear Processes

(b) S2: Shifting Gaussian Distribution

(c) S3: Circular Shifting Gaussian Distribution with Systematic Noise

Figure 5: Values of $\text{XB}_{\lambda}(n)$ and $\text{DB}_{\lambda}(n)$,
$\lambda=0.9$, for sk-means. The vertical lines (red) are times when a change
occurs in the dataset.

###  7.2 Interpretation of the iCVI model output over time

To study the usefulness of XB${}_{\lambda}(n)$ and DB${}_{\lambda}(n)$ in
interpreting the results of online clustering algorithms, we first plot their
values over time in clusters found by sk-means in the three synthetic
datasets. In most cases because of the larger number of clusters and lower
accuracy of sk-means in finding the expected clusters, the scale of the iCVI
values for the sk-means approach is much higher than the values corresponding
to OEC clusters. Therefore, we show the changes of $\text{XB}_{\lambda}(n)$
and $\text{DB}_{\lambda}(n)$ over time for only sk-means in Fig. 5. We will
discuss the noteworthy trends in OEC using separate figures.

In Fig. 5, The indices $\text{XB}_{\lambda}(n)$ and $\text{DB}_{\lambda}(n)$
have almost identical performance for data set $S1$, as can be seen by
comparing graphs of their values in Fig. 5(a). But their performance is quite
different for data sets $S2$ and $S3$, as can be seen in the top and bottom
views in Figs. 5(b) and 5(c), respectively. We highlight two main traits in
values of these indices. These traits relate to the appearance of new clusters
in the stream (spikes in the index) and performance of the clustering
algorithm in modeling the new cluster (the reduction of the index after a
spike).

####  7.2.1 Appearance of new clusters

CVIs assess cohesion and separation so sudden changes in the value of a CVI
indicate changes in the cohesion and separation of the clusters produced by
the clustering algorithm. Thus, a sudden change in an online validity index
usually indicates the appearance of a new cluster in the data stream. This
superficially relates to change-point detection in data streams but the
distinction needs to be made that iCVIs monitor the partitions produced by the
clustering algorithm and do not search for change-points in the streams.
However, change-points that result in appearance of new clusters in the data
streams change the cohesion of the clusters and cause spikes in the iCVI
values.

The vertical lines in Fig. 5 mark the times where new clusters appear in the
synthetic datasets. The sk-means algorithm performs reasonably well on the
$S1$ dataset, completely fails for the $S2$ dataset (Fig. 6(a)), and partially
identifies the evolving clusters in $S3$ (Fig. 6(b)). As we shall see, the OEC
algorithm performs reasonably at detecting times when new distributions are
created in all three datasets.

There is a jump in the values of both indices in Fig. 5 shortly after each new
cluster is introduced in the data stream. As the algorithms learn the
prototype that represents the new distribution in the data, the value of the
index drops (bear in mind that these indices are all min-optimal). This
behaviour is clearly shown in $\text{XB}_{\lambda}(n)$ plots but the
$\text{DB}_{\lambda}(n)$ plot in the top view of Fig. 5(c) does not reflect
this behaviour. We attribute this to partial identification of clusters in
$S3$ and the systematic noise in this dataset which affects
$\text{DB}_{\lambda}(n)$ more than $\text{XB}_{\lambda}(n)$.

####  7.2.2 Distress signals in the clustering algorithm

Another valuable asset of the graphs of incremental validity indices with
forgetting is that a streaming plot of their values can exhibit signs of
failure. The sk-means algorithm performs well in the $S1$ data set and we can
see a sharp peak (sudden increase and decrease) in Fig. 5(a) after the first
time the second cluster appears in the data. The $\text{XB}_{\lambda}$ and
$\text{DB}_{\lambda}$ indices are almost identical in this data set. However,
sk-means fails to identify the expected clusters in both $S2$ and $S3$ as
shown by the end-state partitions in Fig. 6 even when the correct number of
clusters are supplied to the algorithm. Let’s see how the
$\text{XB}_{\lambda}(n)$ plots reflect this problem in the sk-means clustering
algorithm.

Fig. 5(b) shows $\text{XB}_{\lambda}(n)$ over time and the final clusters in
$S2$ shown in Fig. 6(a) by sk-means. The $\text{XB}_{\lambda}(n)$ plot shows
that its values have an increasing trend. In this dataset, we have 11 clusters
that appear one by one in the stream. We know that the appearance of new
clusters results in spikes in the $\text{XB}_{\lambda}(n)$ values. After the
appearance of a new cluster in the data, we expect the clustering algorithm to
create a prototype for the cluster and update the partition to account for the
newly observed cluster. Subsequently, this should bring any min-optimal down
to similar (or lower) values that were observed before the new cluster
appeared. The total failure of sk-means in $S2$ is evident by its inability to
restore the $\text{XB}_{\lambda}(n)$ values after a new cluster enters the
data stream. Indeed, increasing values of these min-optimal indices signals
that things are going awry! The $\text{DB}_{\lambda}(n)$ values show smaller
peaks and decreases but still shows that sk-means cannot find good clusters to
sharply reduce the index. The graph in the top view of Fig. 4 shows the values
of both indices for the OEC algorithm applied to $S2$, which identifies all
the clusters correctly. The sharp peaks in the values verifies the fact that
the clustering algorithm finds all the expected clusters.

(a) S2: Shifting Gaussian Distribution

(b) S3: Circular Shifting Gaussian Distribution with Systematic Noise

Figure 6: The terminal clusters produced by sk-means in $S2$ and $S3$

Figs. 5(c) and 6(b) tell a slightly different story for $S3$. The sk-means
clusters in this dataset provide some separation between the expected
clusters. The separation is achieved by creating two partitions with 3
clusters each, a partition with two clusters and the final cluster is
identified correctly. At two points during the experiment, there are
significant reductions in the $\text{XB}_{\lambda}(n)$ values (around 800 and
1500 samples) which correspond to the times after the two subsets of 3
clusters are created. The decreasing trend of $\text{XB}_{\lambda}(n)$ towards
the end is due to the fact that sk-means performs slightly better in
identifying the last 3 clusters. In particular, correct identification of the
last cluster reduces the index value quite sharply. The
$\text{DB}_{\lambda}(n)$ values in this dataset are very different to the
$\text{XB}_{\lambda}(n)$ values. The index $\text{DB}_{\lambda}(n)$ generates
large peaks around the noisy data and shows higher sensitivity to systematic
noise in the data.

Sharp peaks corresponding to sudden drops in $\text{XB}_{\lambda}(n)$ after
the spikes (appearance of a new cluster) relates to the clustering algorithm
appropriately creating a prototype for the cluster, while a gradual decrease
after a spike is a sign of the failure of the algorithm to identify the new
cluster. The $\text{DB}_{\lambda}(n)$ index is more sensitive to large
amplitude local noise and generates smaller peaks for the new clusters. This
hinders the interpretation of the clustering results in terms of appearance of
new clusters and learning problems in the algorithm.

Figure 7: Terminal OEC clusters on $S3$ dataset.

So far we only looked at problems in the sk-means algorithm because OEC finds
all the expected clusters in these three synthetic datasets. However, In $S3$,
as shown in Fig. 7, one of the cluster prototypes has been expanded and covers
the systematic noise in the centre of the plot (The expanded prototype is a
single prototype that represents the two clusters of data captured by the
horizontally elongated ellipse in Fig. 7). Fig. 8 shows the values of the two
indices for OEC clustering in this dataset. Both indices have similar trends
with $\text{DB}_{\lambda}(n)$ producing larger peaks for noisy data. Neither
graph in Fig. 8 experiences a drop after creating this prototype. This
indicates that OEC does not produce a good model of the data which matches the
visual assessment of the clusters.

The overall conclusion from our synthetic data experiments is that
$\text{XB}_{\lambda}(n)$ is much more effective than $\text{DB}_{\lambda}(n)$,
so the remaining discussion will involve only $\text{XB}_{\lambda}(n)$.

Figure 8: $\text{XB}_{\lambda}(n)$ and $\text{DB}_{\lambda}(n)$ values for OEC
clusters in $S3$

###  7.3 Experiments with real data

Since we have only assumed information about the ground truth in the real
datasets it is harder to interpret the $\text{XB}_{\lambda}(n)$ values. We
will use the insights obtained from our analysis of the synthetic datasets to
understand the behaviour of $\text{XB}_{\lambda}(n)$ when clustering in the
real-life datasets. Fig. 9(a) shows values of $\text{XB}_{\lambda}(n)$
calculated for both OEC and sk-means (with $k=3$ clusters) in the LG dataset.
The major event in this dataset is a snowfall event (the approximate location
of the event between about $n=700$ and $n=900$ is marked in the figure). OEC
identifies three clusters in this dataset [5] which correspond to normal days,
high wind before the snow and the snow. The two very close sharp spikes
followed by drops in both indices values corresponds the spike for the wind
before the snow and the snow fall as identified by OEC. The first spike
corresponds to the time when the high wind had become the dominant feature in
the area. Some level of wind can be observed from step 600.

The index values for sk-means clusters do not correspond to the major events
in the dataset. There is only one major peak around the very start of the high
winds at about step 600. Since sk-means does not consider the correlation in
time between the data points, it fails to account for the temporal nature of
the clusters in the LG data. As shown in Fig. 3(b), we must to look at the
evolution of data in time to see the major events and the scatter plot of the
data does not show a clear cluster tendency. The large overlap between
clusters of data corresponding to major events leads to the fact that aligning
clusters with the events does not guarantee the maximal separation. This is
reflected in the index values for sk-means in Fig. 9(a) that shows no clear
increasing trend or particularly slow reduction in the index values. However,
the index values for sk-means clusters are mainly lower than those for OEC
clusters showing that sk-means clustering creates relatively better clusters
in terms of cohesion and separation.

(a) LG: Sensorscope project in Le Genepi

(b) GSA: Gas Sensor Array. Vertical lines indicate the times where
concentration of gases change

Figure 9: The $\text{XB}_{\lambda}(n)$ index over the real-life datasets. The
y-axis of GSA plot is shown in logarithmic scale.

The GSA dataset is collected in a more controlled environment than LG and
hence, has a more recognizable cluster tendency. In this dataset two different
gas concentrations are introduced during the experiment. When the gas is
introduced (times indicated by the red vertical lines) there is a delay until
the sensors react to the presence of the gas, which is seen as a delay of the
values in Fig. 9(b). The peaks in both indices occur a few seconds after the
introduction of the gas. For the first event at about $n=75$, a gas is
introduced when no other gas is present, while the second event at about
$n=180$ corresponds to removing one gas and introducing another gas. In the
second case, the indices show two peaks close to each other for both
clustering algorithms. The last event at about $n=270$ corresponds to emptying
the chamber and since the empty chamber had been seen by the clustering
algorithms, this results in a much smaller peak in the graphs. OEC has smaller
$\text{XB}_{\lambda}(n)$ values than sk-means, which indicate that in this
dataset, OEC results in a better separation between the expected clusters than
sk-means.

##  8 Conclusions

To the best of our knowledge, the concept of incremental validity indices
(iCVIs) in data stream clustering has not been studied before. Our iCVIs offer
an unsupervised validation mechanism for online clustering algorithms. In this
article, we introduced the novel concept of online iCVIs and derived
forgetting and non-forgetting versions of two well-known internal validity
indices (the Xie-Beni indices (XB($n$) and $\text{XB}_{\lambda}(n)$) and the
Davies-Bouldin indices (DB($n$) and $\text{DB}_{\lambda}(n)$)). Our
experiments used two different styles of algorithms for streaming clustering:
sk-means, which does not account for the historical context of streaming data;
and OEC, which retains a history of time dependency in clusters through the
retention of its cluster statistics (means and covariances). We discussed how
these indices can be used to interpret the performance of online clustering
algorithms with respect to the appearance of new clusters and how the
clustering algorithm reacts to evolving clusters. An increasing trend in any
min-optimal index seems to point to a more questionable belief about cluster
quality in the clustering of the data. Thus, our iCVI models afford a means
for sending distress signals about evolving clusters to real time monitors.

Our experiments indicate that the most reliable of the four incremental
indices {XB($n$),$\text{XB}_{\lambda}(n)$, DB($n$),$\text{DB}_{\lambda}(n)$}
is $\text{XB}_{\lambda}(n)$. And of the two clustering algorithms used, OEC
seems to perform much better than sk-means. But definitive conclusions require
many more tests. There are many, many internal CVIs. Our next focus will be on
deriving other incremental validity indices and conducting comparative studies
among different types of indices.

##  9 Acknowledgment

This research was supported under Australian Research Council’s Discovery
Projects funding scheme (project number DE150100104).

## References

  * [1] S. Guha, A. Meyerson, N. Mishra, R. Motwani, and L. O’Callaghan, “Clustering data streams: Theory and practice,” _IEEE Transaction on Knowledge and Data Engineering_ , vol. 15, no. 3, pp. 515–528, 2003.
  * [2] J. A. Silva, E. R. Faria, R. C. Barros, E. R. Hruschka, A. C. P. L. F. de Carvalho, and J. ao Gama, “Data stream clustering: A survey,” _ACM Computing Surveys_ , vol. 46, no. 1, pp. 13–31, Jul. 2013.
  * [3] M. Ackerman and S. Dasgupta, “Incremental clustering: The case for extra clusters,” in _Neural Information Processing Systems (NIPS)_ , Montreal, Canada, December 2014.
  * [4] P. Angelov, _Evolving Takagi-Sugeno Fuzzy Systems from Streaming Data (eTS+)_. John Wiley & Sons, Inc., 2010, pp. 21–50.
  * [5] M. Moshtaghi, J. Bezdek, and C. Leckie, “Online clustering of multivariate time-series,” in _Proceedings of SIAM Conference on Data Mining_ , Florida, USA, May 2016.
  * [6] O. Arbelaitz, I. Gurrutxaga, J. Muguerza, J. M. Perez, and I. Perona, “An extensive comparative study of cluster validity indices,” _Pattern Recognition_ , vol. 46, no. 1, pp. 243 – 256, 2013.
  * [7] X. L. Xie and G. Beni, “A validity measure for fuzzy clustering,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 13, no. 8, pp. 841–847, Aug 1991.
  * [8] D. L. Davies and D. W. Bouldin, “A cluster separation measure,” _IEEE Transaction on Pattern Analysis and Machine Intelligence_ , vol. 1, no. 2, pp. 224–227, Feb. 1979.
  * [9] C. C. Aggarwal, J. Han, J. Wang, and P. S. Yu, _On Clustering Massive Data Streams: A Summarization Paradigm_. Springer US, 2007, vol. 31, pp. 9–38.
  * [10] N. Ailon, R. Jaiswal, and C. Monteleoni, “Streaming k-means approximation,” in _Neural Information Processing Systems 2009_ , vol. 22, 2009, pp. 10–18.
  * [11] M. Salehi, C. Leckie, M. Moshtaghi, and T. Vaithianathan, “A relevance weighted ensemble model for anomaly detection in switching data streams,” in _Proceedings of PAKDD_ , 2014, pp. 461–473.
  * [12] P. Angelov and X. Zhou, “Evolving fuzzy-rule-based classifiers from data streams,” _IEEE Transactions on Fuzzy Systems_ , vol. 16, no. 6, pp. 1462–1475, 2008.
  * [13] F. Cao, M. Ester, W. Qian, and A. Zhou, “Density-based clustering over an evolving data stream with noise,” in _SIAM Conf. on Data Mining_ , 2006, pp. 328–339.
  * [14] P. Kranen, I. Assent, C. Baldauf, and T. Seid, “The clustree: Indexing micro-clusters for anytime stream mining,” _Knowledge and Information Systems_ , vol. 29, no. 2, pp. 249–272, 2011.
  * [15] N. Mozafari, S. Hashemi, and A. Hamzeh, “A statistical approach for clustering in streaming data,” _Artificial Intelligence Research_ , vol. 3, pp. 38–45, 2014.
  * [16] G. W. Milligan and M. C. Cooper, “An examination of procedures for determining the number of clusters in a data set,” _Psychometrika_ , vol. 50, no. 2, pp. 159–179, 1985.
  * [17] K. Schermelleh-Engel, M. Moosbrugger, and H. Helfried, “Evaluating the fit of structural equation models: Tests of significance and descriptive goodness-of-fit measures,” _Methods of Psychological Research_ , vol. 8, pp. 23 – 74, 2015.
  * [18] J. C. Bezdek, _Pattern Recognition with Fuzzy Objective Function Algorithms_. Norwell, MA, USA: Kluwer Academic Publishers, 1981.
  * [19] N. R. Pal, J. C. Bezdek, and R. J. Hathaway, “Sequential competitive learning and the fuzzy c-Means clustering algorithms,” _Neural Networks_ , vol. 9, no. 5, pp. 787 – 796, 1996.
  * [20] J. MacQueen, “Some methods for classification and analysis of multivariate observations,” in _In 5th Berkeley Symposium on Mathematical Statistics and Probability_ , 1967, pp. 281–297.
  * [21] J. C. Bezdek, _A Primer on Cluster Analysis: Four Basic Methods that (Usually) Work_. Sarasota, FL., USA: First Design Publishing, 2017.
  * [22] N. R. Pal and J. C. Bezdek, “On cluster validity for the fuzzy c-means model,” _IEEE Transactions on Fuzzy Systems_ , vol. 3, no. 3, pp. 370–379, 1995.
  * [23] S. Araki, H. Nomura, and N. Wakami, “Segmentation of thermal images using the fuzzy c-means algorithm,” in _Proceedings of the Second IEEE International Conference on Fuzzy Systems_ , San Francisco, USA, April 1993, pp. 719–724.
  * [24] “SensorScope. http://lcav.epfl.ch/page-86035-en.html,” 2007. [Online]. Available: http://lcav.epfl.ch/page-86035-en.html
  * [25] J. Fonollosa, S. Sheik, R. Huerta, and S. Marco, “Reservoir computing compensates slow response of chemosensor arrays exposed to fast varying gas concentrations in continuous monitoring,” _Sensors and Actuators B: Chemical_ , vol. 215, pp. 618 – 629, 2015.


L_{\textbf{l}}(m)|L_{\textbf{l}}(m)\rangle=1$ (see Fig. 17 (c)):

$\displaystyle\begin{aligned} &k=(N-M)/2+1,...,N/2\\\
&\psi_{2k-1}(0)=-i\psi_{2k}(0),\quad\psi_{2k-1}(\beta_{-})=i\psi_{2k}(\beta_{-}),\\\
&\psi_{2k-1}(\beta_{+})=-i\psi_{2k}(\beta_{+}),\quad\psi_{2k-1}(2\beta)=i\psi_{2k}(2\beta).\end{aligned}$
(148)

(a) $1\leq i\leq N-M$ (b) $N-M<i\leq M$ (c) $M<i\leq N$

Figure 18: Imaginary time contour $\tau\in(0,2\beta)$ for the Euclidean path
integral computing the numerator of equation (143) when $M>N/2$. The twist
operator inserted in the left system is indicated by the dashed line. The
insertion of a measurement operator in the left and right sides is indicated
by the triangles. When the projective measurement operator is inserted in the
left side, the effect of the twist operator is trivial and we therefore omit
the dashed line. (a) Boundary conditions (149) for Majorana fermions
$\psi_{i}$, $1\leq i\leq N-M$. (b) Boundary conditions (150) for Majorana
fermions $\psi_{i}$, $N-M<i\leq M$. (c) Boundary conditions (151) for Majorana
fermions $\psi_{i}$, $M<i\leq N$.

####  M>N/2

Increasing the number of measured Majorana fermions such that $M>N/2$, the
measurement prescription in the right side and left side is the same as
before, but now the measurement occurs in both sides for the Majorana fermions
$\psi_{i}$ with $N-M<i\leq M$.181818For simplicity, we postselect the outcomes
of the right-side measurement and the left-side measurement to be the same so
that the method we developed in Section 2 and Appendix A can be directly
applied. If the measurement outcomes are different, we cannot simply introduce
a new Majorana fermion $\chi_{k}$, because it would not satisfy the
conventional anti-periodic boundary conditions for fermionic fields. A
graphical representation of the resulting boundary conditions is given in Fig.
18.

For $1\leq i\leq N-M$, the boundary conditions must account for the
measurement operator inserted at $\tau=\beta/2,3\beta/2$ as well as the twist
operator inserted in $\tau=0,\beta,2\beta$ (see Fig. 18 (a)):

$\displaystyle\begin{aligned} &k=1,...,(N-M)/2\\\
&\psi_{2k-1}(\frac{\beta_{-}}{2})=i\psi_{2k}(\frac{\beta_{-}}{2}),\quad\psi_{2k-1}(\frac{\beta_{+}}{2})=-i\psi_{2k}(\frac{\beta_{+}}{2}),\\\
&\psi_{2k-1}(\frac{3\beta_{-}}{2})=i\psi_{2k}(\frac{3\beta_{-}}{2}),\quad\psi_{2k-1}(\frac{3\beta_{+}}{2})=-i\psi_{2k}(\frac{3\beta_{+}}{2}),\\\
&\psi_{2k-1}(0)=-\psi_{2k-1}(2\beta),\quad\psi_{2k}(0)=-\psi_{2k}(2\beta).\end{aligned}$
(149)

For $N-M<i\leq M$, there are measurement operators inserted in both
sides—namely at $\tau=0,\beta/2,\beta,3\beta/2,2\beta$—and the action of the
twist operator is trivial. This leads to the boundary conditions (see Fig. 18
(b))

$\displaystyle\begin{aligned} &k=(N-M)/2+1,...,M/2\\\
&\psi_{2k-1}(0)=-i\psi_{2k}(0),\quad\psi_{2k-1}(\beta_{-})=i\psi_{2k}(\beta_{-}),\\\
&\psi_{2k-1}(\beta_{+})=-i\psi_{2k}(\beta_{+}),\quad\psi_{2k-1}(2\beta)=i\psi_{2k}(2\beta),\\\
&\psi_{2k-1}(\frac{\beta_{-}}{2})=i\psi_{2k}(\frac{\beta_{-}}{2}),\quad\psi_{2k-1}(\frac{\beta_{+}}{2})=-i\psi_{2k}(\frac{\beta_{+}}{2}),\\\
&\psi_{2k-1}(\frac{3\beta_{-}}{2})=i\psi_{2k}(\frac{3\beta_{-}}{2}),\quad\psi_{2k-1}(\frac{3\beta_{+}}{2})=-i\psi_{2k}(\frac{3\beta_{+}}{2}).\end{aligned}$
(150)

Finally, for $M<i\leq N$, the boundary conditions must take into account the
measurement operator inserted at $\tau=0,\beta,2\beta$, while the twist
operator at the same point is trivial (see Fig. 18 (c)):

$\displaystyle\begin{aligned} &k=M/2+1,...,N/2\\\
&\psi_{2k-1}(0)=-i\psi_{2k}(0),\quad\psi_{2k-1}(\beta_{-})=i\psi_{2k}(\beta_{-}),\\\
&\psi_{2k-1}(\beta_{+})=-i\psi_{2k}(\beta_{+}),\quad\psi_{2k-1}(2\beta)=i\psi_{2k}(2\beta).\end{aligned}$
(151)

The denominator of equation (143) can also represented by a path integral with
the same action (145). The difference is in the boundary conditions, because
there is no twist operator. More explicitly, the boundary conditions for the
denominator are given by equations (146)-(151) after replacing the last line
in equation (146), equation (147), and the last line in equation (149) by

$\displaystyle\psi_{i}(0)=-\psi_{i}(\beta_{-}),\quad\psi_{i}(\beta_{+})=-\psi_{i}(2\beta),\quad
i=1,...,N-M.$ (152)

Having identifying the boundary conditions, we can now proceed to introduce
bilocal fields, derive and numerically solve the Schwinger-Dyson equations,
evaluate the numerator and denominator on-shell actions $I_{num}(m)$ and
$I_{den}(m)$, and finally compute the Renyi-2 mutual information—which is
given in the saddle point approximation by
$I_{LR}^{(2)}(\Phi(m))=2(I_{num}(m)-I_{den}(m))$. This analysis is completely
analogous to the one carried out for the one-sided measurement in Section 2.3
and Appendix A, and we will therefore omit it here.

Our numerical results for the Renyi-2 mutual information
$\mathcal{I}^{(2)}_{LR}(\Phi(m))=I_{LR}^{(2)}(\Phi(m))/N$ are shown in Fig.
19. As in the one-sided measurement case, the mutual information remains
finite for any $m<1$ and only vanishes for $m=1$. In other words, the left and
right systems are completely disentangled only when all Majorana fermions are
measured. One way to understand this result is to remember that at infinite
temperature the TFD state reduces to EPR pairs between the left and right
Majorana fermions. In this case the mutual information is

$\displaystyle\mathcal{I}^{\Phi}_{LR}(m)=\begin{cases}(1-2m)\log 2,&m<1/2,\\\
0,&m\geq 1/2.\end{cases}$ (153)

The transition between the entangled phase and the completely disentangled
phase occurs at $m=1/2$. This can be readily understood if we recall that for
$m=1/2$ (which is $M=N/2$) one fermion is measured for each one of the $N$ EPR
pairs. In particular, the right fermion of each EPR pair is measured for
$i=1,...,N/2$, while the left fermion of each EPR pair is measured for
$i=N/2+1,...,N$. Now imagine we turn on a small imaginary time evolution, i.e.
we consider the TFD state at high but finite temperature. Then the all-to-all
interaction in the Hamiltonian will generate entanglement between the two
sides even for $m>1/2$. This is particularly evident for example in the
$\beta\mathcal{J}=2$ numerical result in Fig. 19.

Finally, we would like to remark that the holographic dual analysis of the
two-sided measurement setup described in the present Appendix is complicated
by the fact that when large subsets of fermions are measured on both sides, it
is not clear whether a geometric bulk dual description still exists.
Therefore, when $m$ is small we expect the bulk dual analysis to remain
substantially similar to the one described in Section 3 for the one-sided
measurement case. But when $m$ is large, it is not clear what bulk geometry is
left behind by the measurement, if any. For intermediate values of $m$, the
transition from a semiclassical geometry to a non-geometric description could
be radical and it is not well-understood. We thus leave the bulk dual analysis
of the two-sided measurement case for future work.

Figure 19: Renyi-2 mutual information $\mathcal{I}_{LR}^{(2)}(\Phi(m))$
between unmeasured Majorana fermions on the left and right sides for the two-
sided measurement case. $m=M/N$ is the ratio of measured Majorana fermions on
the left and right sides. The mutual information is always non-vanishing for
$m<1$, and vanishes only when all the Majorana fermions are measured ($m=1$).

##  Appendix D Entanglement entropy in JT gravity with matter CFT

Varying the dilaton and metric in the JT plus matter action (equation (47)),
respectively, yield the equations of motion

$\displaystyle R=-2,$ (154) $\displaystyle\frac{1}{8\pi
G_{N}}(\nabla_{\mu}\nabla_{\nu}\phi-
g_{\mu\nu}\nabla^{2}\phi+g_{\mu\nu}\phi)=T_{\mu\nu},$ (155)

where $T_{\mu\nu}$ is the stress tensor of the matter CFT. We first consider
the solution at zero temperature, setting $T_{\mu\nu}$ to zero. The metric in
Euclidean Poincaré coordinates is

$\displaystyle ds^{2}=\frac{dt^{2}+dz^{2}}{z^{2}},$ (156)

and we implement the following boundary conditions at $z=\epsilon$
maldacena2016conformal ,

$\displaystyle
g|_{\text{bdy}}=\frac{1}{\epsilon^{2}},\quad\phi(\epsilon)=\frac{1}{\epsilon}.$
(157)

The solution for the dilaton is then given by

$\displaystyle\phi=\frac{\alpha+\gamma t+\delta(t^{2}+z^{2})}{z},$ (158)

where $\alpha,\gamma,\delta$ are constants determined by the boundary
conditions, $\alpha=1$, $\gamma=\delta=0$.

To describe a thermal state with temperature $1/\beta$, we can make a series
of coordinate transformations. To this end, we first introduce $x=z+it$,
$\bar{x}=z-it$ with $ds^{2}=\frac{4dxd\bar{x}}{(x+\bar{x})^{2}}$. Now then
take

$\displaystyle
x=\frac{\beta}{2\pi}\tanh\frac{\pi}{\beta}y,\quad\bar{x}=\frac{\beta}{2\pi}\tanh\frac{\pi}{\beta}\bar{y}.$
(159)

Finally, with $y=\sigma+i\tau$, $\bar{y}=\sigma-i\tau$, the metric becomes

$\displaystyle
ds^{2}=\frac{4\pi^{2}}{\beta^{2}}\frac{d\sigma^{2}+d\tau^{2}}{\sinh^{2}\frac{2\pi}{\beta}\sigma}.$
(160)

Implementing the boundary condition at $\sigma=\epsilon$,
$g_{\tau\tau}=\frac{1}{\epsilon^{2}}$,
$\phi(\epsilon)=\frac{\phi_{r}}{\epsilon}$, yields the dilaton solution

$\displaystyle\phi(\sigma)=\phi_{r}\frac{2\pi}{\beta}\frac{1}{\tanh\frac{2\pi}{\beta}\sigma}.$
(161)

Figure 20: Two copies of Poincare disk with one of them filling outside
$|z|>1$.

We can now turn to understanding the bulk matter entropy, computed via a
“doubling trick.” We introduce a second Euclidean copy of the spacetime, glued
together as in Fig. 20. To this end, we employ a (final) coordinate
transformation,

$\displaystyle\begin{aligned} &\text{first copy}\\\
&z=e^{\frac{2\pi}{\beta}(\sigma+i\tau)},\quad\bar{z}=e^{\frac{2\pi}{\beta}(\sigma-i\tau)},\quad
ds^{2}=\frac{4dzd\bar{z}}{(1-|z|^{2})^{2}},\quad|z|<1,\\\ &\text{second
copy}\\\
&z=e^{-\frac{2\pi}{\beta}(\sigma+i\tau)},\quad\bar{z}=e^{-\frac{2\pi}{\beta}(\sigma-i\tau)},\quad
ds^{2}=\frac{4dzd\bar{z}}{(1-|z|^{2})^{2}},\quad|z|>1.\end{aligned}$ (162)

The metric is identical in the two copies, except one copy is mapped to
$|z|<1$, and the other one to $|z|>1$. In other words, the two copies are
related by the transformation $z\rightarrow\frac{1}{z}$. The free Dirac
fermion can propagate in the entire plane, as shown in Fig. 20. Note that the
UV cutoff of the two copies is located at
$|z|=1\pm\frac{2\pi\epsilon}{\beta}$, respectively.

Since the metric differs from a flat space by only a Weyl factor
$ds^{2}=\frac{\sum_{i}dx^{i}dx^{i}}{\Omega(x)^{2}}$—and as noted above, the
fermions propagate in the entire plane—we know the entanglement entropy
between two twist operators at $\bf x_{1},x_{2}$ near the vacuum state is
given by

$\displaystyle S_{\text{CFT}}({\bf x}_{1},{\bf
x}_{2})=\frac{c}{6}\log\frac{{\bf x_{12}^{2}}}{\Omega({\bf x_{1}})\Omega({\bf
x_{2}})},$ (163)

where ${\bf x}_{12}={\bf x}_{1}-{\bf x}_{2}$. Finally, the dilaton field in
this coordinate is given by

$\displaystyle\phi(z,\bar{z})=\phi_{r}\frac{2\pi}{\beta}\frac{1+|z|^{2}}{|1-|z|^{2}|}.$
(164)

##  Appendix E Haar random unitary model for TFDs with measurement

In this appendix, we construct a simple model that shows the essential physics
of measurement-induced entanglement wedge phase transitions in the context of
bulk reconstruction. Our simple model consists of two random unitaries $U_{L}$
and $U_{R}$ with equal dimension $2^{N}$, drawn from a Haar distribution. To
allow for a more general setup, we take the bulk state of the system to be
$|\psi\rangle$. As shown in Fig. 21, these random unitaries map the state
$|\psi\rangle$ from bulk to boundary (we take the incoming legs at the bottom
of each unitary to be the bulk Hilbert space, and the outgoing legs on top of
each unitary to be the boundary Hilbert space, which we refer to as the left
and right sides for simplicity).

Figure 21: A model with two random unitaries, $U_{L}$ and $U_{R}$,
representing the two sides each with $N$ qubits. $|\psi\rangle$ denotes a bulk
state. $M$ out of $N$ qubits in the left side are projected to state
$|0\rangle$. $W$ denotes an operator acting on the left bulk Hilbert space.

In analogy with the measurements considered in the main text, we measure $M$
qubits in the left boundary Hilbert space. In the main text, we found a phase
transition occurs when the fraction $m=M/N$ of measured fermions exceeds a
critical value $m^{*}$ (see e.g. Fig. 1). When $m>m^{*}$, the information that
was stored in the left side before measurement is teleported to the right
side. Here, in this simple model, we ask a similar question. Suppose we have
an operator $W$ in the left bulk Hilbert space acting on the state
$|\psi\rangle$. How many qubits $M$ need to be measured in the left side such
that this operator can be successfully reconstructed in the right boundary
Hilbert space? In Fig. 21, $M$ out of $N$ qubits in the left side are
projected to state $|0\rangle$. The measurement outcome is not important here
since we are interested in the typical behavior of random unitaries drawn from
a Haar distribution.

To answer this question, we employ the technique introduced in Ref.
akers2022black ; Antonini:2022sfm  to give a bound on the reconstruction
error. This calculation is straightforward, and we refer to Ref.
Antonini:2022sfm  for details. Here we just present the final result. Let
$W_{R}$ be an operator acting on the right boundary Hilbert space, and
$V=\frac{1}{2^{M/2}}U_{L}\otimes U_{R}$ be the linear map from the bulk to the
boundary for both sides (here the prefactor comes from measurements). The
reconstruction error is defined by

$\displaystyle||VW|\psi\rangle-W_{R}V|\psi\rangle||\leq\epsilon_{1},$ (165)

where $||.||$ denotes the norm of the state. The reconstruction error for a
typical error modeled by a random unitary drawn from the Haar distribution is
bounded by $\epsilon_{1}\leq\sqrt{\epsilon_{2}}$, where

$\displaystyle\epsilon_{2}\leq 2\sqrt{\frac{2^{N-M}}{e^{S^{(2)}_{L}(\Psi)}}}.$
(166)

$S^{(2)}_{L}(\Psi)=\text{Tr}_{R}\left[\text{Tr}_{L}(\Psi)^{2}\right]$ is the
pre-measurement Renyi-2 entanglement entropy between left and right side
evaluated for the state $\Psi=|\psi\rangle\langle\psi|$. This bound
intuitively states that the reconstruction would be successful if the
measurement is implemented on a subset of fermions large enough that the
entanglement of the remaining $N-M$ unmeasured qubits is much smaller than the
pre-measurement entanglement between the two sides, i.e. the entanglement
resource. With this understanding, we can give an estimate of the critical
value of $m$ by solving $\frac{2^{N-M}}{e^{S^{(2)}_{L}(\Psi)}}=1$. The
transition point is then given by

$\displaystyle m^{*}=1-\frac{S^{(2)}_{L}(\Psi)}{N\log 2}.$ (167)

Roughly speaking, this means that the critical point is determined by the pre-
measurement entanglement between the two sides: the greater the entanglement,
the fewer number of qubits needs to be measured to achieve the transition.

##  Appendix F Bound on teleportation fidelity

In this appendix, we are going to relate the bound on teleportation fidelity
to the left-right correlation functions studied in Section 4. Our discussion
is similar to that of Ref. schuster2022many . We will use diagrams to
represent states and amplitudes. Assuming two copies of the same Hilbert
space, each consisting of $N/2$ qubits, we define some relevant states and
projection operators as follows,

$\displaystyle\begin{aligned} {\begin{array}[]{c}\vbox{\vskip
1.19998pt\hbox{\hskip
0.79999pt\includegraphics[scale={0.4}]{Figures/fig_TFD_1.png}}}\end{array}}&=\frac{1}{(\text{Tr}[e^{-\beta
H}])^{1/2}}\sum_{n=1}^{N/2}e^{-\beta
E_{n}}|E_{n}\rangle\otimes|E_{n}\rangle,\\\ {\begin{array}[]{c}\vbox{\vskip
1.19998pt\hbox{\hskip
0.79999pt\includegraphics[scale={0.4}]{Figures/fig_EPR.png}}}\end{array}}&=\frac{1}{d^{1/2}}\sum_{i=1}^{d}|i\rangle\otimes|i\rangle,\\\
{\begin{array}[]{c}\vbox{\vskip 1.05002pt\hbox{\hskip
0.70001pt\includegraphics[scale={0.35}]{Figures/fig_projection.png}}}\end{array}}&=\bigotimes_{k=1}^{N/2}|s_{k}\rangle\langle
s_{k}|,\end{aligned}$ (168)

where $H$ denotes the Hamiltonian and $|E_{n}\rangle$ denotes the energy
eigenstates. $|s_{k}\rangle,s_{k}=\pm$ are two possible states of a qubit, and
$|s_{k}\rangle\langle s_{k}|$ denotes a projection operator. In this appendix,
we use $s_{k}$ to denote the measurement outcome (note the change in notation
from Section 4, where we used $l_{k}$ to denote the measurement outcome).

The state teleportation protocol conditioned on measurement outcome $s_{k}$ is
given by

$\displaystyle{\begin{array}[]{c}\vbox{\vskip 1.19998pt\hbox{\hskip
0.79999pt\includegraphics[scale={0.4}]{Figures/fig_teleport.png}}}\end{array}}$
(170)

where $U$ is a time evolution operator and $U_{s}$ is a unitary operator
dependent on the measurement outcome, which we refer to as a decoding
operator. $D$ is another decoding operator at the end of the teleportation.
Thus, the teleportation protocol consists of three parts:

  1. 1.

At time $t=-t_{0}$, a state $|\psi\rangle$ is inserted in the left side;

  2. 2.

At time $t=0$, a projective measurement is performed in the left side and a
decoding operator is implemented in the right side according to the
measurement outcome;

  3. 3.

At time $t=t_{0}$, another decoding operator $D$ is performed in the right
side, and the state is teleported to the right side.

At time $t=0$, the measurement-decoding protocol should be understood as a
quantum channel, which can be defined by a Kraus operator,

$\displaystyle K_{\textbf{s}}=\left(\otimes_{k}|s_{k}\rangle\langle
s_{k}|\right)\otimes
U_{\textbf{s}},\quad\sum_{\textbf{s}}K_{\textbf{s}}^{\dagger}K_{\textbf{s}}=1,$
(171)

where the decoding operator we perform on the right side depends on the
measurement outcome on the left. Conditioned on the measurement outcome, the
full protocol is given in (170).

The state teleportation fidelity is related to the fidelity of the
distillation of EPR pairs yoshida2019disentangling ,

$\displaystyle{\begin{array}[]{c}\vbox{\vskip 1.10999pt\hbox{\hskip
0.73999pt\includegraphics[scale={0.37}]{Figures/fig_distill.png}}}\end{array}}$
(173)

Here, instead of inserting a state, a swap between two qubits is performed in
the right side, and finally a distillation of EPR states between the swap
qubit is implemented. The average fidelity of state teleportation is related
to the fidelity of EPR distillation by,

$\displaystyle\langle F_{\psi}\rangle=\frac{d}{d+1}F_{EPR}+\frac{1}{d+1},$
(174)

where $d$ is the dimension of the state Hilbert space, and
$\langle\cdot\rangle$ implies an average over all states.

The fidelity for distillation of an EPR state is

$\displaystyle
F_{EPR}=\sum_{\textbf{s}_{1},\textbf{s}_{2}}{\begin{array}[]{c}\vbox{\vskip
1.80002pt\hbox{\hskip
1.20001pt\includegraphics[scale={0.6}]{Figures/fig_F_EPR.png}}}\end{array}}=\frac{1}{d^{4}}\sum_{\textbf{s}_{1},\textbf{s}_{2}}\sum_{Q_{1},Q_{2}}{\begin{array}[]{c}\vbox{\vskip
1.80002pt\hbox{\hskip
1.20001pt\includegraphics[scale={0.6}]{Figures/fig_F_EPR_operator.png}}}\end{array}}$
(177)

where in the first diagram, the projection operator enforces
$\textbf{s}_{1}=\textbf{s}_{2}$. The sum over all possible measurement
outcomes occurs because the full fidelity is given by the sum of fidelities
conditioned on each measurement outcome. To get the second diagram, the swap
operator is replaced by the following identity,

$\displaystyle{\begin{array}[]{c}\vbox{\vskip 1.5pt\hbox{\hskip
1.0pt\includegraphics[scale={0.5}]{Figures/fig_swap.png}}}\end{array}}=\frac{1}{d}\sum_{Q}{\begin{array}[]{c}\vbox{\vskip
1.5pt\hbox{\hskip
1.0pt\includegraphics[scale={0.5}]{Figures/fig_swap_operator.png}}}\end{array}}$
(180)

where $Q$ is a basis operator living in the $d$ dimensional Hilbert space. For
instance for $d=2$, $Q=\mathds{1},X,Y,Z$ the identity and three Pauli
matrices. Then we move one of the operators to the right side, and define
$\tilde{Q}=D^{\dagger}QD$. In the case of a qubit, the decoding operator is
$D=Y$ brown2019quantum ; nezami2021quantum ; schuster2022many . For the SYK
model we consider in the main text, the decoding operator is trivial
$D=\mathds{1}$, and so we neglect it in the main text.

Thus, the distillation fidelity is lower bounded by

$\displaystyle
F_{EPR}=\frac{1}{d^{4}}\sum_{\textbf{s}_{1},\textbf{s}_{2}}\sum_{Q_{1},Q_{2}}{\begin{array}[]{c}\vbox{\vskip
1.85999pt\hbox{\hskip
1.23999pt\includegraphics[scale={0.62}]{Figures/fig_F_EPR_operator.png}}}\end{array}}\geq\frac{1}{d^{4}}\sum_{\textbf{s}_{1},\textbf{s}_{2}}\sum_{Q_{1},Q_{2}}{\begin{array}[]{c}\vbox{\vskip
1.80002pt\hbox{\hskip
1.20001pt\includegraphics[scale={0.6}]{Figures/fig_F_bound.png}}}\end{array}}.$
(183)

Namely,

$\displaystyle
F_{EPR}\geq\left|\sum_{\textbf{s},Q}\frac{1}{d^{2}}C_{\textbf{s},Q}\right|^{2},$
(184)

where the left-right correlation function $C_{\textbf{s},Q}$ is given by

$\displaystyle\begin{aligned} &C_{\textbf{s},Q}=\langle
TFD|\tilde{Q}_{R}^{\dagger}(t_{0})\left[\left(\otimes_{k}|s_{k}\rangle\langle
s_{k}|\right)\otimes U_{\textbf{s}}\right]Q_{L}(-t_{0})|TFD\rangle,\\\
&Q_{L}(-t_{0})=(U_{L}QU^{\dagger}_{L})\otimes\mathds{1},\quad\tilde{Q}_{R}(t_{0})=\mathds{1}\otimes(U_{R}^{\dagger}\tilde{Q}U_{R}),\end{aligned}$
(185)

where $U_{L,R}$ denotes the time evolution operator in the left and right
side. Thus, the left-right correlation function is closely related to the
fidelity bound.

##  Appendix G Symmetries of the twisted correlation function

In this appendix, we discuss the symmetries of the twisted correlation
function, which are used to simplify the calculations in the main text by
reducing to a fundamental domain. Recall the twisted correlation function is
defined as

$\displaystyle G_{\chi}(s_{1},s_{2})=$ (186) $\displaystyle\frac{\langle
L_{\textbf{l}}|e^{-\beta
H}\mathcal{T}\left[O(\frac{\beta}{2})\left(\begin{array}[]{cccc}\psi_{2k-1}(s_{1})\psi_{2k-1}(s_{2})&il_{k}\psi_{2k-1}(s_{1})\psi_{2k}(2\beta-
s_{2})\\\ il_{k}\psi_{2k}(2\beta-s_{1})\psi_{2k-1}(s_{2})&-\psi_{2k}(2\beta-
s_{1})\psi_{2k}(2\beta-
s_{2})\end{array}\right)\right]|L_{\textbf{l}}\rangle}{\langle
L_{\textbf{l}}|e^{-\beta
H}\mathcal{T}\left[O(\frac{\beta}{2})\right]|L_{\textbf{l}}\rangle}$ (189)

where $\psi(\tau)=e^{H\tau}\psi e^{-H\tau}$ is the imaginary time evolution,
and
$O(\frac{\beta}{2})=e^{-i\frac{\theta}{q}\sum_{k}l_{k}i\psi_{2k-1}(\frac{\beta}{2})\psi_{2k}(\frac{\beta}{2})}$
is the decoding operator. Because the Majorana field $\chi_{k}$ satisfies
conventional anti-periodic boundary conditions, we have

$\displaystyle G_{\chi}(s_{1},s_{2})=-G_{\chi}(s_{2},s_{1}),\quad
G_{\chi}(s_{1},s_{2})=-G_{\chi}(2\beta+s_{1},s_{2}).$ (190)

The twisted correlation function is real if $\theta$ is a purely imaginary
number, $G_{\chi}(s_{1},s_{2})=G_{\chi}(s_{1},s_{2})^{*}$. This is because we
can find a real representation for all operators, including
$O(\frac{\beta}{2})$. Thus, we take $\theta$ imaginary, then analytically
continue it to a real number. If we take the complex conjugate of the twisted
correlation function, we get

$\displaystyle\begin{aligned} &G_{\chi}(s_{1},s_{2})^{*}\\\ &=\frac{\langle
L_{\textbf{l}}|\mathcal{T}\left[O(-\frac{\beta}{2})\left(\begin{array}[]{cccc}\psi_{2k-1}(-s_{2})\psi_{2k-1}(-s_{1})&-il_{k}\psi_{2k}(-2\beta+s_{2})\psi_{2k-1}(-s_{1})\\\
-il_{k}\psi_{2k-1}(-s_{2})\psi_{2k}(-2\beta+s_{1})&-\psi_{2k}(-2\beta+s_{2})\psi_{2k}(-2\beta+s_{1})\end{array}\right)\right]e^{-\beta
H}|L_{\textbf{l}}\rangle}{\langle
L_{\textbf{l}}|\mathcal{T}\left[O(-\frac{\beta}{2})\right]e^{-\beta
H}|L_{\textbf{l}}\rangle}\\\ &=\frac{\langle L_{\textbf{l}}|e^{-\beta
H}\mathcal{T}\left[O(\frac{\beta}{2})\left(\begin{array}[]{cccc}\psi_{2k-1}(\beta-
s_{2})\psi_{2k-1}(\beta-
s_{1})&-il_{k}\psi_{2k}(-\beta+s_{2})\psi_{2k-1}(\beta-s_{1})\\\
-il_{k}\psi_{2k-1}(\beta-
s_{2})\psi_{2k}(-\beta+s_{1})&-\psi_{2k}(-\beta+s_{2})\psi_{2k}(-\beta+s_{1})\end{array}\right)\right]|L_{\textbf{l}}\rangle}{\langle
L_{\textbf{l}}|e^{-\beta
H}\mathcal{T}\left[O(\frac{\beta}{2})\right]|L_{\textbf{l}}\rangle}\\\
&=G_{\chi}(\beta-s_{2},\beta-s_{1}),\end{aligned}$ (191)

where we have used $[\psi(\tau)]^{\dagger}=\psi(-\tau)$ to get the first line.
In the second line, we bring the operator $e^{-\beta H}$ from the right to the
left by shifting the time by $\beta$ for all operators. The last line can be
obtained by recalling the definition of $\chi_{k}$ and its anti-periodic
boundary conditions.

Moreover, the twisted correlation function is invariant under
$\psi_{2k-1}\rightarrow\psi_{2k}$, $\psi_{2k}\rightarrow-\psi_{2k-1}$. The
large-$N$ action is invariant under $O(N)$ transformations, so it is invariant
under this subgroup. Now we should check the measurement-decoding channel
under this transformation. Both the projected state $|L_{\textbf{l}}\rangle$
and the decoding operator $O(\frac{\beta}{2})$ are invariant because
$-2i\psi_{2k-1}\psi_{2k}$ is invariant. Therefore, the twisted correlation
function is invariant, which leads to the following properties

$\displaystyle G_{\chi}(s_{1},s_{2})=-G_{\chi}(2\beta-s_{1},2\beta-
s_{2})=G_{\chi}(2\beta-s_{2},2\beta-s_{1}).$ (192)

In summary, (190), (191), and (192) lead to the following (anti-)reflection
conditions

$\displaystyle\begin{aligned} &\text{anti-reflection about}\quad
s_{1}-s_{2}=0:\quad G_{\chi}(s_{1},s_{2})=-G_{\chi}(s_{2},s_{1}),\\\
&\text{reflection about}\quad s_{1}-s_{2}=\beta:\quad
G_{\chi}(s_{1},s_{2})=G_{\chi}(s_{2},s_{1}-\beta),\\\ &\text{reflection
about}\quad s_{1}-s_{2}=-\beta:\quad
G_{\chi}(s_{1},s_{2})=G_{\chi}(s_{2}-\beta,s_{1}+\beta),\\\ &\text{reflection
about}\quad s_{1}+s_{2}=2\beta:\quad G_{\chi}(s_{1},s_{2})=G_{\chi}(2\beta-
s_{2},2\beta-s_{1}),\\\ &\text{reflection about}\quad s_{1}+s_{2}=2\beta:\quad
G_{\chi}(s_{1},s_{2})=G_{\chi}(\beta-s_{2},\beta-s_{1}),\\\ &\text{reflection
about}\quad s_{1}+s_{2}=\frac{3}{2}\beta:\quad
G_{\chi}(s_{1},s_{2})=G_{\chi}(\frac{3}{2}\beta-s_{2},\frac{3}{2}\beta-
s_{1}).\end{aligned}$ (193)

They are illustrated in Fig. 14. With the help of these symmetries, we can
reduce the calculation from the full domain to the fundamental domain $A$,
$B$, $C$ and $D$, as discussed in the main text.

## References

  * (1) J. M. Maldacena, _The Large N limit of superconformal field theories and supergravity_ , _Adv. Theor. Math. Phys._ 2 (1998) 231–252, [hep-th/9711200].
  * (2) E. Witten, _Anti-de Sitter space and holography_ , _Adv. Theor. Math. Phys._ 2 (1998) 253–291, [hep-th/9802150].
  * (3) S. S. Gubser, I. R. Klebanov and A. M. Polyakov, _Gauge theory correlators from noncritical string theory_ , _Phys. Lett. B_ 428 (1998) 105–114, [hep-th/9802109].
  * (4) O. Aharony, S. S. Gubser, J. M. Maldacena, H. Ooguri and Y. Oz, _Large N field theories, string theory and gravity_ , _Phys. Rept._ 323 (2000) 183–386, [hep-th/9905111].
  * (5) S. Ryu and T. Takayanagi, _Aspects of Holographic Entanglement Entropy_ , _JHEP_ 08 (2006) 045, [hep-th/0605073].
  * (6) S. Ryu and T. Takayanagi, _Holographic derivation of entanglement entropy from the anti–de sitter space/conformal field theory correspondence_ , _Physical review letters_ 96 (2006) 181602.
  * (7) V. E. Hubeny, M. Rangamani and T. Takayanagi, _A Covariant holographic entanglement entropy proposal_ , _JHEP_ 07 (2007) 062, [0705.0016].
  * (8) B. Swingle, _Entanglement Renormalization and Holography_ , _Phys. Rev. D_ 86 (2012) 065007, [0905.1317].
  * (9) M. Van Raamsdonk, _Building up spacetime with quantum entanglement_ , _Gen. Rel. Grav._ 42 (2010) 2323–2329, [1005.3035].
  * (10) J. Maldacena and L. Susskind, _Cool horizons for entangled black holes_ , _Fortsch. Phys._ 61 (2013) 781–811, [1306.0533].
  * (11) N. Engelhardt and A. C. Wall, _Quantum Extremal Surfaces: Holographic Entanglement Entropy beyond the Classical Regime_ , _JHEP_ 01 (2015) 073, [1408.3203].
  * (12) X. Dong, D. Harlow and A. C. Wall, _Reconstruction of Bulk Operators within the Entanglement Wedge in Gauge-Gravity Duality_ , _Phys. Rev. Lett._ 117 (2016) 021601, [1601.05416].
  * (13) D. Harlow, _The Ryu–Takayanagi Formula from Quantum Error Correction_ , _Commun. Math. Phys._ 354 (2017) 865–912, [1607.03901].
  * (14) S. Antonini, G. Bentsen, C. Cao, J. Harper, S.-K. Jian and B. Swingle, _Holographic measurement and bulk teleportation_ , 2209.12903.
  * (15) T. Numasawa, N. Shiba, T. Takayanagi and K. Watanabe, _Epr pairs, local projections and quantum teleportation in holography_ , _Journal of High Energy Physics_ 2016 (2016) 1–52.
  * (16) S. Sachdev and J. Ye, _Gapless spin fluid ground state in a random, quantum heisenberg magnet_ , _Phys. Rev. Lett._ 70 (1993) 3339, [cond-mat/9212030].
  * (17) A. Kitaev, “A simple model of quantum holography.” http://online.kitp.ucsb.edu/online/entangled15/, Feb. 12, April 7, and May 27, 2015.
  * (18) J. Maldacena and D. Stanford, _Remarks on the sachdev-ye-kitaev model_ , _Physical Review D_ 94 (2016) 106002.
  * (19) J. Maldacena, D. Stanford and Z. Yang, _Conformal symmetry and its breaking in two-dimensional nearly anti-de sitter space_ , _Progress of Theoretical and Experimental Physics_ 2016 (2016) .
  * (20) J. Maldacena, _Eternal black holes in anti-de sitter_ , _Journal of High Energy Physics_ 2003 (2003) 021.
  * (21) J. Maldacena and X.-L. Qi, _Eternal traversable wormhole_ , 1804.00491.
  * (22) V. P. Su, _Variational preparation of the thermofield double state of the Sachdev-Ye-Kitaev model_ , _Phys. Rev. A_ 104 (2021) 012427, [2009.04488].
  * (23) C. Teitelboim, _Gravitation and Hamiltonian Structure in Two Space-Time Dimensions_ , _Phys. Lett. B_ 126 (1983) 41–45.
  * (24) R. Jackiw, _Lower Dimensional Gravity_ , _Nucl. Phys. B_ 252 (1985) 343–356.
  * (25) I. Kourkoulou and J. Maldacena, _Pure states in the syk model and nearly- $ads\\_2$ gravity_, _arXiv preprint arXiv:1707.02325_ (2017) .
  * (26) S. Antonini and B. Swingle, _Holographic boundary states and dimensionally reduced braneworld spacetimes_ , _Phys. Rev. D_ 104 (2021) 046023, [2105.02912].
  * (27) J. L. Cardy, _Boundary conditions, fusion rules and the verlinde formula_ , _Nuclear Physics B_ 324 (1989) 581–596.
  * (28) S. Nezami, H. W. Lin, A. R. Brown, H. Gharibyan, S. Leichenauer, G. Salton et al., _Quantum Gravity in the Lab: Teleportation by Size and Traversable Wormholes, Part II_ , 2102.01064.
  * (29) V. Chandrasekaran and A. Levine, _Quantum error correction in syk and bulk emergence_ , _arXiv preprint arXiv:2203.05058_ (2022) .
  * (30) S. Antonini, B. Grado-White, S.-K. Jian and B. Swingle, _Quantum extremal surface formula in the syk model_ , _In preparation_ .
  * (31) P. Gao, D. L. Jafferis and A. C. Wall, _Traversable wormholes via a double trace deformation_ , _Journal of High Energy Physics_ 2017 (2017) 1–25.
  * (32) J. Maldacena, D. Stanford and Z. Yang, _Diving into traversable wormholes_ , _Fortschritte der Physik_ 65 (2017) 1700034.
  * (33) A. R. Brown, H. Gharibyan, S. Leichenauer, H. W. Lin, S. Nezami, G. Salton et al., _Quantum gravity in the lab: teleportation by size and traversable wormholes_ , _arXiv preprint arXiv:1911.06314_ (2019) .
  * (34) P. Gao and D. L. Jafferis, _A traversable wormhole teleportation protocol in the syk model_ , _Journal of High Energy Physics_ 2021 (2021) 1–44.
  * (35) P. Zhang, _Entanglement entropy and its quench dynamics for pure states of the sachdev-ye-kitaev model_ , _Journal of High Energy Physics_ 2020 (2020) 1–20.
  * (36) S. Sachdev, _Bekenstein-Hawking Entropy and Strange Metals_ , _Phys. Rev. X_ 5 (2015) 041025, [1506.05111].
  * (37) C. Liu, X. Chen and L. Balents, _Quantum entanglement of the sachdev-ye-kitaev models_ , _Physical Review B_ 97 (2018) 245126.
  * (38) P. Zhang, C. Liu and X. Chen, _Subsystem rényi entropy of thermal ensembles for syk-like models_ , _SciPost Physics_ 8 (2020) 094.
  * (39) I. Affleck and A. W. Ludwig, _Universal noninteger “ground-state degeneracy”in critical quantum systems_ , _Physical Review Letters_ 67 (1991) 161.
  * (40) A. Hamilton, D. Kabat, G. Lifschytz and D. A. Lowe, _Holographic representation of local bulk operators_ , _Physical Review D_ 74 (2006) 066009.
  * (41) Y. Huang, Y. Gu et al., _Eigenstate entanglement in the sachdev-ye-kitaev model_ , _Physical Review D_ 100 (2019) 041901.
  * (42) A. M. García-García and J. J. Verbaarschot, _Analytical spectral density of the sachdev-ye-kitaev model at finite n_ , _Physical Review D_ 96 (2017) 066012.
  * (43) T. Schuster, B. Kobrin, P. Gao, I. Cong, E. T. Khabiboulline, N. M. Linke et al., _Many-body quantum teleportation via operator spreading in the traversable wormhole protocol_ , _Physical Review X_ 12 (2022) 031013\.
  * (44) A. Milekhin and F. K. Popov, _Measurement-induced phase transition in teleportation and wormholes_ , 2210.03083.
  * (45) X.-L. Qi and A. Streicher, _Quantum epidemiology: operator growth, thermal effects, and syk_ , _Journal of High Energy Physics_ 2019 (2019) 1–28.
  * (46) B. Yoshida and N. Y. Yao, _Disentangling scrambling and decoherence via quantum teleportation_ , _Physical Review X_ 9 (2019) 011006.
  * (47) M. Miyaji, S. Ryu, T. Takayanagi and X. Wen, _Boundary states as holographic duals of trivial spacetimes_ , _Journal of High Energy Physics_ 2015 (May, 2015) 152, [1412.6226].
  * (48) T. Takayanagi, _Holographic dual of a boundary conformal field theory_ , _Physical review letters_ 107 (2011) 101602.
  * (49) M. Fujita, T. Takayanagi and E. Tonni, _Aspects of ads/bcft_ , _Journal of High Energy Physics_ 2011 (2011) 1–40.
  * (50) G. Penington, _Entanglement wedge reconstruction and the information paradox_ , _Journal of High Energy Physics_ 2020 (2020) 1–84.
  * (51) G. Penington, S. H. Shenker, D. Stanford and Z. Yang, _Replica wormholes and the black hole interior_ , _Journal of High Energy Physics_ 2022 (2022) 1–87.
  * (52) A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian and A. Tajdini, _Replica wormholes and the entropy of hawking radiation_ , _Journal of High Energy Physics_ 2020 (may, 2020) .
  * (53) G. T. Horowitz and J. M. Maldacena, _The Black hole final state_ , _JHEP_ 02 (2004) 008, [hep-th/0310281].
  * (54) C. Akers and G. Penington, _Quantum minimal surfaces from quantum error correction_ , _SciPost Physics_ 12 (2022) 157.
  * (55) C. Akers, N. Engelhardt, D. Harlow, G. Penington and S. Vardhan, _The black hole interior from non-isometric codes and complexity_ , _arXiv preprint arXiv:2207.06536_ (2022) .
  * (56) Y. Li, X. Chen and M. P. Fisher, _Measurement-driven entanglement transition in hybrid quantum circuits_ , _Physical Review B_ 100 (2019) 134306.
  * (57) B. Skinner, J. Ruhman and A. Nahum, _Measurement-induced phase transitions in the dynamics of entanglement_ , _Physical Review X_ 9 (2019) 031009\.
  * (58) A. Chan, R. M. Nandkishore, M. Pretko and G. Smith, _Unitary-projective entanglement dynamics_ , _Physical Review B_ 99 (2019) 224307.
  * (59) S.-K. Jian, C. Liu, X. Chen, B. Swingle and P. Zhang, _Measurement-induced phase transition in the monitored sachdev-ye-kitaev model_ , _Phys. Rev. Lett._ 127 (Sep, 2021) 140601.
  * (60) G. S. Bentsen, S. Sahu and B. Swingle, _Measurement-induced purification in large-N hybrid Brownian circuits_ , _Phys. Rev. B_ 104 (2021) 094304, [2104.07688].
  * (61) K. Kawabata, A. Kulkarni, J. Li, T. Numasawa and S. Ryu, _Dynamical quantum phase transitions in syk lindbladians_ , _arXiv preprint arXiv:2210.04093_ (2022) .
  * (62) A. M. García-García, L. Sá, J. J. Verbaarschot and J. P. Zheng, _Keldysh wormholes and anomalous relaxation in the dissipative sachdev-ye-kitaev model_ , _arXiv preprint arXiv:2210.01695_ (2022) .
  * (63) K. Goto, M. Nozaki, K. Tamaoka and M. T. Tan, _Entanglement dynamics of the non-unitary holographic channel_ , _arXiv preprint arXiv:2211.03944_ (2022) .